C
C
C
C
      SUBROUTINE DIAG(N,H,U)
C        N = ORDER OF H
C        H = MATRIX TO BE DIAGONALIZED
C        H IN OUTPUT CONTAINS THE EIGENVALUES ON DIAGONAL
C        U CONTAINS THE EIGENVECTORS (?)
C        IEGEN is a control parameter whose value is zero
C             or non-zero. Here IEGEN = 0
C        MAXIMUM ORDER = 20
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION H(42,42),U(42,42),X(42),IQ(42)
      IEGEN=0
      IF(IEGEN)15,10,15
   10 DO 14 I=1,N
      DO 14 J=1,N
      IF(I-J)12,11,12
   11 U(I,J)=1.
      GO TO 14
   12 U(I,J)=0.
   14 CONTINUE
   15 IF(N-1)1000,1000,17
   17 NMI1=N-1
      DO 30 I=1,NMI1
      X(I)=0.
      IPL1=I+1
      DO 30 J=IPL1,N
      IF(X(I)-DABS(H(I,J)))20,20,30
   20 X(I)=DABS(H(I,J))
      IQ(I)=J
   30 CONTINUE
      RAP=.745058059D-8
      HDTEST=1.D38
   40 DO 70 I=1,NMI1
      IF(I-1)60,60,45
   45 continue
      critic=1.d-10
      abxmax=dabs(xmax)
      if(abxmax.ge.critic) goto 701
      abx=dabs(x(i))
      if(abx.ge.critic) goto 701 
      goto 70
  701 IF(XMAX-X(I))60,70,70
   60 XMAX=X(I)
      IPIV=I
      JPIV=IQ(I)
   70 CONTINUE
      abxmax=dabs(xmax)
      if(abxmax.ge.critic) goto 702
      goto 1000
  702 IF(XMAX)1000,1000,80
   80 continue
      abhdte=dabs(hdtest)
      if(abhdte.ge.critic) goto 703
      goto 90
  703 IF(HDTEST)90,90,85
   85 continue
      abxmax=dabs(xmax)
      if(abxmax.ge.critic) goto 704
      abhdte=dabs(hdtest)
      if(abhdte.ge.critic) goto 704
      goto 90
  704 IF(XMAX - HDTEST)90,90,148
   90 HDMIN=DABS(H(1,1))
      DO 110 I=2,N
      abhdim=dabs(hdimin)
      if(abhdim.ge.critic) goto 705
      if(dabs(h(i,i)).ge.critic) goto 705
      goto 110
  705 IF(HDIMIN-DABS(H(I,I)))110,110,100
  100 HDIMIN=DABS(H(I,I))
  110 CONTINUE
      HDTEST = HDIMIN*RAP
      IF(HDTEST-XMAX) 148,1000,1000
  148 XDIF=H(IPIV,IPIV)-H(JPIV,JPIV)
C
C        AMENDMENT !!!
C
      IF(XDIF.GT.1.D-20) GO TO 9999
      IF(XDIF.LT.-1.D-20) GO TO 9999
      IF(XDIF)9997,9998,9998
 9997 XDIF=-1.D-20
      GO TO 9999
 9998 XDIF=1.D-20
C
C        AMENDMENT !!!
C
 9999 XO=DSIGN(2.D0,XDIF)*H(IPIV,JPIV)
      XS=XDIF**2+4.D0*H(IPIV,JPIV)**2
  150 TANG=XO/(DABS(XDIF)+DSQRT(XS))
      COSINE=1.D0/DSQRT(1.D0+TANG**2)
      SINE=TANG*COSINE
      HII=H(IPIV,IPIV)
      H(IPIV,IPIV)= COSINE**2*(HII+TANG*(2.D0*H(IPIV,JPIV)+TANG*H(JPIV,
     *    JPIV)))
      H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.D0*H(IPIV,JPIV)-TANG*
     *   HII))
      H(IPIV,JPIV)=0.D0
      IF(H(IPIV,IPIV) - H(JPIV,JPIV))152,153,153
  152 HTEMP = H(IPIV,IPIV)
      H(IPIV,IPIV) = H(JPIV,JPIV)
      H(JPIV,JPIV) =HTEMP
      HTEMP=DSIGN(1.D0,-SINE)*COSINE
      COSINE=DABS(SINE)
      SINE =HTEMP
  153 CONTINUE
      DO 350 I=1,NMI1
      IF(I-IPIV)210,350,200
  200 IF(I-JPIV)210,350,210
  210 IF(IQ(I)-IPIV)230,240,230
  230 IF(IQ(I)-JPIV)350,240,350
  240 K=IQ(I)
  250 HTEMP=H(I,K)
      H(I,K)=0.
      IPL1=I+1
      X(I)=0.
      DO 320 J=IPL1,N
      IF(X(1)-DABS(H(I,J)))300,300,320
  300 X(I)=DABS(H(I,J))
      IQ(I)=J
  320 CONTINUE
      H(I,K)=HTEMP
  350 CONTINUE
      X(IPIV)=0.
      X(JPIV)=0.
      DO 530 I=1,N
      IF(I-IPIV)370,530,420
  370 HTEMP = H(I,IPIV)
      H(I,IPIV)= COSINE*HTEMP +SINE*H(I,JPIV)
      IF(X(I)-DABS(H(I,IPIV)))380,390,390
  380 X(I)= DABS(H(I,IPIV))
      IQ(I) = IPIV
  390 H(I,JPIV) = -SINE*HTEMP + COSINE*H(I,JPIV)
      IF(X(I)-DABS(H(I,JPIV)))400,530,530
  400 X(I)=DABS(H(I,JPIV))
      IQ(I) = JPIV
      GO TO 530
  420 IF(I-JPIV) 430,530,480
  430 HTEMP = H(IPIV,I)
      H(IPIV,I) = COSINE*HTEMP +  SINE*H(I,JPIV)
      IF(X(IPIV)-DABS(H(IPIV,I)))440,450,450
  440 X(IPIV)=DABS(H(IPIV,I))
      IQ(IPIV) = I
  450 H(I,JPIV) = - SINE*HTEMP + COSINE*H(I,JPIV)
      IF(X(I)-DABS(H(I,JPIV)))400,530,530
  480 HTEMP = H(IPIV,I)
      H(IPIV,I) = COSINE*HTEMP + SINE*H(JPIV,I)
      IF(X(IPIV)-DABS(H(IPIV,I)))490,500,500
  490 X(IPIV)=DABS(H(IPIV,I))
      IQ(IPIV) = I
  500 H(JPIV,I) = -SINE*HTEMP + COSINE*H(JPIV,I)
      IF(X(IPIV)-DABS(H(JPIV,I)))510,530,530
  510 X(JPIV)=DABS(H(JPIV,I))
      IQ(JPIV) = I
  530 CONTINUE
      IF(IEGEN) 40,540,40
  540 DO 550 I=1,N
      HTEMP=U(I,IPIV)
      U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV)
  550 U(I,JPIV)= -SINE*HTEMP+COSINE*U(I,JPIV)
      GO TO 40
 1000 CONTINUE
      I=0
  560 I=I+1
      IF(I.EQ.N) GO TO 580
      IF(H(I,I).GE.H(I+1,I+1)) GO TO 560
      HH=H(I,I)
      H(I,I)=H(I+1,I+1)
      H(I+1,I+1)=HH
  570 IF(I.EQ.1) GO TO 560
      IF(H(I-1,I-1).GE.H(I,I)) GO TO 560
      HH=H(I-1,I-1)
      H(I-1,I-1)=H(I,I)
      H(I,I)=HH
      I=I-1
      GO TO 570
  580 RETURN
      END

