home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / mateigen / ateig.for next >
Text File  |  1985-11-29  |  7KB  |  300 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE ATEIG
  5. C
  6. C        PURPOSE
  7. C           COMPUTE THE EIGENVALUES OF A REAL ALMOST TRIANGULAR MATRIX
  8. C
  9. C        USAGE
  10. C           CALL ATEIG(M,A,RR,RI,IANA,IA)
  11. C
  12. C        DESCRIPTION OF THE PARAMETERS
  13. C           M      ORDER OF THE MATRIX
  14. C           A      THE INPUT MATRIX, M BY M
  15. C           RR     VECTOR CONTAINING THE REAL PARTS OF THE EIGENVALUES
  16. C                  ON RETURN
  17. C           RI     VECTOR CONTAINING THE IMAGINARY PARTS OF THE EIGEN-
  18. C                  VALUES ON RETURN
  19. C           IANA   VECTOR WHOSE DIMENSION MUST BE GREATER THAN OR EQUAL
  20. C                  TO M, CONTAINING ON RETURN INDICATIONS ABOUT THE WAY
  21. C                  THE EIGENVALUES APPEARED (SEE MATH. DESCRIPTION)
  22. C           IA     SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A
  23. C                  IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE
  24. C                  SUBSCRIPTED DATA STORAGE MODE.
  25. C                  IA=M WHEN THE MATRIX IS IN SSP VECTOR STORAGE MODE.
  26. C
  27. C        REMARKS
  28. C           THE ORIGINAL MATRIX IS DESTROYED
  29. C           THE DIMENSION OF RR AND RI MUST BE GREATER OR EQUAL TO M
  30. C
  31. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  32. C           NONE
  33. C
  34. C        METHOD
  35. C           QR DOUBLE ITERATION
  36. C
  37. C        REFERENCES
  38. C           J.G.F. FRANCIS - THE QR TRANSFORMATION---THE COMPUTER
  39. C           JOURNAL, VOL. 4, NO. 3, OCTOBER 1961, VOL. 4, NO. 4, JANUARY
  40. C           1962.  J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM -
  41. C           CLARENDON PRESS, OXFORD, 1965.
  42. C
  43. C     ..................................................................
  44. C
  45.       SUBROUTINE ATEIG(M,A,RR,RI,IANA,IA)
  46.       DIMENSION A(1),RR(1),RI(1),PRR(2),PRI(2),IANA(1)
  47.       INTEGER P,P1,Q
  48. C
  49.       E7=1.0E-8
  50.       E6=1.0E-6
  51.       E10=1.0E-10
  52.       DELTA=0.5
  53.       MAXIT=30
  54. C
  55. C        INITIALIZATION
  56. C
  57.       N=M
  58.    20 N1=N-1
  59.       IN=N1*IA
  60.       NN=IN+N
  61.       IF(N1) 30,1300,30
  62.    30 NP=N+1
  63. C
  64. C        ITERATION COUNTER
  65. C
  66.       IT=0
  67. C
  68. C        ROOTS OF THE 2ND ORDER MAIN SUBMATRIX AT THE PREVIOUS
  69. C        ITERATION
  70. C
  71.       DO 40 I=1,2
  72.       PRR(I)=0.0
  73.    40 PRI(I)=0.0
  74. C
  75. C        LAST TWO SUBDIAGONAL ELEMENTS AT THE PREVIOUS ITERATION
  76. C
  77.       PAN=0.0
  78.       PAN1=0.0
  79. C
  80. C        ORIGIN SHIFT
  81. C
  82.       R=0.0
  83.       S=0.0
  84. C
  85. C        ROOTS OF THE LOWER MAIN 2 BY 2 SUBMATRIX
  86. C
  87.       N2=N1-1
  88.       IN1=IN-IA
  89.       NN1=IN1+N
  90.       N1N=IN+N1
  91.       N1N1=IN1+N1
  92.    60 T=A(N1N1)-A(NN)
  93.       U=T*T
  94.       V=4.0*A(N1N)*A(NN1)
  95.       IF(ABS(V)-U*E7) 100,100,65
  96.    65 T=U+V
  97.       IF(ABS(T)-AMAX1(U,ABS(V))*E6) 67,67,68
  98.    67 T=0.0
  99.    68 U=(A(N1N1)+A(NN))/2.0
  100.       V=SQRT(ABS(T))/2.0
  101.       IF(T)140,70,70
  102.    70 IF(U) 80,75,75
  103.    75 RR(N1)=U+V
  104.       RR(N)=U-V
  105.       GO TO 130
  106.    80 RR(N1)=U-V
  107.       RR(N)=U+V
  108.       GO TO 130
  109.   100 IF(T)120,110,110
  110.   110 RR(N1)=A(N1N1)
  111.       RR(N)=A(NN)
  112.       GO TO 130
  113.   120 RR(N1)=A(NN)
  114.       RR(N)=A(N1N1)
  115.   130 RI(N)=0.0
  116.       RI(N1)=0.0
  117.       GO TO 160
  118.   140 RR(N1)=U
  119.       RR(N)=U
  120.       RI(N1)=V
  121.       RI(N)=-V
  122.   160 IF(N2)1280,1280,180
  123. C
  124. C        TESTS OF CONVERGENCE
  125. C
  126.   180 N1N2=N1N1-IA
  127.       RMOD=RR(N1)*RR(N1)+RI(N1)*RI(N1)
  128.       EPS=E10*SQRT(RMOD)
  129.       IF(ABS(A(N1N2))-EPS)1280,1280,240
  130.   240 IF(ABS(A(NN1))-E10*ABS(A(NN))) 1300,1300,250
  131.   250 IF(ABS(PAN1-A(N1N2))-ABS(A(N1N2))*E6) 1240,1240,260
  132.   260 IF(ABS(PAN-A(NN1))-ABS(A(NN1))*E6)1240,1240,300
  133.   300 IF(IT-MAXIT) 320,1240,1240
  134. C
  135. C        COMPUTE THE SHIFT
  136. C
  137.   320 J=1
  138.       DO 360 I=1,2
  139.       K=NP-I
  140.       IF(ABS(RR(K)-PRR(I))+ABS(RI(K)-PRI(I))-DELTA*(ABS(RR(K))
  141.      1    +ABS(RI(K)))) 340,360,360
  142.   340 J=J+I
  143.   360 CONTINUE
  144.       GO TO (440,460,460,480),J
  145.   440 R=0.0
  146.       S=0.0
  147.       GO TO 500
  148.   460 J=N+2-J
  149.       R=RR(J)*RR(J)
  150.       S=RR(J)+RR(J)
  151.       GO TO 500
  152.   480 R=RR(N)*RR(N1)-RI(N)*RI(N1)
  153.       S=RR(N)+RR(N1)
  154. C
  155. C        SAVE THE LAST TWO SUBDIAGONAL TERMS AND THE ROOTS OF THE
  156. C        SUBMATRIX BEFORE ITERATION
  157. C
  158.   500 PAN=A(NN1)
  159.       PAN1=A(N1N2)
  160.       DO 520 I=1,2
  161.       K=NP-I
  162.       PRR(I)=RR(K)
  163.   520 PRI(I)=RI(K)
  164. C
  165. C        SEARCH FOR A PARTITION OF THE MATRIX, DEFINED BY P AND Q
  166. C
  167.       P=N2
  168.       IF (N-3)600,600,525
  169.   525 IPI=N1N2
  170.       DO 580 J=2,N2
  171.       IPI=IPI-IA-1
  172.       IF(ABS(A(IPI))-EPS) 600,600,530
  173.   530 IPIP=IPI+IA
  174.       IPIP2=IPIP+IA
  175.       D=A(IPIP)*(A(IPIP)-S)+A(IPIP2)*A(IPIP+1)+R
  176.       IF(D)540,560,540
  177.   540 IF(ABS(A(IPI)*A(IPIP+1))*(ABS(A(IPIP)+A(IPIP2+1)-S)+ABS(A(IPIP2+2)
  178.      1 )) -ABS(D)*EPS) 620,620,560
  179.   560 P=N1-J
  180.   580 CONTINUE
  181.   600 Q=P
  182.       GO TO 680
  183.   620 P1=P-1
  184.       Q=P1
  185.       IF (P1-1) 680,680,650
  186.   650 DO 660 I=2, P1
  187.       IPI=IPI-IA-1
  188.       IF(ABS(A(IPI))-EPS)680,680,660
  189.   660 Q=Q-1
  190. C
  191. C        QR DOUBLE ITERATION
  192. C
  193.   680 II=(P-1)*IA+P
  194.       DO 1220 I=P,N1
  195.       II1=II-IA
  196.       IIP=II+IA
  197.       IF(I-P)720,700,720
  198.   700 IPI=II+1
  199.       IPIP=IIP+1
  200. C
  201. C        INITIALIZATION OF THE TRANSFORMATION
  202. C
  203.       G1=A(II)*(A(II)-S)+A(IIP)*A(IPI)+R
  204.       G2=A(IPI)*(A(IPIP)+A(II)-S)
  205.       G3=A(IPI)*A(IPIP+1)
  206.       A(IPI+1)=0.0
  207.       GO TO 780
  208.   720 G1=A(II1)
  209.       G2=A(II1+1)
  210.       IF(I-N2)740,740,760
  211.   740 G3=A(II1+2)
  212.       GO TO 780
  213.   760 G3=0.0
  214.   780 CAP=SQRT(G1*G1+G2*G2+G3*G3)
  215.       IF(CAP)800,860,800
  216.   800 IF(G1)820,840,840
  217.   820 CAP=-CAP
  218.   840 T=G1+CAP
  219.       PSI1=G2/T
  220.       PSI2=G3/T
  221.       ALPHA=2.0/(1.0+PSI1*PSI1+PSI2*PSI2)
  222.       GO TO 880
  223.   860 ALPHA=2.0
  224.       PSI1=0.0
  225.       PSI2=0.0
  226.   880 IF(I-Q)900,960,900
  227.   900 IF(I-P)920,940,920
  228.   920 A(II1)=-CAP
  229.       GO TO 960
  230.   940 A(II1)=-A(II1)
  231. C
  232. C        ROW OPERATION
  233. C
  234.   960 IJ=II
  235.       DO 1040 J=I,N
  236.       T=PSI1*A(IJ+1)
  237.       IF(I-N1)980,1000,1000
  238.   980 IP2J=IJ+2
  239.       T=T+PSI2*A(IP2J)
  240.  1000 ETA=ALPHA*(T+A(IJ))
  241.       A(IJ)=A(IJ)-ETA
  242.       A(IJ+1)=A(IJ+1)-PSI1*ETA
  243.       IF(I-N1)1020,1040,1040
  244.  1020 A(IP2J)=A(IP2J)-PSI2*ETA
  245.  1040 IJ=IJ+IA
  246. C
  247. C        COLUMN OPERATION
  248. C
  249.       IF(I-N1)1080,1060,1060
  250.  1060 K=N
  251.       GO TO 1100
  252.  1080 K=I+2
  253.  1100 IP=IIP-I
  254.       DO 1180 J=Q,K
  255.       JIP=IP+J
  256.       JI=JIP-IA
  257.       T=PSI1*A(JIP)
  258.       IF(I-N1)1120,1140,1140
  259.  1120 JIP2=JIP+IA
  260.       T=T+PSI2*A(JIP2)
  261.  1140 ETA=ALPHA*(T+A(JI))
  262.       A(JI)=A(JI)-ETA
  263.       A(JIP)=A(JIP)-ETA*PSI1
  264.       IF(I-N1)1160,1180,1180
  265.  1160 A(JIP2)=A(JIP2)-ETA*PSI2
  266.  1180 CONTINUE
  267.       IF(I-N2)1200,1220,1220
  268.  1200 JI=II+3
  269.       JIP=JI+IA
  270.       JIP2=JIP+IA
  271.       ETA=ALPHA*PSI2*A(JIP2)
  272.       A(JI)=-ETA
  273.       A(JIP)=-ETA*PSI1
  274.       A(JIP2)=A(JIP2)-ETA*PSI2
  275.  1220 II=IIP+1
  276.       IT=IT+1
  277.       GO TO 60
  278. C
  279. C        END OF ITERATION
  280. C
  281.  1240 IF(ABS(A(NN1))-ABS(A(N1N2))) 1300,1280,1280
  282. C
  283. C        TWO EIGENVALUES HAVE BEEN FOUND
  284. C
  285.  1280 IANA(N)=0
  286.       IANA(N1)=2
  287.       N=N2
  288.       IF(N2)1400,1400,20
  289. C
  290. C        ONE EIGENVALUE HAS BEEN FOUND
  291. C
  292.  1300 RR(N)=A(NN)
  293.       RI(N)=0.0
  294.       IANA(N)=1
  295.       IF(N1)1400,1400,1320
  296.  1320 N=N1
  297.       GO TO 20
  298.  1400 RETURN
  299.       END
  300.