C     MAIN PROGRAM
      INTEGER LUNIT
C     ALLOW 5000 UNDERFLOWS.
      CALL TRAPS(0,0,5001,0,0)
C
C     OUTPUT UNIT NUMBER
C
      LUNIT = 6
C
      CALL DCHTS(LUNIT)
      STOP
      END
      SUBROUTINE DCHTS(LUNIT)
C     LUNIT IS THE OUTPUT UNIT NUMBER
C
C     TESTS
C        DCHDC
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     J.J. DONGARRA AND G.W. STEWART, ARGONNE NATIONAL LABORATORY,
C     UNIVERSITY OF MARYLAND
C
C     SUBROUTINES AND FUNCTIONS
C
C     LINPACK DCHDC
C     EXTERNAL DAGEN,DCHXX,DMACH,DARRAY
C     FORTRAN DMAX1,DABS,DBLE,FLOAT
C
C     INTERNAL VARIABLES
C
      INTEGER LUNIT
      INTEGER CASE,LDA,P,JPVT(25),KPVT(25),KD
      DOUBLE PRECISION RESDUL
      DOUBLE PRECISION A(25,25),AA(25,25),ASAVE(25,25),D(25),WORK(25)
      DOUBLE PRECISION DMACH,TINY,HUGE,EPS
      LOGICAL NOTWRT
      LDA = 25
      NCASE = 7
      ERRLVL = 100.0D0
      NOTWRT = .TRUE.
      EPS = DMACH(1)
      TINY = DMACH(2)*10000.0D0
      HUGE = DMACH(3)
      WRITE (LUNIT,300)
      DO 290 CASE = 1, NCASE
         WRITE (LUNIT,310) CASE
         GO TO (10,50,80,120,160,200,240), CASE
   10    CONTINUE
C
C           5 X 5 CASE NO PIVOTING.
C
            WRITE (LUNIT,320)
            P = 5
C
C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
C
            DO 20 I = 1, P
               D(I) = DBLE(FLOAT(I))
               JPVT(I) = 0
   20       CONTINUE
            CALL DAGEN(A,LDA,P,D,ASAVE)
            IF (NOTWRT) GO TO 30
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
   30       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 0
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 40
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
   40       CONTINUE
         GO TO 280
   50    CONTINUE
C
C           1 X 1 CASE WITH PIVOTING.
C
            WRITE (LUNIT,330)
            P = 1
C
C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
C
            A(1,1) = 1.0D0
            ASAVE(1,1) = A(1,1)
            JPVT(1) = 0
            IF (NOTWRT) GO TO 60
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
   60       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 1
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 70
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
   70       CONTINUE
         GO TO 280
   80    CONTINUE
C
C           8 X 8 CASE PIVOT LOGIC TEST
C
            WRITE (LUNIT,340)
            P = 8
C
C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
C
            DO 90 I = 1, P
               D(I) = DBLE(FLOAT(I))
               JPVT(I) = 1 - 2*MOD(I+1,2)
   90       CONTINUE
            CALL DAGEN(A,LDA,P,D,ASAVE)
            IF (NOTWRT) GO TO 100
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  100       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 1
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 110
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  110       CONTINUE
         GO TO 280
  120    CONTINUE
C
C           6 X 6 CASE PIVOTING WITH NEGATIVE EIGENVALUE.
C
            WRITE (LUNIT,350)
            P = 6
C
C           GENERATE MATRIX
C
            D(1) = 1.0D0
            D(2) = 1.0D0/2.0D0
            D(3) = 1.0D0/4.0D0
            D(4) = -1.0D0
            D(5) = 1.0D0/2.0D0
            D(6) = 1.0D0/4.0D0
            DO 130 I = 1, P
               JPVT(I) = 0
  130       CONTINUE
            CALL DAGEN(A,LDA,P,D,ASAVE)
            IF (NOTWRT) GO TO 140
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  140       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 1
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 150
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  150       CONTINUE
         GO TO 280
  160    CONTINUE
C
C           25 X 25 CASE WITH PIVOTING.
C
            WRITE (LUNIT,360)
            P = 25
C
C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
C
            DO 170 I = 1, P
               D(I) = DBLE(FLOAT(I))
               JPVT(I) = 0
  170       CONTINUE
            CALL DAGEN(A,LDA,P,D,ASAVE)
            IF (NOTWRT) GO TO 180
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  180       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 1
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 190
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  190       CONTINUE
         GO TO 280
  200    CONTINUE
C
C           5 X 5 CASE WITH PIVOTING AND UNDERFLOW CHECK.
C
            WRITE (LUNIT,370)
            P = 5
C
C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
C
            DO 210 I = 1, P
               D(I) = DBLE(FLOAT(I))*TINY
               JPVT(I) = 0
  210       CONTINUE
            CALL DAGEN(A,LDA,P,D,ASAVE)
            IF (NOTWRT) GO TO 220
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  220       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 1
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 230
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  230       CONTINUE
         GO TO 280
  240    CONTINUE
C
C           5 X 5 CASE WITH PIVOTING AND OVERFLOW CHECK.
C
            WRITE (LUNIT,380)
            P = 5
C
C           GENERATE MATRIX WITH POSITIVE EIGENVALUES.
C
            DO 250 I = 1, P
               D(I) = DBLE(FLOAT(I))*HUGE
               JPVT(I) = 0
  250       CONTINUE
            CALL DAGEN(A,LDA,P,D,ASAVE)
            IF (NOTWRT) GO TO 260
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  260       CONTINUE
C
C           DECOMPOSE THE MATRIX.
C
            JOB = 1
            WRITE (LUNIT,420) JOB,(JPVT(I), I = 1, P)
            CALL DCHDC(A,LDA,P,WORK,JPVT,JOB,KD)
            WRITE (LUNIT,430) KD
            CALL DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
            WRITE (LUNIT,400) (JPVT(I), I = 1, P)
            RESDUL = RESDUL/EPS
            WRITE (LUNIT,410) RESDUL
            IF (NOTWRT) GO TO 270
               WRITE (LUNIT,390)
               CALL DARRAY(A,LDA,P,P,P,LUNIT)
  270       CONTINUE
  280    CONTINUE
  290 CONTINUE
      WRITE (LUNIT,440)
      RETURN
  300 FORMAT (22H1LINPACK TESTER, DCH** /
     *        29H THIS VERSION DATED 08/14/78.)
  310 FORMAT ( // 5H1CASE, I3)
  320 FORMAT ( / 19H 5 X 5 NO PIVOTING.)
  330 FORMAT ( / 22H MONOELEMENTAL MATRIX.)
  340 FORMAT ( / 24H 8 X 8 PIVOT LOGIC TEST.)
  350 FORMAT ( / 32H 6 X 6 NEGATIVE EIGENVALUE TEST.)
  360 FORMAT ( / 16H 25 X 25 MATRIX.)
  370 FORMAT ( / 32H 5 X 5 PIVOT AND UNDERFLOW TEST.)
  380 FORMAT ( / 31H 5 X 5 PIVOT AND OVERFLOW TEST.)
  390 FORMAT ( / 2H A)
  400 FORMAT ( / 5H JPVT // (1H , 10I5))
  410 FORMAT ( // 18H A - TRANS(R)*R = , 1PD16.8)
  420 FORMAT ( / 39H JOB AND JPVT BEFORE THE DECOMPOSITION. / I5 /
     *         (10I5))
  430 FORMAT ( / 20H THE VALUE OF  KD  =, I5)
  440 FORMAT ( /// 12H END OF TEST / 1H1)
      END
      SUBROUTINE DCHXX(A,LDA,P,AA,JPVT,KPVT,ASAVE,RESDUL,JOB)
C
      INTEGER LDA,P,JPVT(1),KPVT(1),JOB
      DOUBLE PRECISION A(LDA,1),AA(LDA,1),ASAVE(LDA,1)
      DOUBLE PRECISION RESDUL,ANORM,DASUM
C
      DOUBLE PRECISION TEMP,DDOT
C
C     FORM TRANS(R)*R
C
      DO 20 I = 1, P
         DO 10 J = I, P
            AA(I,J) = DDOT(I,A(1,I),1,A(1,J),1)
            AA(J,I) = AA(I,J)
   10    CONTINUE
         KPVT(I) = JPVT(I)
   20 CONTINUE
      IF (JOB .EQ. 0) GO TO 70
C
C        UNSCRAMBLE TRANS(R)*R
C
         DO 60 J = 1, P
   30       IF (KPVT(J) .EQ. J) GO TO 50
               IK = KPVT(J)
               CALL DSWAP(P,AA(1,J),1,AA(1,IK),1)
               DO 40 I = 1, P
                  TEMP = AA(J,I)
                  AA(J,I) = AA(IK,I)
                  AA(IK,I) = TEMP
   40          CONTINUE
               KPVT(J) = KPVT(IK)
               KPVT(IK) = IK
            GO TO 30
   50       CONTINUE
   60    CONTINUE
   70 CONTINUE
      ANORM = 0.0D0
      RESDUL = 0.0D0
      DO 90 J = 1, P
         ANORM = DMAX1(ANORM,DASUM(P,ASAVE(1,J),1))
         DO 80 I = 1, P
            ASAVE(I,J) = ASAVE(I,J) - AA(I,J)
   80    CONTINUE
         RESDUL = DMAX1(RESDUL,DASUM(P,ASAVE(1,J),1))
   90 CONTINUE
      RESDUL = RESDUL/ANORM
      RETURN
      END
      SUBROUTINE DARRAY(A,LDA,M,N,NNL,LUNIT)
      INTEGER LDA,M,N,NNL,LUNIT
      DOUBLE PRECISION A(LDA,1)
C
C     FORTRAN IABS,MIN0
C
      INTEGER I,J,K,KU,NL
      NL = IABS(NNL)
      IF (NNL .LT. 0) GO TO 30
         DO 20 I = 1, M
            WRITE (LUNIT,70)
            DO 10 K = 1, N, NL
               KU = MIN0(K+NL-1,N)
               WRITE (LUNIT,70) (A(I,J), J = K, KU)
   10       CONTINUE
   20    CONTINUE
      GO TO 60
   30 CONTINUE
         DO 50 J = 1, N
            WRITE (LUNIT,70)
            DO 40 K = 1, M, NL
               KU = MIN0(K+NL-1,M)
               WRITE (LUNIT,70) (A(I,J), I = K, KU)
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
      RETURN
   70 FORMAT (1H , 8D13.6)
      END
      SUBROUTINE DAGEN(A,LDA,P,D,ASAVE)
      INTEGER LDA,P
      DOUBLE PRECISION A(LDA,1),D(1),ASAVE(LDA,1)
C
      DOUBLE PRECISION EN,ES,FDPTP,S,TDP
C
      S = 0.0D0
      DO 10 I = 1, P
         S = S + D(I)
   10 CONTINUE
      TDP = 2.0D0/DBLE(FLOAT(P))
      FDPTP = 4.0D0*S/DBLE(FLOAT(P*P))
      DO 30 I = 1, P
         ES = -1.0D0
         DO 20 J = I, P
            ES = (-1.0D0)*ES
            A(I,J) = -ES*TDP*(D(J) + D(I)) + ES*FDPTP
            A(I,J) = A(I,J)*DBLE(FLOAT(I+1))*DBLE(FLOAT(J+1))
     *               /DBLE(FLOAT((J+1)**2+J**2))
            ASAVE(I,J) = A(I,J)
            A(J,I) = A(I,J)
            ASAVE(J,I) = A(J,I)
   20    CONTINUE
         A(I,I) = A(I,I) + D(I)
         ASAVE(I,I) = A(I,I)
   30 CONTINUE
      RETURN
      END