C SOLUTION OF THE MATRIX EQUATION AX + XB = C. C C THIS IS A SAMPLE MAIN PROGRAM WHICH SOLVES RANDOMLY GENERATED C PROBLEMS. C C A IS REDUCED TO SCHUR FORM AND B TO UPPER-HESSENBERG FORM. C THE RESULTING BLOCK-SCHUR SYSTEM OF EQUATIONS IS C SOLVED AND THE SOLUTION TRANSFORMED BACK TO THAT OF THE C ORIGINAL EQUATION. C C APRIL 9, 1978. WRITTEN BY STEPHEN NASH C C REAL XX(10),YY(10),ZZ(10) REAL*8 EPS,ORT(10),A(10,10),B(10,10),C(10,10),V(10,10) * ,AA(10,10),BB(10,10),TT(10,10),DABS,T(10,10),DMAX,WORK(5310) c DATA EPS/Z3410000000000000/ C C WITHIN THE MAIN PROGRAM, A,B AND C ARE STORED NORMALLY. C AT THE END OF THE COMPUTATION C HOLDS THE SOLUTION X. THE ARRAY C WORK IS USED FOR TEMPORARY STORAGE WITHIN THE SUBROUTINES. C V IS USED TO STORE THE TRANSFORMATION MATRIX FOR A. C THE TRANSFORMATION MATRIX FOR B IS STORED IN ORT AND B. C C AND X ARE M * N MATRICES. A IS M * M, B IS N * N. C INTEGER IPR(200) eps = d1mach(4) NDIM = 10 MDIM = 10 MODE = 1 ISEED = 14365725 MAXIIJ = 2 DO 90 IIJJ = 1,MAXIIJ C C N AND M ARE RANDOMLY GENERATED IN THE RANGE U10,50e. C AFTERWARDS, THE MATRICES A,B, AND X (STORED IN T) ARE C ALSO GENERATED RANDOMLY WITH ENTRIES IN THE INTERVAL U0,1). C ISEED IS A SEED FOR THE RANDOM NUMBER GENERATOR. C IF IIJJ=1 WE ARE GENERATING A NEW PROBLEM. WHEN C IIJJ=2, THE RIGHT HAND SIDE IS CHANGED, BUT THE C COEFFICIENT MATRICES ARE LEFT UNCHANGED. C CALL RANDK(ISEED,XX(1),0) CALL RANDK(ISEED,XX(2),0) M = 10 N = 5 IF (IIJJ .EQ. 2) MODE = 2 NM = MAX0(N,M) WRITE(6,10) M,N 10 FORMAT('-SOLUTION OF',I3,' BY',I3,' MATRIX EQUATION') DO 20 I = 1,NM DO 15 K = 1,NM CALL RANDK(ISEED,XX(K),0) CALL RANDK(ISEED,YY(K),0) CALL RANDK(ISEED,ZZ(K),0) 15 CONTINUE DO 20 J = 1,NM IF (IIJJ .EQ. 2) GO TO 17 A(I,J) = XX(J) B(I,J) = YY(J) AA(I,J) = A(I,J) BB(I,J) = B(I,J) 17 T(I,J) = ZZ(J) 20 CONTINUE DO 30 I = 1,NM DO 30 J = 1,NM TT(I,J) = T(I,J) 30 CONTINUE C C GENERATE C FROM A, B, AND X C CALL SLVCHK(AA,BB,T,C,NDIM,N,MDIM,M) C C NOW SOLVE THE MATRIX EQUATION C CALL AXXBC(MODE,NDIM,N,MDIM,M,A,V,B,C,ORT,WORK,IPR,EPS,IERR) IF (IERR.NE.0) GO TO 60 C C CHECK THE SIZE OF THE ERROR. C DMAX = 0.D0 DO 40 I = 1,M DO 40 J = 1,N IF (DMAX.LT.DABS(C(I,J)-TT(I,J))) DMAX = DABS(C(I,J)-TT(I,J)) 40 CONTINUE WRITE(6,50) DMAX 50 FORMAT(' MAXIMAL ERROR = ', D14.6) GO TO 90 60 WRITE(6,70) 70 FORMAT('0SINGULAR SYSTEM OF EQUATIONS') 90 CONTINUE STOP END C C C C C SUBROUTINE RANDK (IY, YFL, INDEX) C C THIS IS A UNIFORM RANDOM NUMBER GENERATOR WRITTEN BY G. E. C FORSYTHE IN SPRING 1969, FOLLOWING D. KNUTH, THE ART OF COMPUTER C PROGRAMMING, VOL. 2, PP. 155-156. IT IS MUCH SUPERIOR TO RANDU, C THE RANDOM NUMBER GENERATOR FOUND IN IBM'S SCIENTIFIC SUBROUTINE C PACKAGE. C C BEFORE THE FIRST CALL OF RANDK, IY SHOULD BE SET OUTSIDE RANDK C TO AN ARBITRARY INTEGER VALUE. (IN WATFOR THIS IS ESSENTIAL.) C FOR PROGRAM CHECKOUT, THE INITIAL VALUE OF IY SHOULD BE A FIXED C INTEGER. FOR RANDOM NUMBERS DIFFERENT ON EVERY RUN (AND HENCE C NOT REPRODUCIBLE), DECLARE INTEGER CLOCK1 AND THEN INITIALIZE C IY TO CLOCK1(4). C C IF RANDK IS CALLED WITH AN INTEGER INDEX = 1, THEN THE OUTPUT C VALUE OF IY IS A PSEUDORANDOM INTEGER UNIFORMLY DISTRIBUTED IN THE C RANGE 0 <= IY < 2**31. C C IF RANDK IS CALLED WITH INDEX = 0, THEN NOT ONLY IS IY PRODUCED, C BUT ALSO (AT SOME EXTRA COST IN TIME) A FLOATING NUMBER YFL, UNI- C FORMLY DISTRIBUTED IN THE INTERVAL 0.0 <= YFL < 1.0. C IY = IY*314159269 + 453806245 4 IF (IY .GE. 0) GO TO 6 C C CAUTION: THE STATEMENT LABEL 4 IS ESSENTIAL IN ORDER TO PREVENT C CERTAIN COMPILERS (E.G., FORTRAN H WITH OPT 0) FROM PERFORMING C UNWANTED "OPTIMIZATIONS." IT SHOULD NOT BE REMOVED. C 5 IY = IY + 2147483647 + 1 C STATEMENT 5 ADDS 2**31 TO NEGATIVE VALUES OF IY C 6 IF (INDEX) 7, 7, 8 C 7 YFL = IY YFL = YFL*.4656613E-9 C 8 RETURN END C C C C C SUBROUTINE SLVCHK(A,B,X,C,NDIM,N,MDIM,M) C INTEGER I,J,K,M,MDIM,N,NDIM REAL*8 A(MDIM,M),B(NDIM,N),X(MDIM,N),C(MDIM,N) C C THIS SUBROUTINE FORMS AX + XB IN ORDER TO CHECK THE SOLUTION C TO THE PROBLEM. THE RESULT IS STORED IN THE MATRIX C. C PARAMETERS: C A,B,X - INPUT MATRICES C C - OUTPUT MATRIX C NDIM - DECLARED ROW DIMENSION OF THE MATRIX B C N - ACTUAL ROW DIMENSION OF B AND COLUMN DIMENSION C OF B, C, AND X C MDIM - DECLARED ROW DIMENSION OF THE MATRICES A, X, AND C C M - ACTUAL COLUMN DIMESION OF A AND ROW DIMENSION OF C A, X, AND C. C DO 10 I = 1,M DO 10 J = 1,N C(I,J) = 0.D0 DO 10 K = 1,M C(I,J) = C(I,J) + A(I,K) * X(K,J) 10 CONTINUE DO 20 I = 1,M DO 20 J = 1,N DO 20 K = 1,N C(I,J) = C(I,J) + X(I,K) * B(K,J) 20 CONTINUE RETURN END C C------------------------------------------------------------------------ C SUBROUTINE AXXBC(MODE,NDIM,N,MDIM,M,A,V,B,C,ORT,D,IPR,EPS,IERR) C INTEGER I,IERR,IND,IPR(1),J,M,MDIM,MODE,N,NDIM REAL*8 A(MDIM,M),B(NDIM,N),V(NDIM,N), * C(MDIM,N),ORT(1),EPS,REPS,D(1) C C THIS ROUTINE SOLVES THE MATRIX EQUATION AX + XB = C WHERE C C A - M * M MATRIX C X - M * N MATRIX C B - N * N MATRIX C C - M * N MATRIX. C C THE PARAMETERS FOR THIS SUBROUTINE ARE: C C MODE - INDICATES WHETHER A NEW MATRIX EQUATION IS BEING SOLVED, C OR THAT JUST THE RIGHT HAND SIDE HAS BEEN CHANGED. C MODE=1 IMPLIES THAT THE EQUATION IS NEW C MODE=2 IMPLIES THAT ONLY THE RIGHT HAND SIDE IS NEW C NDIM - DECLARED ROW DIMENSION OF THE MATRIX B C N - ROW DIMENSION OF B; COLUMN DIMENSION OF B AND C C MDIM - DECLARED ROW DIMENSION OF A, C, AND V C M - ROW DIMENSION OF A, C, AND V; COLUMN DIMENSION OF A AND V C A - M * M MATRIX (DOUBLE PRECISION) C B - N * N MATRIX (DOUBLE PRECISION) C V - N * N MATRIX (DOUBLE PRECISION) USED TO STORE THE MATRIX C WHICH TRANSFORMS B(T) TO UPPER SCHUR FORM C C - M * N MATRIX WHICH HOLDS C ON INPUT AND X ON OUTPUT C ORT - DOUBLE PRECISION VECTOR USED BY SUBROUTINE RG C OF LENGTH AT LEAST MAX(N,M) C D - DOUBLE PRECISION VECTOR FOR INTERNAL USE OF LENGTH C AT LEAST 2M**2 + 7M C IPR - INTEGER VECTOR FOR INTERNAL USE OF LENGTH AT LEAST 4N C EPS - ERROR TOLERANCE. EPS SHOULD EQUAL THE SMALLEST NUMBER C WHICH SATISFIES FL(1 + EPS) > 1. C IERR - INTEGER VARIABLE USED TO SET ERROR CODES. C IERR=0 => PROPER RETURN. C IERR=J => J-TH EIGENVALUE OF B HAS NOT BEEN DETERMINED C AFTER 30 QR ITERATIONS. C IERR=-J => A SINGULAR MATRIX WAS ENCOUNTERED WHEN C SOLVING FOR THE J-TH COLUMN OF X. C C METHOD: C C FIRST OF ALL, B(TRANSPOSE) IS FORMED AND IS USED INSTEAD OF C B THROUGHOUT THE PROGRAM. (IN THE FOLLOWING DISCUSSION, B(T) IS C USED TO REPRESENT B(TRANSPOSE)). C B(T) IS TRANSFORMED TO REAL UPPER SCHUR FORM AND A TO UPPER C HESSENBERG FORM USING ORTHOGONAL TRANSFORMATIONS. C THE RIGHT HAND SIDE C IS MULTIPLIED BY THESE TRANSFORMATION C MATRICES AND THE SOLUTION OF THIS TRANSFORMED SYSTEM IS C COMPUTED. THIS SOLUTION IS THEN MULTIPLIED BY THE TRANS- C FORMATION MATRICES IN ORDER TO OBTAIN THE SOLUTION TO THE C ORIGINAL PROBLEM. C C NOTE: C C IN ORDER TO GET OPTIMAL EFFICIENCY FROM THIS ROUTINE, N SHOULD C BE <= M. IF THIS IS NOT TRUE FOR YOUR PROBLEM, TRANSPOSE THE C EQUATION AND SOLVE THE RESULTING PROBLEM. C C ------------------------------------------------------------- C C FORM B(TRANSPOSE) FOR INTERNAL USE C IF (MODE .EQ. 2) GO TO 7 DO 5 I = 1,N DO 5 J = I,N D(1) = B(I,J) B(I,J) = B(J,I) B(J,I) = D(1) 5 CONTINUE C C A AND B(T) WILL BE TRANSFORMED INTO THE APPROPRIATE FORMS USING C MODIFICATIONS OF EISPACK ROUTINES (MODIFIED BY S. NASH). C C WARNING: DO NOT CHANGE THE ORDER OF THE NEXT TWO STATEMENTS C CALL RG(NDIM,N,B,1,V,ORT,EPS,IERR) CALL RG(MDIM,M,A,0,A,ORT,EPS,IERR) C C NOW, TRANSFORM THE RIGHT HAND SIDE. C 7 CALL TRANSF(A,ORT,1,C,V,0,M,N,MDIM,NDIM,D) C C NOW SOLVE THE SYSTEM OF EQUATIONS (BLOCK BACK SUBSTITUTION) C THE SYSTEM IS BLOCK UPPER HESSENBERG WITH UPPER HESSENBERG C MATRICES ON THE DIAGONAL (A + B(I,I) * I) AND MULTIPLES OF C THE IDENTITY OFF THE DIAGONAL. C TESTS ARE MADE TO SEE WHETHER A SINGLE OR A DOUBLE BLOCK C MUST BE HANDLED AT THE CURRENT STAGE AND APPROPRIATE C SUBROUTINES ARE CALLED. C C SET RELATIVE ERROR TOLERANCE C REPS = EPS*N*N*M*M IND = N - 1 IF (IND.EQ.0) GO TO 40 10 IF (DABS(B(IND + 1,IND)).LE.REPS) GO TO 20 CALL N2SOLV(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR) IF (IERR.NE.0) RETURN GO TO 30 20 CALL NSOLVE(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR) IF (IERR.NE.0) RETURN 30 IF (IND.GT.0) GO TO 10 40 IF (IND.EQ.0) CALL NSOLVE(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR) C C THE SOLUTION OF THE TRANSFORMED SYSTEM OF EQUATIONS C HAS NOW BEEN OBTAINED. IT IS TRANSFORMED BACK INTO C THE SOLUTION OF THE ORIGINAL SYSTEM OF EQUATIONS. C CALL TRANSF(A,ORT,0,C,V,1,M,N,MDIM,NDIM,D) RETURN END C C C C C SUBROUTINE RG(NDIM,N,A,MATZ,Z,ORT,EPS,IERR) C INTEGER N,NDIM,IERR,MATZ REAL*8 A(NDIM,N),Z(NDIM,N),ORT(N),EPS C C THIS SUBROUTINE REDUCES A TO EITHER UPPER-HESSENBERG C OR SCHUR FORM DEPENDING ON WHETHER MATZ IS ZERO OR NOT. C THE ROUTINES USED ARE FROM EISPACK EXCEPT THAT HQR2 HAS C BEEN MODIFIED SO AS NOT TO COMPUTE THE EIGENVECTORS OF C A AND ONLY REDUCES A TO SCHUR FORM. ORTHOGONAL TRANSFORMATIONS C ARE USED THROUGHOUT. C PARAMETERS: C NDIM - DECLARED ROW DIMENSION OF THE MATRIX A C N - ACTUAL DIMENSION OF THE MATRIX A C A - MATRIX TO BE TRANSFORMED (ON INPUT) AND TRANSFORMED C MATRIX (ON OUTPUT) C MATZ - INTEGER VARIABLE; MATZ = 0 => TRANSFORM A TO UPPER- C HESSENBERG FORM; OTHERWISE COMPUTE THE SCHUR FORM OF A C Z - STORES THE TRANSFORMATION MATRIX IF A IS TRANSFORMED TO C - SCHUR FORM (IF MATZ=0, Z CAN BE THE SAME AS A). C ORT - VECTOR USED FOR INTERNAL STORAGE; IF MATZ=0, ORT HOLDS C INFORMATION ABOUT THE TRANSFORMATION MATRIX C EPS - ERROR TOLERANCE, AS IN SUBROUTINE AXXBC. C IERR - ERROR CODE; SEE ROUTINE AXXBC FOR DETAILS C IERR=0 CALL ORTHES(NDIM,N,A,ORT) IF (MATZ.EQ.0) RETURN CALL ORTRAN(NDIM,N,A,ORT,Z) CALL HQR2(NDIM,N,A,Z,EPS,IERR) 10 RETURN END C C C C C SUBROUTINE ORTHES(NDIM,N,A,ORT) C INTEGER I,J,M,N,II,JJ,LA,MP,NDIM,KP1 REAL*8 A(NDIM,N),ORT(N) REAL*8 F,G,H,SCALE REAL*8 DSQRT,DABS,DSIGN C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTHES, C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). C C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS C 1 THROUGH N TO UPPER HESSENBERG FORM BY C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT: C C NDIM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C A CONTAINS THE INPUT MATRIX. C C ON OUTPUT: C C A CONTAINS THE HESSENBERG MATRIX. INFORMATION ABOUT C THE ORTHOGONAL TRANSFORMATIONS USED IN THE REDUCTION C IS STORED IN THE REMAINING TRIANGLE UNDER THE C HESSENBERG MATRIX; C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS. C ONLY ELEMENTS 1 THROUGH N ARE USED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ----------------------------------------------------------------- C LA = N - 1 KP1 = 2 IF (LA .LT. KP1) GO TO 200 C DO 180 M = KP1, LA H = 0.0D0 ORT(M) = 0.0D0 SCALE = 0.0D0 C C SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) C DO 90 I = M, N 90 SCALE = SCALE + DABS(A(I,M-1)) C IF (SCALE .EQ. 0.0D0) GO TO 180 MP = M + N C C FOR I=N STEP -1 UNTIL M DO -- C DO 100 II = M, N I = MP - II ORT(I) = A(I,M-1) / SCALE H = H + ORT(I) * ORT(I) 100 CONTINUE C G = -DSIGN(DSQRT(H),ORT(M)) H = H - ORT(M) * G ORT(M) = ORT(M) - G C C FORM (I-(U*UT)/H) * A C DO 130 J = M, N F = 0.0D0 C C FOR I=N STEP -1 UNTIL M DO -- C DO 110 II = M, N I = MP - II F = F + ORT(I) * A(I,J) 110 CONTINUE C F = F / H C DO 120 I = M, N 120 A(I,J) = A(I,J) - F * ORT(I) C 130 CONTINUE C C FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) C DO 160 I = 1, N F = 0.0D0 C C FOR J=N STEP -1 UNTIL M DO -- C DO 140 JJ = M, N J = MP - JJ F = F + ORT(J) * A(I,J) 140 CONTINUE C F = F / H C DO 150 J = M, N 150 A(I,J) = A(I,J) - F * ORT(J) C 160 CONTINUE C ORT(M) = SCALE * ORT(M) A(M,M-1) = SCALE * G 180 CONTINUE 200 RETURN END C C C C C SUBROUTINE ORTRAN(NDIM,N,A,ORT,Z) C INTEGER I,J,N,KL,MM,MP,NDIM,MP1 REAL*8 A(NDIM,N),ORT(N),Z(NDIM,N) REAL*8 G C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ORTRANS, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C C THIS SUBROUTINE ACCUMULATES THE ORTHOGONAL SIMILARITY C TRANSFORMATIONS USED IN THE REDUCTION OF A REAL GENERAL C MATRIX TO UPPER HESSENBERG FORM BY ORTHES. C C ON INPUT: C C NDIM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION BY ORTHES C IN ITS STRICT LOWER TRIANGLE; C C ORT CONTAINS FURTHER INFORMATION ABOUT THE TRANS- C FORMATIONS USED IN THE REDUCTION BY ORTHES. C ONLY ELEMENTS 1 THROUGH N ARE USED. C C ON OUTPUT: C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY ORTHES; C C ORT HAS BEEN ALTERED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ----------------------------------------------------------- C C INITIALIZE Z TO IDENTITY MATRIX C DO 80 I = 1, N C DO 60 J = 1, N 60 Z(I,J) = 0.0D0 C Z(I,I) = 1.0D0 80 CONTINUE C KL = N - 2 IF (KL .LT. 1) GO TO 200 C C FOR MP=N-1 STEP -1 UNTIL 2 DO -- C DO 140 MM = 1, KL MP = N - MM IF (A(MP,MP-1) .EQ. 0.0D0) GO TO 140 MP1 = MP + 1 C DO 100 I = MP1, N 100 ORT(I) = A(I,MP-1) C DO 130 J = MP, N G = 0.0D0 C DO 110 I = MP, N 110 G = G + ORT(I) * Z(I,J) C C DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES. C DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW C G = (G / ORT(MP)) / A(MP,MP-1) C DO 120 I = MP, N 120 Z(I,J) = Z(I,J) + G * ORT(I) C 130 CONTINUE C 140 CONTINUE C 200 RETURN END C C C C C SUBROUTINE HQR2(NDIM,N,H,Z,EPS,IERR) C INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NDIM,NN, X ITS,MP2,ENM2,IERR REAL*8 H(NDIM,N),Z(NDIM,N) REAL*8 P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,EPS REAL*8 DSQRT,DABS,DSIGN INTEGER MIN0 LOGICAL NOTLAS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). C IT HAS BEEN FURTHER MODIFIED BY STEPHEN NASH SO THAT IT ONLY C COMPUTES THE SCHUR FORM OF THE UPPER-HESSENBERG MATRIX H. C C ON INPUT: C C NDIM MUST BE SET TO THE ROW DIMENSION OF TWO - DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C H CONTAINS THE UPPER HESSENBERG MATRIX; C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ORTRAN C AFTER THE REDUCTION BY ORTHES. C C ON OUTPUT: C C H HAS BEEN DESTROYED; C C Z CONTAINS THE TRANSFORMATION MATRIX. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C ----------------------------------------------------------------- C C EPS IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C IERR = 0 NORM = 0.0D0 K = 1 C C COMPUTE MATRIX NORM C DO 50 I = 1, N C DO 40 J = K, N 40 NORM = NORM + DABS(H(I,J)) C K = I 50 CONTINUE C EN = N T = 0.0D0 C C SEARCH FOR NEXT EIGENVALUES C 60 IF (EN .LT. 1) GO TO 340 ITS = 0 NA = EN - 1 ENM2 = NA - 1 C C LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT C FOR L=EN STEP -1 UNTIL 1 DO -- C 70 DO 80 LL = 1, EN L = EN + 1 - LL IF (L .EQ. 1) GO TO 100 S = DABS(H(L-1,L-1)) + DABS(H(L,L)) IF (S .EQ. 0.0D0) S = NORM IF (DABS(H(L,L-1)) .LE. EPS * S) GO TO 100 80 CONTINUE C C FORM SHIFT C 100 X = H(EN,EN) IF (L .EQ. EN) GO TO 270 Y = H(NA,NA) W = H(EN,NA) * H(NA,EN) IF (L .EQ. NA) GO TO 280 IF (ITS .EQ. 30) GO TO 1000 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 C C FORM EXCEPTIONAL SHIFT C T = T + X C DO 120 I = 1, EN 120 H(I,I) = H(I,I) - X C S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) X = 0.75D0 * S Y = X W = -0.4375D0 * S * S 130 ITS = ITS + 1 C C LOOK FOR TWO CONSECUTIVE SMALL C SUB-DIAGONAL ELEMENTS. C FOR M=EN-2 STEP -1 UNTIL L DO -- C DO 140 MM = L, ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R * S - W) / H(M+1,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = DABS(P) + DABS(Q) + DABS(R) P = P / S Q = Q / S R = R / S IF (M .EQ. L) GO TO 150 IF (DABS(H(M,M-1)) * (DABS(Q) + DABS(R)) .LE. EPS * DABS(P) X * (DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I = MP2, EN H(I,I-2) = 0.0D0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.0D0 160 CONTINUE C C DOUBLE QR STEP INVOLVING ROWS L TO EN AND C COLUMNS M TO EN C DO 260 K = M, NA NOTLAS = K .NE. NA IF (K .EQ. M) GO TO 170 P = H(K,K-1) Q = H(K+1,K-1) R = 0.0D0 IF (NOTLAS) R = H(K+2,K-1) X = DABS(P) + DABS(Q) + DABS(R) IF (X .EQ. 0.0D0) GO TO 260 P = P / X Q = Q / X R = R / X 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) IF (K .EQ. M) GO TO 180 H(K,K-1) = -S * X GO TO 190 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 190 P = P + S X = P / S Y = Q / S ZZ = R / S Q = Q / P R = R / P C C ROW MODIFICATION C DO 210 J = K, N P = H(K,J) + Q * H(K+1,J) IF (.NOT. NOTLAS) GO TO 200 P = P + R * H(K+2,J) H(K+2,J) = H(K+2,J) - P * ZZ 200 H(K+1,J) = H(K+1,J) - P * Y H(K,J) = H(K,J) - P * X 210 CONTINUE C J = MIN0(EN,K+3) C C COLUMN MODIFICATION C DO 230 I = 1, J P = X * H(I,K) + Y * H(I,K+1) IF (.NOT. NOTLAS) GO TO 220 P = P + ZZ * H(I,K+2) H(I,K+2) = H(I,K+2) - P * R 220 H(I,K+1) = H(I,K+1) - P * Q H(I,K) = H(I,K) - P 230 CONTINUE C C ACCUMULATE TRANSFORMATIONS C DO 250 I = 1, N P = X * Z(I,K) + Y * Z(I,K+1) IF (.NOT. NOTLAS) GO TO 240 P = P + ZZ * Z(I,K+2) Z(I,K+2) = Z(I,K+2) - P * R 240 Z(I,K+1) = Z(I,K+1) - P * Q Z(I,K) = Z(I,K) - P 250 CONTINUE C 260 CONTINUE C GO TO 70 C C ONE ROOT FOUND C 270 H(EN,EN) = X + T EN = NA GO TO 60 C C TWO ROOTS FOUND C 280 P = (Y - X) / 2.0D0 Q = P * P + W ZZ = DSQRT(DABS(Q)) H(EN,EN) = X + T X = H(EN,EN) H(NA,NA) = Y + T IF (Q .LT. 0.0D0) GO TO 320 C C REAL PAIR C ZZ = P + DSIGN(ZZ,P) X = H(EN,NA) S = DABS(X) + DABS(ZZ) P = X / S Q = ZZ / S R = DSQRT(P*P+Q*Q) P = P / R Q = Q / R C C ROW MODIFICATION C DO 290 J = NA, N ZZ = H(NA,J) H(NA,J) = Q * ZZ + P * H(EN,J) H(EN,J) = Q * H(EN,J) - P * ZZ 290 CONTINUE C C COLUMN MODIFICATION C DO 300 I = 1, EN ZZ = H(I,NA) H(I,NA) = Q * ZZ + P * H(I,EN) H(I,EN) = Q * H(I,EN) - P * ZZ 300 CONTINUE C C ACCUMULATE TRANSFORMATIONS C DO 310 I = 1, N ZZ = Z(I,NA) Z(I,NA) = Q * ZZ + P * Z(I,EN) Z(I,EN) = Q * Z(I,EN) - P * ZZ 310 CONTINUE C GO TO 330 C C COMPLEX PAIR C 320 CONTINUE 330 EN = ENM2 GO TO 60 C C ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND C VECTORS OF UPPER TRIANGULAR FORM C 1000 IERR = EN 340 RETURN END C C C C C SUBROUTINE TRANSF(A,ORT,IT1,C,V,IT2,M,N,MDIM,NDIM,D) C INTEGER I,IT1,IT2,J,K,K1,K2,KK,M,MDIM,N,N2 REAL*8 V(NDIM,N),C(MDIM,N),A(MDIM,M),ORT(M),D(1),G C C SUBROUTINE TO FORM C(NEW) = U * C * V. THE INTEGER VARIABLES C IT1 AND IT2 INDICATE WHETHER U(TRANSPOSE) OR V(TRANSPOSE) C ARE REQUIRED INSTEAD OF U OR V. C IF ITX = 0 THEN THE ACTUAL MATRIX IS USED, IF ITX = 1, THEN C ITS TRANSPOSE IS USED. C IN ORDER TO SAVE STORAGE, THE MATRIX U (THE TRANSFORMATION C MATRIX FOR THE MATRIX A) IS STORED IN THE LOWER PART OF A C AND IN THE VECTOR ORT IN CODED FORM. (SEE SUBROUTINE ORTHES) C THE PARAMETERS FOR THIS SUBROUTINE ARE AS FOR THE SUBROUTINE C AXXBC EXCEPT THAT HERE ORT TAKES THE PLACE OF ORT. C C -------------------------------------------------------------- C C FORM U * C AND STORE IN C (AS IN SUBROUTINE ORTRAN) C M2 = M - 2 C C CHECK FOR SMALL MATRIX C IF (M2.LE.0) GO TO 45 DO 40 KK = 1,M2 K = M2 - KK + 1 IF (IT1.EQ.1) K = KK K1 = K + 1 IF (A(K1,K).EQ.0.D0) GO TO 40 D(K1) = ORT(K1) K2 = K + 2 DO 10 I = K2,M D(I) = A(I,K) 10 CONTINUE DO 30 J = 1,N G = 0.D0 DO 20 I = K1,M G = G + D(I) * C(I,J) 20 CONTINUE G = G / ORT(K1) / A(K1,K) DO 30 I = K1,M C(I,J) = C(I,J) + G * D(I) 30 CONTINUE 40 CONTINUE C C FORM C(NEW) = (U * C) * V C 45 DO 60 I = 1,M DO 50 J = 1,N D(J) = 0.D0 DO 50 K = 1,N IF (IT2.EQ.0) D(J) = D(J) + C(I,K) * V(K,J) IF (IT2.EQ.1) D(J) = D(J) + C(I,K) * V(J,K) 50 CONTINUE DO 60 J = 1,N C(I,J) = D(J) 60 CONTINUE RETURN END C C C C C SUBROUTINE NSOLVE(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR) C INTEGER I,I1,IERR,IND,IPR(1),IROW1,J,M,M1,MDIM,N,NDIM,MFIN REAL*8 A(MDIM,M),B(NDIM,N),C(MDIM,N),D(1),REPS C C THIS SUBROUTINE SETS UP AND SOLVES (WITH THE HELP OF C THE SUBROUTINE HESOLV) A SINGLE SYSTEM OF EQUATIONS C C ( A + B(IND+1,IND+1) * I )(X(IND+1)) = (C(IND+1)) C C (X(I),C(I) - COLUMN I OF THE MATRICES X AND C) C AND STORES THE RESULT IN THE APPROPRIATE ROW OF THE C MATRIX C. THE UPPER-HESSENBERG MATRIX IS STORED IN THE C VECTOR D IN ORDER TO SAVE STORAGE (ENTRIES BELOW THE SUB- C DIAGONAL ARE IGNORED). C THE RIGHT HAND SIDE IS ALSO STORED IN THE VECTOR D IN ORDER C TO REDUCE THE NUMBER OF VECTORS USED IN THE SUBROUTINES. C THE PARAMETERS ARE AS IN THE SUBROUTINE AXXBC. C C ---------------------------------------------------------- C C PERFORM BLOCK BACK-SUBSTITUTION IF NECESSARY C IF (IND.LT.N - 1) CALL BACKSB(C,B,IND,N,M,MDIM,NDIM) C C SET UP THE SYSTEM OF EQUATIONS C MFIN = (M * (M + 1)) / 2 + M DO 20 I = 1,M C C FIND INDEX OF BEGINNING OF ROW I WITHIN THE VECTOR D C M1 = IROW1(I,M) I1 = I - 1 IF (I.EQ.1) I1 = 1 DO 10 J = I1,M D(M1+J-I1+1) = A(I,J) 10 CONTINUE J = 2 IF (I.EQ.1) J = 1 D(M1+J) = D(M1+J) + B(IND+1,IND+1) C C STORE THE RIGHT HAND SIDE C D(MFIN+I) = C(I,IND+1) 20 CONTINUE C C SOLVE THE SYSTEM OF EQUATIONS AND STORE THE RESULT IN C. C CALL HESOLV(D,IPR,M,REPS,IERR) IF (IERR.NE.0) GO TO 40 DO 30 I = 1,M C(I,IND+1) = D(IPR(I)) 30 CONTINUE IND = IND - 1 RETURN 40 IERR = -IND - 1 RETURN END C C C C C SUBROUTINE HESOLV(D,IPR,M,REPS,IERR) C INTEGER I,I1,I2,IERR,IPR(1),IROW1,J,J1,K,N,N1,MFIN REAL*8 D(1),DABS,MULT,REPS C C THIS SUBROUTINE SOLVES AN UPPER-HESSENBERG SYSTEM OF EQUATIONS C WHICH IS STORED IN PACKED FORM IN THE VECTOR D. THE METHOD USED C IS GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING. THE PARAMETERS C IN THE SUBROUTINE ARE: C D - MATRIX OF COEFFICIENTS, RIGHT HAND SIDE (ON INPUT), C AND SOLUTION (ON OUTPUT) C IPR - INTEGER VECTOR WHICH RECORDS INTERCHANGES C M - SIZE OF SYSTEM C REPS - ERROR TOLERANCE C IERR - ERROR CODE (IERR=0 => NORMAL RETURN) C (IERR=-1 => SINGULAR MATRIX) C C ---------------------------------------------------------------- C C INITIALIZE INTERCHANGE VECTORS AND PARAMETERS. C IERR = 0 MFIN = (M * (M + 1)) / 2 + M DO 10 I = 1,M IPR(M+I) = IROW1(I,M) IPR(I) = I + MFIN 10 CONTINUE M1 = M - 1 C C REDUCE THE MATRIX TO UPPER TRIANGULAR FORM C IF (M.EQ.1) GO TO 35 DO 30 I = 1,M1 IF (DABS(D(IPR(M+I)+1)).GT.DABS(D(IPR(M+I+1)+1))) GO TO 20 C C INTERCHANGE ROWS IF NECESSARY C K = IPR(M+I) IPR(M+I) = IPR(M+I+1) IPR(M+I+1) = K K = IPR(I) IPR(I) = IPR(I+1) IPR(I+1) = K C C CHECK FOR COMPUTATIONALLY SINGULAR MATRIX C 20 IF (DABS(D(IPR(M+I)+1)).LT.REPS) GO TO 60 IPR(M+I+1) = IPR(M+I+1) + 1 C C ELIMINATE SUBDIAGONAL ELEMENT OF A C MULT = D(IPR(M+I+1)) / D(IPR(M+I)+1) D(IPR(I+1)) = D(IPR(I+1)) - MULT * D(IPR(I)) I1 = I + 1 DO 30 J = I1,M D(IPR(M+I+1)+J-I) = D(IPR(M+I+1)+J-I) - * MULT*D(IPR(M+I)+J+1-I) 30 CONTINUE 35 IF (DABS(D(IPR(M+M)+1)).LT.REPS) GO TO 60 C C PERFORM BACK - SUBSTITUTION C D(IPR(M)) = D(IPR(M)) / D(IPR(M+M)+1) IF (M1.EQ.0) RETURN DO 50 I1 = 1,M1 I = M - I1 I2 = I + 1 MULT = 0.D0 DO 40 J1 = I2,M J = J1 - I2 + 2 MULT = MULT + D(IPR(J1)) * D(IPR(M+I)+J) 40 CONTINUE D(IPR(I)) = (D(IPR(I)) - MULT) / D(IPR(M+I)+1) 50 CONTINUE RETURN 60 IERR = -1 RETURN END C C C C SUBROUTINE BACKSB(C,B,IND,N,M,MDIM,NDIM) C INTEGER I,IND,IND1,IND2,J,M,MDIM,N,NDIM REAL*8 B(NDIM,N),C(MDIM,N) C C BLOCK BACK - SUBSTITUTION FOR ONE 'ROW'. C PARAMETERS ARE AS IN SUBROUTINE AXXBC. C IND1 = IND + 1 IND2 = IND + 2 DO 10 I = IND2,N DO 10 J = 1,M C(J,IND1) = C(J,IND1) - B(IND1,I) * C(J,I) 10 CONTINUE RETURN END C C C C C INTEGER FUNCTION IROW1(I,M) C INTEGER I,M C C THIS FUNCTION COMPUTES THE INDEX (FOR THE VECTOR D) ONE BEFORE C THE BEGINNING OF ROW I IN THE M*M MATRIX. C IROW1 = (I-1) * M - (I-2) * (I-3) / 2 IF (I.EQ.1) IROW1 = 0 RETURN END C C C C C SUBROUTINE N2SOLV(A,B,C,D,NDIM,N,MDIM,M,IND,IPR,REPS,IERR) C INTEGER I,I1,IERR,IND,IPR(1),IROW2,J,J1,J2,K,LROW2,M,M1, * MDIM,N,NDIM,MFIN REAL*8 A(MDIM,M),B(NDIM,N),C(MDIM,N),D(1),REPS C C THIS SUBROUTINE SETS UP A DOUBLE SYSTEM OF EQUATIONS. C IT SOLVES (USING SUBROUTINE H2SOLV) THE SET OF EQUATIONS C C ( A+B(IND,IND)*I B(IND,IND+1) )( X(IND) ) ( C(IND) ) C ( )( ) = ( ) C ( B(IND+1,IND)*I A+B(IND+1,IND+1)*I )(X(IND+1)) (C(IND+1)) C C BY INTERCHANGING ROWS AND COLUMNS (TO INSURE THAT THERE ARE C ZEROES BELOW THE SUB-SUB-DIAGONAL) AND THEN USING GAUSSIAN C ELIMINATION. IN ORDER TO TAKE ADVANTAGE OF THE NUMBER OF C ZEROES IN THIS MATRIX, THE MATRIX IS STORED IN THE VECTOR C D IN PACKED FORM. (ENTRIES BELOW THE SUB-SUB-DIAGONAL ARE C IGNORED) THE RIGHT HAND SIDE AND SOLUTION ARE ALSO STORED IN D. C THE PARAMETERS ARE AS IN THE SUBROUTINE AXXBC. C C ------------------------------------------------------------ C C PERFORM BLOCK BACK-SUBSTITUTION IF NECESSARY C IF (IND.LT.N-1) CALL BACKS2(C,B,IND,N,M,MDIM,NDIM) C C SET UP THE SYSTEM OF EQUATIONS FOR H2SOLV C M2 = 2 * M MFIN = (M2 * (M2 + 1)) / 2 + 4 * M DO 20 I = 1,M C C FIND BEGINNING AND LENGTH OF ROW 2I-1 IN THE VECTOR D. C M1 = IROW2(2*I-1,M) K = LROW2(2*I-1,M) I1 = I - 1 IF (I.EQ.1) I1 = 1 DO 10 J = I1,M J1 = 2 * (J - I1 + 1) J2 = 1 IF (M1.EQ.0) J2 = 0 D(M1+J1-1) = A(I,J) D(M1+J1) = 0.D0 D(M1+K+J1-J2) = A(I,J) D(M1+K+J1-1-J2) = 0.D0 10 CONTINUE J1 = 3 IF (I.EQ.1) J1 = 1 D(J1+M1) = D(J1+M1) + B(IND,IND) D(J1+M1+1) = D(J1+M1+1) + B(IND,IND+1) IF (I.NE.1) J1 = 2 D(J1+M1+K) = D(J1+M1+K) + B(IND+1,IND) D(J1+M1+K+1) = D(J1+M1+K+1) + B(IND+1,IND+1) C C STORE RIGHT HAND SIDE C D(2*I+MFIN) = C(I,IND+1) D(2*I-1+MFIN) = C(I,IND) 20 CONTINUE C C SOLVE THE SYSTEM OF EQUATIONS AND STORE THE RESULT BACK IN C C CALL H2SOLV(D,IPR,M,REPS,IERR) IF (IERR.NE.0) GO TO 40 DO 30 I = 1,M C(I,IND) = D(IPR(2*I-1)) C(I,IND+1) = D(IPR(2*I)) 30 CONTINUE IND = IND - 2 RETURN 40 IERR = -IND - 1 RETURN END C C C C C SUBROUTINE H2SOLV(D,IPR,M,REPS,IERR) C INTEGER I,I1,I2,IERR,IPR(1),IROW2,J,J1,K,K1,L,N,N2,N21,MFIN REAL*8 REPS,D(1),MAX,DABS C C THIS SUBROUTINE SOLVES A SYSTEM OF 2M * 2M EQUATIONS C WHERE THE MATRIX HAS ZEROES BELOW THE SUB - SUB - DIAGONAL C AND WHERE THE MATRIX IS STORED AS A VECTOR IN PACKED FORM. C THE METHOD USED IS GAUSSIAN ELIMINATION. C PARAMETERS: C D - MATRIX OF COEFFICIENTS, RIGHT HAND SIDE (ON INPUT), C AND SOLUTION (ON OUTPUT). C IPR - VECTOR USED TO RECORD ROW PERMUTATIONS C M - SIZE OF SYSTEM C REPS - ERROR TOLERANCE C IERR - ERROR CODE (IERR=0 => NORMAL RETURN) C (IERR=-1 => SINGULAR MATRIX) C C -------------------------------------------------------------- C C INITIALIZE PIVOT VECTORS AND SIZE PARAMETERS. C IERR = 0 M2 = 2 * M MFIN = (M2 * (M2 + 1)) / 2 + 4 * M DO 10 I = 1,M2 IPR(M2+I) = IROW2(I,M) IPR(I) = I + MFIN 10 CONTINUE M21 = M2 - 1 C C TRANSFORM MATRIX TO TRIANGULAR FORM C DO 40 I = 1,M21 I1 = 2 IF (I.EQ.M21) I1 = 1 L = 0 C C INTERCHANGE ROWS IF NECESSARY C MAX = DABS(D(IPR(M2+I) + 1)) DO 20 J = 1,I1 IF (DABS(D(IPR(M2+J+I) + 1)).LE.MAX) GO TO 20 MAX = DABS(D(IPR(M2+J+I) + 1)) L = J 20 CONTINUE C C CHECK FOR COMPUTATIONALLY SINGULAR MATRIX C IF (MAX.LE.REPS) GO TO 80 IF (L.EQ.0) GO TO 30 K = IPR(M2+I) IPR(M2+I) = IPR(M2+L+I) IPR(M2+L+I) = K K = IPR(I) IPR(I) = IPR(I+L) IPR(I+L) = K 30 CONTINUE C C ADJUST POINTERS TO BEGINNINGS OF ROWS C IPR(M2+I+1) = IPR(M2+I+1) + 1 IF (I.NE.M21) IPR(M2+I+2) = IPR(M2+I+2) + 1 IP1 = I + 1 C C ELIMINATE SUBDIAGONAL ELEMENTS IN THE MATRIX C DO 40 J = 1,I1 MAX = D(IPR(M2+I+J)) / D(IPR(M2+I)+1) D(IPR(I+J)) = D(IPR(I+J)) - MAX * D(IPR(I)) DO 40 K1 = IP1,M2 K = K1 - I D(IPR(M2+I+J)+K) = D(IPR(M2+I+J)+K) - * MAX * D(IPR(M2+I)+1+K) 40 CONTINUE IF (DABS(D(IPR(M2+M2)+1)).LE.REPS) GO TO 80 C C PERFORM BACK SUBSTITUTION C D(IPR(M2)) = D(IPR(M2)) / D(IPR(M2+M2)+1) DO 60 I1 = 1,M21 I = M2 - I1 I2 = I + 1 MAX = 0.D0 DO 50 J1 = I2,M2 J = J1 - I2 + 2 MAX = MAX + D(IPR(J1)) * D(IPR(M2+I)+J) 50 CONTINUE D(IPR(I)) = (D(IPR(I)) - MAX) / D(IPR(M2+I)+1) 60 CONTINUE 70 RETURN 80 IERR = -1 GO TO 70 END C C C C C SUBROUTINE BACKS2(C,B,IND,N,M,MDIM,NDIM) C INTEGER I,IND,IND1,IND2,J,M,MDIM,N,NDIM REAL*8 B(NDIM,N),C(MDIM,N) C C BLOCK BACK - SUBSTITUTION FOR TWO 'ROWS'. C PARAMETERS ARE AS IN SUBROUTINE AXXBC. C IND1 = IND + 1 IND2 = IND + 2 DO 10 I = IND2,N DO 10 J = 1,M C(J,IND1) = C(J,IND1) - B(IND1,I) * C(J,I) C(J,IND) = C(J,IND) - B(IND,I) * C(J,I) 10 CONTINUE RETURN END C C C C C INTEGER FUNCTION IROW2(I,M) C INTEGER I,K,M C C THIS FUNCTION FINDS THE INDEX (FOR THE VECTOR D) ONE C BEFORE THE BEGINNING OF ROW I IN THE 2N*2N MATRIX. C IROW2 = 2 * (I - 1) * M K = (I - 4) * (I - 3) / 2 IF (I.LE.2) K = 0 IROW2 = IROW2 - K RETURN END C C C C C INTEGER FUNCTION LROW2(I,M) C INTEGER I,K,M C C THIS FUNCTION FINDS THE LENGTH OF ROW I IN THE C DOUBLE MATRIX (WHEN IT IS STORED IN PACKED FORM IN THE C VECTOR D) C K = I - 3 IF (I.LE.2) K = 0 LROW2 = 2 * M - K IF (I.LE.2) LROW2 = 2 * M RETURN END