C
C     ..................................................................
C
C        SUBROUTINE ATEIG
C
C        PURPOSE
C           COMPUTE THE EIGENVALUES OF A REAL ALMOST TRIANGULAR MATRIX
C
C        USAGE
C           CALL ATEIG(M,A,RR,RI,IANA,IA)
C
C        DESCRIPTION OF THE PARAMETERS
C           M      ORDER OF THE MATRIX
C           A      THE INPUT MATRIX, M BY M
C           RR     VECTOR CONTAINING THE REAL PARTS OF THE EIGENVALUES
C                  ON RETURN
C           RI     VECTOR CONTAINING THE IMAGINARY PARTS OF THE EIGEN-
C                  VALUES ON RETURN
C           IANA   VECTOR WHOSE DIMENSION MUST BE GREATER THAN OR EQUAL
C                  TO M, CONTAINING ON RETURN INDICATIONS ABOUT THE WAY
C                  THE EIGENVALUES APPEARED (SEE MATH. DESCRIPTION)
C           IA     SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A
C                  IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE
C                  SUBSCRIPTED DATA STORAGE MODE.
C                  IA=M WHEN THE MATRIX IS IN SSP VECTOR STORAGE MODE.
C
C        REMARKS
C           THE ORIGINAL MATRIX IS DESTROYED
C           THE DIMENSION OF RR AND RI MUST BE GREATER OR EQUAL TO M
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           QR DOUBLE ITERATION
C
C        REFERENCES
C           J.G.F. FRANCIS - THE QR TRANSFORMATION---THE COMPUTER
C           JOURNAL, VOL. 4, NO. 3, OCTOBER 1961, VOL. 4, NO. 4, JANUARY
C           1962.  J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM -
C           CLARENDON PRESS, OXFORD, 1965.
C
C     ..................................................................
C
      SUBROUTINE ATEIG(M,A,RR,RI,IANA,IA)
      DIMENSION A(1),RR(1),RI(1),PRR(2),PRI(2),IANA(1)
      INTEGER P,P1,Q
C
      E7=1.0E-8
      E6=1.0E-6
      E10=1.0E-10
      DELTA=0.5
      MAXIT=30
C
C        INITIALIZATION
C
      N=M
   20 N1=N-1
      IN=N1*IA
      NN=IN+N
      IF(N1) 30,1300,30
   30 NP=N+1
C
C        ITERATION COUNTER
C
      IT=0
C
C        ROOTS OF THE 2ND ORDER MAIN SUBMATRIX AT THE PREVIOUS
C        ITERATION
C
      DO 40 I=1,2
      PRR(I)=0.0
   40 PRI(I)=0.0
C
C        LAST TWO SUBDIAGONAL ELEMENTS AT THE PREVIOUS ITERATION
C
      PAN=0.0
      PAN1=0.0
C
C        ORIGIN SHIFT
C
      R=0.0
      S=0.0
C
C        ROOTS OF THE LOWER MAIN 2 BY 2 SUBMATRIX
C
      N2=N1-1
      IN1=IN-IA
      NN1=IN1+N
      N1N=IN+N1
      N1N1=IN1+N1
   60 T=A(N1N1)-A(NN)
      U=T*T
      V=4.0*A(N1N)*A(NN1)
      IF(ABS(V)-U*E7) 100,100,65
   65 T=U+V
      IF(ABS(T)-AMAX1(U,ABS(V))*E6) 67,67,68
   67 T=0.0
   68 U=(A(N1N1)+A(NN))/2.0
      V=SQRT(ABS(T))/2.0
      IF(T)140,70,70
   70 IF(U) 80,75,75
   75 RR(N1)=U+V
      RR(N)=U-V
      GO TO 130
   80 RR(N1)=U-V
      RR(N)=U+V
      GO TO 130
  100 IF(T)120,110,110
  110 RR(N1)=A(N1N1)
      RR(N)=A(NN)
      GO TO 130
  120 RR(N1)=A(NN)
      RR(N)=A(N1N1)
  130 RI(N)=0.0
      RI(N1)=0.0
      GO TO 160
  140 RR(N1)=U
      RR(N)=U
      RI(N1)=V
      RI(N)=-V
  160 IF(N2)1280,1280,180
C
C        TESTS OF CONVERGENCE
C
  180 N1N2=N1N1-IA
      RMOD=RR(N1)*RR(N1)+RI(N1)*RI(N1)
      EPS=E10*SQRT(RMOD)
      IF(ABS(A(N1N2))-EPS)1280,1280,240
  240 IF(ABS(A(NN1))-E10*ABS(A(NN))) 1300,1300,250
  250 IF(ABS(PAN1-A(N1N2))-ABS(A(N1N2))*E6) 1240,1240,260
  260 IF(ABS(PAN-A(NN1))-ABS(A(NN1))*E6)1240,1240,300
  300 IF(IT-MAXIT) 320,1240,1240
C
C        COMPUTE THE SHIFT
C
  320 J=1
      DO 360 I=1,2
      K=NP-I
      IF(ABS(RR(K)-PRR(I))+ABS(RI(K)-PRI(I))-DELTA*(ABS(RR(K))
     1    +ABS(RI(K)))) 340,360,360
  340 J=J+I
  360 CONTINUE
      GO TO (440,460,460,480),J
  440 R=0.0
      S=0.0
      GO TO 500
  460 J=N+2-J
      R=RR(J)*RR(J)
      S=RR(J)+RR(J)
      GO TO 500
  480 R=RR(N)*RR(N1)-RI(N)*RI(N1)
      S=RR(N)+RR(N1)
C
C        SAVE THE LAST TWO SUBDIAGONAL TERMS AND THE ROOTS OF THE
C        SUBMATRIX BEFORE ITERATION
C
  500 PAN=A(NN1)
      PAN1=A(N1N2)
      DO 520 I=1,2
      K=NP-I
      PRR(I)=RR(K)
  520 PRI(I)=RI(K)
C
C        SEARCH FOR A PARTITION OF THE MATRIX, DEFINED BY P AND Q
C
      P=N2
      IF (N-3)600,600,525
  525 IPI=N1N2
      DO 580 J=2,N2
      IPI=IPI-IA-1
      IF(ABS(A(IPI))-EPS) 600,600,530
  530 IPIP=IPI+IA
      IPIP2=IPIP+IA
      D=A(IPIP)*(A(IPIP)-S)+A(IPIP2)*A(IPIP+1)+R
      IF(D)540,560,540
  540 IF(ABS(A(IPI)*A(IPIP+1))*(ABS(A(IPIP)+A(IPIP2+1)-S)+ABS(A(IPIP2+2)
     1 )) -ABS(D)*EPS) 620,620,560
  560 P=N1-J
  580 CONTINUE
  600 Q=P
      GO TO 680
  620 P1=P-1
      Q=P1
      IF (P1-1) 680,680,650
  650 DO 660 I=2, P1
      IPI=IPI-IA-1
      IF(ABS(A(IPI))-EPS)680,680,660
  660 Q=Q-1
C
C        QR DOUBLE ITERATION
C
  680 II=(P-1)*IA+P
      DO 1220 I=P,N1
      II1=II-IA
      IIP=II+IA
      IF(I-P)720,700,720
  700 IPI=II+1
      IPIP=IIP+1
C
C        INITIALIZATION OF THE TRANSFORMATION
C
      G1=A(II)*(A(II)-S)+A(IIP)*A(IPI)+R
      G2=A(IPI)*(A(IPIP)+A(II)-S)
      G3=A(IPI)*A(IPIP+1)
      A(IPI+1)=0.0
      GO TO 780
  720 G1=A(II1)
      G2=A(II1+1)
      IF(I-N2)740,740,760
  740 G3=A(II1+2)
      GO TO 780
  760 G3=0.0
  780 CAP=SQRT(G1*G1+G2*G2+G3*G3)
      IF(CAP)800,860,800
  800 IF(G1)820,840,840
  820 CAP=-CAP
  840 T=G1+CAP
      PSI1=G2/T
      PSI2=G3/T
      ALPHA=2.0/(1.0+PSI1*PSI1+PSI2*PSI2)
      GO TO 880
  860 ALPHA=2.0
      PSI1=0.0
      PSI2=0.0
  880 IF(I-Q)900,960,900
  900 IF(I-P)920,940,920
  920 A(II1)=-CAP
      GO TO 960
  940 A(II1)=-A(II1)
C
C        ROW OPERATION
C
  960 IJ=II
      DO 1040 J=I,N
      T=PSI1*A(IJ+1)
      IF(I-N1)980,1000,1000
  980 IP2J=IJ+2
      T=T+PSI2*A(IP2J)
 1000 ETA=ALPHA*(T+A(IJ))
      A(IJ)=A(IJ)-ETA
      A(IJ+1)=A(IJ+1)-PSI1*ETA
      IF(I-N1)1020,1040,1040
 1020 A(IP2J)=A(IP2J)-PSI2*ETA
 1040 IJ=IJ+IA
C
C        COLUMN OPERATION
C
      IF(I-N1)1080,1060,1060
 1060 K=N
      GO TO 1100
 1080 K=I+2
 1100 IP=IIP-I
      DO 1180 J=Q,K
      JIP=IP+J
      JI=JIP-IA
      T=PSI1*A(JIP)
      IF(I-N1)1120,1140,1140
 1120 JIP2=JIP+IA
      T=T+PSI2*A(JIP2)
 1140 ETA=ALPHA*(T+A(JI))
      A(JI)=A(JI)-ETA
      A(JIP)=A(JIP)-ETA*PSI1
      IF(I-N1)1160,1180,1180
 1160 A(JIP2)=A(JIP2)-ETA*PSI2
 1180 CONTINUE
      IF(I-N2)1200,1220,1220
 1200 JI=II+3
      JIP=JI+IA
      JIP2=JIP+IA
      ETA=ALPHA*PSI2*A(JIP2)
      A(JI)=-ETA
      A(JIP)=-ETA*PSI1
      A(JIP2)=A(JIP2)-ETA*PSI2
 1220 II=IIP+1
      IT=IT+1
      GO TO 60
C
C        END OF ITERATION
C
 1240 IF(ABS(A(NN1))-ABS(A(N1N2))) 1300,1280,1280
C
C        TWO EIGENVALUES HAVE BEEN FOUND
C
 1280 IANA(N)=0
      IANA(N1)=2
      N=N2
      IF(N2)1400,1400,20
C
C        ONE EIGENVALUE HAS BEEN FOUND
C
 1300 RR(N)=A(NN)
      RI(N)=0.0
      IANA(N)=1
      IF(N1)1400,1400,1320
 1320 N=N1
      GO TO 20
 1400 RETURN
      END
