C ALGORITHM 598, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2, C JUN., 1983, P. 246-254. COMPLEX A(50,50), B(50,50), C(50,50), S(50,50) MAN 10 COMPLEX WORK(17550) MAN 20 COMPLEX CMPLX MAN 30 REAL REAL MAN 40 REAL FXNORM, TOL MAN 50 INTEGER NM, N, I, ICASE, IERR, IGUESS, J, MAXITS, NW MAN 60 NM = 50 MAN 70 NW = 17550 MAN 80 C MAN 90 C CASE 1: DENNIS, TRAUB AND WEBER. MAN 100 C MAN 110 ICASE = 1 MAN 120 NOUT = 6 MAN 130 WRITE (NOUT,99999) ICASE MAN 140 N = 2 MAN 150 A(1,1) = CMPLX(1.0,0.0) MAN 160 A(1,2) = CMPLX(0.0,0.0) MAN 170 A(2,1) = CMPLX(0.0,0.0) MAN 180 A(2,2) = CMPLX(1.0,0.0) MAN 190 B(1,1) = CMPLX(7.0,0.0) MAN 200 B(1,2) = CMPLX(8.0,0.0) MAN 210 B(2,1) = CMPLX(8.0,0.0) MAN 220 B(2,2) = CMPLX(10.0,0.0) MAN 230 C(1,1) = CMPLX(9.0,0.0) MAN 240 C(1,2) = CMPLX(3.0,0.0) MAN 250 C(2,1) = CMPLX(4.0,0.0) MAN 260 C(2,2) = CMPLX(4.0,0.0) MAN 270 IGUESS = 0 MAN 280 TOL = 0.0 MAN 290 MAXITS = 0 MAN 300 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 310 * IERR) MAN 320 IF (IERR.NE.0) GO TO 10 MAN 330 FXNORM = REAL(WORK(1)) MAN 340 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 350 CALL PRINT(NM, N, S) MAN 360 GO TO 20 MAN 370 10 WRITE (NOUT,99998) ICASE, IERR MAN 380 C MAN 390 C CASE 2: ONE-BY-ONE TEST. MAN 400 C MAN 410 20 ICASE = 2 MAN 420 WRITE (NOUT,99999) ICASE MAN 430 N = 1 MAN 440 A(1,1) = CMPLX(1.0,0.0) MAN 450 B(1,1) = CMPLX(-3.0,0.0) MAN 460 C(1,1) = CMPLX(2.0,0.0) MAN 470 IGUESS = 0 MAN 480 TOL = 0.0 MAN 490 MAXITS = 0 MAN 500 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 510 * IERR) MAN 520 IF (IERR.NE.0) GO TO 30 MAN 530 FXNORM = REAL(WORK(1)) MAN 540 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 550 CALL PRINT(NM, N, S) MAN 560 GO TO 40 MAN 570 30 WRITE (NOUT,99998) ICASE, IERR MAN 580 C MAN 590 C CASE 3: COMPLEX SOLVENT. MAN 600 C MAN 610 40 ICASE = 3 MAN 620 WRITE (NOUT,99999) ICASE MAN 630 N = 3 MAN 640 A(1,1) = CMPLX(1.0,0.0) MAN 650 A(1,2) = CMPLX(2.0,0.0) MAN 660 A(1,3) = CMPLX(-1.0,0.0) MAN 670 A(2,1) = CMPLX(3.0,0.0) MAN 680 A(2,2) = CMPLX(4.0,0.0) MAN 690 A(2,3) = CMPLX(0.0,0.0) MAN 700 A(3,1) = CMPLX(-1.0,0.0) MAN 710 A(3,2) = CMPLX(2.0,0.0) MAN 720 A(3,3) = CMPLX(1.0,0.0) MAN 730 B(1,1) = CMPLX(-2.0,0.0) MAN 740 B(1,2) = CMPLX(0.0,1.0) MAN 750 B(1,3) = CMPLX(0.0,0.0) MAN 760 B(2,1) = CMPLX(0.0,1.0) MAN 770 B(2,2) = CMPLX(2.0,0.0) MAN 780 B(2,3) = CMPLX(0.0,0.0) MAN 790 B(3,1) = CMPLX(0.0,0.0) MAN 800 B(3,2) = CMPLX(0.0,0.0) MAN 810 B(3,3) = CMPLX(0.0,1.0) MAN 820 C(1,1) = CMPLX(2.0,0.0) MAN 830 C(1,2) = CMPLX(1.0,-7.0) MAN 840 C(1,3) = CMPLX(4.0,0.0) MAN 850 C(2,1) = CMPLX(1.0,-7.0) MAN 860 C(2,2) = CMPLX(-8.0,-17.0) MAN 870 C(2,3) = CMPLX(0.0,0.0) MAN 880 C(3,1) = CMPLX(0.0,2.0) MAN 890 C(3,2) = CMPLX(2.0,-2.0) MAN 900 C(3,3) = CMPLX(-4.0,-2.0) MAN 910 IGUESS = 0 MAN 920 TOL = 0.001 MAN 930 MAXITS = 0 MAN 940 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 950 * IERR) MAN 960 IF (IERR.NE.0) GO TO 50 MAN 970 FXNORM = REAL(WORK(1)) MAN 980 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 990 CALL PRINT(NM, N, S) MAN 1000 GO TO 60 MAN 1010 50 WRITE (NOUT,99998) ICASE, IERR MAN 1020 C MAN 1030 C CASE 4: NO SOLVENT. MAN 1040 C MAN 1050 60 ICASE = 4 MAN 1060 WRITE (NOUT,99999) ICASE MAN 1070 N = 2 MAN 1080 A(1,1) = CMPLX(1.0,0.0) MAN 1090 A(1,2) = CMPLX(0.0,0.0) MAN 1100 A(2,1) = CMPLX(0.0,0.0) MAN 1110 A(2,2) = CMPLX(0.0,0.0) MAN 1120 B(1,1) = CMPLX(0.0,0.0) MAN 1130 B(1,2) = CMPLX(0.0,0.0) MAN 1140 B(2,1) = CMPLX(0.0,0.0) MAN 1150 B(2,2) = CMPLX(0.0,0.0) MAN 1160 C(1,1) = CMPLX(1.0,0.0) MAN 1170 C(1,2) = CMPLX(0.0,0.0) MAN 1180 C(2,1) = CMPLX(0.0,0.0) MAN 1190 C(2,2) = CMPLX(1.0,0.0) MAN 1200 IGUESS = 0 MAN 1210 TOL = 0.0 MAN 1220 MAXITS = 0 MAN 1230 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 1240 * IERR) MAN 1250 IF (IERR.NE.0) GO TO 70 MAN 1260 FXNORM = REAL(WORK(1)) MAN 1270 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 1280 CALL PRINT(NM, N, S) MAN 1290 GO TO 80 MAN 1300 70 WRITE (NOUT,99998) ICASE, IERR MAN 1310 C MAN 1320 C CASE 5: EXAMPLE CITED IN PAPER. MAN 1330 C MAN 1340 80 ICASE = 5 MAN 1350 WRITE (NOUT,99999) ICASE MAN 1360 N = 2 MAN 1370 A(1,1) = CMPLX(1.0,0.0) MAN 1380 A(1,2) = CMPLX(0.0,0.0) MAN 1390 A(2,1) = CMPLX(0.0,0.0) MAN 1400 A(2,2) = CMPLX(1.0,0.0) MAN 1410 B(1,1) = CMPLX(-1.0,0.0) MAN 1420 B(1,2) = CMPLX(-1.0,0.0) MAN 1430 B(2,1) = CMPLX(1.0,0.0) MAN 1440 B(2,2) = CMPLX(-1.0,0.0) MAN 1450 C(1,1) = CMPLX(0.0,0.0) MAN 1460 C(1,2) = CMPLX(1.0,0.0) MAN 1470 C(2,1) = CMPLX(-1.0,0.0) MAN 1480 C(2,2) = CMPLX(0.0,0.0) MAN 1490 IGUESS = 0 MAN 1500 TOL = 0.0 MAN 1510 MAXITS = 0 MAN 1520 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 1530 * IERR) MAN 1540 IF (IERR.NE.0) GO TO 90 MAN 1550 FXNORM = REAL(WORK(1)) MAN 1560 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 1570 CALL PRINT(NM, N, S) MAN 1580 GO TO 100 MAN 1590 90 WRITE (NOUT,99998) ICASE, IERR MAN 1600 C MAN 1610 C CASE 6: COMPLEX GUESS FOR CASE 5. MAN 1620 C MAN 1630 100 ICASE = 6 MAN 1640 WRITE (NOUT,99999) ICASE MAN 1650 N = 2 MAN 1660 A(1,1) = CMPLX(1.0,0.0) MAN 1670 A(1,2) = CMPLX(0.0,0.0) MAN 1680 A(2,1) = CMPLX(0.0,0.0) MAN 1690 A(2,2) = CMPLX(1.0,0.0) MAN 1700 B(1,1) = CMPLX(-1.0,0.0) MAN 1710 B(1,2) = CMPLX(-1.0,0.0) MAN 1720 B(2,1) = CMPLX(1.0,0.0) MAN 1730 B(2,2) = CMPLX(-1.0,0.0) MAN 1740 C(1,1) = CMPLX(0.0,0.0) MAN 1750 C(1,2) = CMPLX(1.0,0.0) MAN 1760 C(2,1) = CMPLX(-1.0,0.0) MAN 1770 C(2,2) = CMPLX(0.0,0.0) MAN 1780 IGUESS = 1 MAN 1790 DO 120 I=1,N MAN 1800 DO 110 J=1,N MAN 1810 S(I,J) = CMPLX(0.0,0.0) MAN 1820 110 CONTINUE MAN 1830 S(I,I) = CMPLX(0.0,1.0) MAN 1840 120 CONTINUE MAN 1850 TOL = 0.0 MAN 1860 MAXITS = 0 MAN 1870 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 1880 * IERR) MAN 1890 IF (IERR.NE.0) GO TO 130 MAN 1900 FXNORM = REAL(WORK(1)) MAN 1910 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 1920 CALL PRINT(NM, N, S) MAN 1930 GO TO 140 MAN 1940 130 WRITE (NOUT,99998) ICASE, IERR MAN 1950 C MAN 1960 C CASE 7: BAD GUESS FOR CASE 5. MAN 1970 C MAN 1980 140 ICASE = 7 MAN 1990 WRITE (NOUT,99999) ICASE MAN 2000 N = 2 MAN 2010 A(1,1) = CMPLX(1.0,0.0) MAN 2020 A(1,2) = CMPLX(0.0,0.0) MAN 2030 A(2,1) = CMPLX(0.0,0.0) MAN 2040 A(2,2) = CMPLX(1.0,0.0) MAN 2050 B(1,1) = CMPLX(-1.0,0.0) MAN 2060 B(1,2) = CMPLX(-1.0,0.0) MAN 2070 B(2,1) = CMPLX(1.0,0.0) MAN 2080 B(2,2) = CMPLX(-1.0,0.0) MAN 2090 C(1,1) = CMPLX(0.0,0.0) MAN 2100 C(1,2) = CMPLX(1.0,0.0) MAN 2110 C(2,1) = CMPLX(-1.0,0.0) MAN 2120 C(2,2) = CMPLX(0.0,0.0) MAN 2130 IGUESS = 1 MAN 2140 DO 160 I=1,N MAN 2150 DO 150 J=1,N MAN 2160 S(I,J) = CMPLX(500.0,0.0) MAN 2170 150 CONTINUE MAN 2180 160 CONTINUE MAN 2190 TOL = 0.0 MAN 2200 MAXITS = 90 MAN 2210 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 2220 * IERR) MAN 2230 IF (IERR.NE.0) GO TO 170 MAN 2240 FXNORM = REAL(WORK(1)) MAN 2250 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 2260 CALL PRINT(NM, N, S) MAN 2270 GO TO 180 MAN 2280 170 WRITE (NOUT,99998) ICASE, IERR MAN 2290 C MAN 2300 C CASE 8: SINGULAR A,C,SOLVENT AND GUESS. MAN 2310 C MAN 2320 180 ICASE = 8 MAN 2330 WRITE (NOUT,99999) ICASE MAN 2340 N = 2 MAN 2350 A(1,1) = CMPLX(1.0E-8,0.0) MAN 2360 A(1,2) = CMPLX(0.0,0.0) MAN 2370 A(2,1) = CMPLX(0.0,0.0) MAN 2380 A(2,2) = CMPLX(1.0,0.0) MAN 2390 B(1,1) = CMPLX(1.0,0.0) MAN 2400 B(1,2) = CMPLX(0.0,0.0) MAN 2410 B(2,1) = CMPLX(0.0,0.0) MAN 2420 B(2,2) = CMPLX(1.0,0.0) MAN 2430 C(1,1) = CMPLX(0.0,0.0) MAN 2440 C(1,2) = CMPLX(0.0,0.0) MAN 2450 C(2,1) = CMPLX(0.0,0.0) MAN 2460 C(2,2) = CMPLX(-2.0,0.0) MAN 2470 IGUESS = 1 MAN 2480 S(1,1) = CMPLX(1.0,0.0) MAN 2490 S(1,2) = CMPLX(1.0,0.0) MAN 2500 S(2,1) = CMPLX(1.0,0.0) MAN 2510 S(2,2) = CMPLX(1.0,0.0) MAN 2520 TOL = 0.0 MAN 2530 MAXITS = 0 MAN 2540 CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 2550 * IERR) MAN 2560 IF (IERR.NE.0) GO TO 190 MAN 2570 FXNORM = REAL(WORK(1)) MAN 2580 WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 2590 CALL PRINT(NM, N, S) MAN 2600 GO TO 200 MAN 2610 190 WRITE (NOUT,99998) ICASE, IERR MAN 2620 C MAN 2630 200 CONTINUE MAN 2640 99999 FORMAT (//6H CASE , I2, 18H FOLLOWS..........) MAN 2650 99998 FORMAT (5H CASE, I2, 8H IERR IS, I3, 1H.) MAN 2660 99997 FORMAT (15H CONVERGENCE IN, I3, 23H ITERATIONS. NORM(F(X(, I2, MAN 2670 * 6H))) = , E15.6) MAN 2680 STOP MAN 2690 END MAN 2700 SUBROUTINE PRINT(NM, N, A) PRI 10 C C C SUBROUTINE PRINT PRINTS OUT AN INPUT MATRIX. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF A IN THE MAIN PROGRAM. C C N IS THE ORDER OF THE MATRIX A. C C A IS THE MATRIX TO BE PRINTED. C C C ON RETURN, C C NM, N AND A ARE UNALTERED. C C INTEGER I, J, JJ, JP1, N, NM COMPLEX A(NM,N) DO 20 J=1,N,2 JP1 = J + 1 IF (JP1.GT.N) JP1 = N NOUT = 6 WRITE (NOUT,99999) J, JP1 DO 10 I=1,N WRITE (NOUT,99998) (A(I,JJ),JJ=J,JP1) 10 CONTINUE IF (JP1.EQ.N) RETURN 20 CONTINUE 99999 FORMAT (//8H COLUMNS, I3, 8H THROUGH, I3, 11H ..........) 99998 FORMAT (/3(E15.6, 2H +, E15.6, 3H *I, 2X)) RETURN END SUBROUTINE SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, SQU 10 * MAXITS, IERR) C C C SUBROUTINE SQUINT BREAKS DOWN THE WORK ARRAY 'WORK' INTO C SMALLER PIECES. THE ACTUAL SOLUTION TO AX**2 + BX + C = 0 IS C DONE IN SUBROUTINE SQUIN2. THIS SUBROUTINE MERELY RELIEVES C THE USER FROM A LONG CALLING SEQUENCE. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF ALL THE MATRICES C IN THE CALLING PROGRAM. C C N IS THE ORDER OF THE MATRICES A, B AND C. C C A IS THE MATRIX COEFFICIENT OF X**2. C C B IS THE MATRIX COEFFICIENT OF X. C C C IS THE CONSTANT MATRIX. C C IGUESS IS AN INTEGER. IF IGUESS.NE.0, THE USER SUPPLIES AN C INITIAL GUESS AT A SOLVENT. THIS GUESS IS STORED IN C ARRAY S. IF IGUESS.EQ.0, THE SUBROUTINE PROVIDES ITS C OWN INITIAL GUESS. C C S CONTAINS THE USERS INITIAL GUESS AT A SOLVENT, IF IGUESS C HAS BEEN SET TO A NONZERO QUANTITY. OTHERWISE THE INPUT C CONTENTS IN S ARE IGNORED. C C WORK IS A WORK VECTOR. IT MUST BE DIMENSIONED AT LEAST C (7N**2 + N), WHERE N IS THE ORDER OF A, B, C AND S. C C NW IS THE DIMENSION OF THE ARRAY WORK IN THE CALLING C PROGRAM. C C TOL IS A USER-SUPPLIED ACCURACY TOLERANCE. SETTING TOL = 0.0 C CAUSES ITERATION TO PROCEED UNTIL FULL MACHINE PRECISION C IS ATTAINED. OTHERWISE, EXECUTION TERMINATES WHEN C NORM(AS**2+BS+C).LT.TOL. C C MAXITS IS AN INTEGER. IF MAXITS.NE.0, THE USER SPECIFIES THE C MOST INTERATIONS THE ALGORITHM IS TO TAKE. IF MAXITS. C .LE.0, IT IS RESET TO 30. C C C ON RETURN, C C A,B,C ARE DESTROYED. C C IGUESS CONTAINS THE NUMBER OF ITERATIONS PERFORMED TO C COMPUTE S. C C S CONTAINS THE RIGHT SOLVENT. C C WORK(1) IS A COMPLEX NUMBER WITH REAL PART EQUAL TO THE NORM C OF AS**2+BS+C. C C IERR IS AN INTEGER ERROR RETURN. C C IERR = 0 FOR A NORMAL RETURN. C C IERR = 1 INDICATES FAILURE OF SQUINT TO CONVERGE TO C A SOLVENT IN THE MAXIMUM NUMBER OF ITERATIONS. C C IERR = 2 INDICATES FAILURE IN THE UPPER REDUCTION IN C CQZIT. C C IERR = 3 INDICATES FAILURE IN THE LOWER REDUCTION IN C CQZIT. C C IERR = 10 + N INDICATES AN ERROR RETURN FROM TRISLV C ON ITERATION N, DESIGNATING INCONSISTENCY OF C THE TRIANGULAR SYSTEM. C C IERR = 999 INDICATES IMPROPER DIMENSIONING. THE C CONDITIONS NM.GE.N.GT.0 AND C NW.GE.(7*N*N + N) MUST HOLD. C INTEGER NM, N, IGUESS, NW, MAXITS, IERR, I1, I2, I3, I4, I5, I6, * I7 COMPLEX WORK(NW) COMPLEX A(NM,N), B(NM,N), C(NM,N), S(NM,N) REAL TOL I1 = N*N + 1 I2 = N*N + I1 I3 = N*N + I2 I4 = N*N + I3 I5 = N*N + I4 I6 = N*N + I5 I7 = N*N + I6 CALL SQUIN2(NM, N, A, B, C, IGUESS, S, WORK(1), WORK(I1), * WORK(I2), WORK(I3), WORK(I4), WORK(I5), WORK(I6), WORK(I7), NW, * TOL, MAXITS, IERR) RETURN END C SQU 10 SUBROUTINE SQUIN2(NM, N, A, B, C, IGUESS, S, L, U, V, Z, R, XOLD, SQU 20 * EYE, TEMP, NW, TOL, MAXITS, IERR) C C SUBROUTINE SQUIN2 FINDS A RIGHT SOLVENT OF THE MATRIX C EQUATION AX**2 + BX + C = 0. C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF ALL THE MATRICES IN C THE CALLING PROGRAM. C C N IS THE ORDER OF THE MATRICES A, B AND C. C C A IS THE MATRIX COEFFICIENT OF X**2. C C B IS THE MATRIX COEFFICIENT OF X. C C C IS THE CONSTANT MATRIX. C C IGUESS IS AN INTEGER SET IN THE CALL TO SQUINT. C C S IS A MATRIX SET IN THE CALL TO SQUINT. C C NW, TOL, AND MAXITS ARE INTEGER AND REAL PARAMETERS SET IN C THE CALL TO SQUINT. C C C THE FOLLOWING ARE INTERNAL VARIABLES: C C C L IS A MATRIX CONTAINING THE ITERATE X(I) FOR C REDUCTION TO LOWER TRIANGULAR FORM. C C U IS A MATRIX CONTAINING AX(I) + B FOR REDUCTION C TO UPPER TRIANGULAR FORM. C C V IS A MATRIX CONTAINING A FOR REDUCTION TO UPPER C TRIANGULAR FORM. C C Z,R ARE MATRICES CONTAINING THE HISTORY OF THE C TRANSFORMATIONS IN THE REDUCTIONS. C C XOLD IS A MATRIX HOLDING THE CURRENT ITERATE X(I). C C EYE CONTAINS AN IDENTITY MATRIX FOR THE LOWER REDUCTION C STEP. C C TEMP IS A WORK VECTOR. C C C ON RETURN, C C A, B, C, IGUESS, S AND IERR HAVE THE SAME PROPERTIES AS C DESCRIBED IN THE RETURN FROM SUBROUTINE SQUINT. C C L(1,1) IS A COMPLEX NUMBER WITH REAL PART EQUAL TO THE NORM C OF AS**2+BS+C. C C INTEGER NM, N, IGUESS, NW, MAXITS, IERR INTEGER ITS, MATS, I, J, K COMPLEX A(NM,N), B(NM,N), C(NM,N), S(NM,N), L(NM,N), U(NM,N) COMPLEX V(NM,N), Z(NM,N), R(NM,N), XOLD(NM,N), EYE(NM,N), TEMP(N) COMPLEX CMPLX, CONJG REAL ANORM, ANI, BNORM, BNI, CNORM, CNI, XNORM, XNI REAL FXNORM, FXNI, GNORM, TNORM, T, TOL REAL SQRT, CABS, FLOAT K = 7*N*N + N IF (NW.LT.K) GO TO 460 IF (NM.LT.N) GO TO 460 IF (N.LE.0) GO TO 460 IF (MAXITS.LE.0) MAXITS = 30 C ********** INITIALIZE ARRAYS ********** DO 20 I=1,N DO 10 J=1,N L(I,J) = CMPLX(0.0,0.0) U(I,J) = CMPLX(0.0,0.0) V(I,J) = CMPLX(0.0,0.0) Z(I,J) = CMPLX(0.0,0.0) XOLD(I,J) = CMPLX(0.0,0.0) EYE(I,J) = CMPLX(0.0,0.0) 10 CONTINUE TEMP(I) = CMPLX(0.0,0.0) 20 CONTINUE C ********** SET INITIAL GUESS(ES) ********** ANORM = 0.0 BNORM = 0.0 CNORM = 0.0 DO 40 I=1,N ANI = 0.0 BNI = 0.0 CNI = 0.0 DO 30 J=1,N ANI = ANI + CABS(A(I,J)) BNI = BNI + CABS(B(I,J)) CNI = CNI + CABS(C(I,J)) 30 CONTINUE IF (ANI.GT.ANORM) ANORM = ANI IF (BNI.GT.BNORM) BNORM = BNI IF (CNI.GT.CNORM) CNORM = CNI 40 CONTINUE GNORM = (BNORM+SQRT(BNORM**2+4.0*ANORM*CNORM))/(2.0*ANORM) IF (IGUESS.EQ.0) GO TO 70 DO 60 I=1,N DO 50 J=1,N XOLD(I,J) = S(I,J) 50 CONTINUE 60 CONTINUE GO TO 100 C 70 DO 90 I=1,N DO 80 J=1,N XOLD(I,J) = CMPLX(0.0,0.0) 80 CONTINUE XOLD(I,I) = CMPLX(GNORM,0.0) 90 CONTINUE C 100 DO 360 ITS=1,MAXITS IF (ITS.NE.31) GO TO 130 DO 120 I=1,N DO 110 J=1,N XOLD(I,J) = CMPLX(0.0,0.0) 110 CONTINUE XOLD(I,I) = CMPLX(0.0,GNORM) 120 CONTINUE C 130 IF (ITS.NE.61) GO TO 160 DO 150 I=1,N DO 140 J=1,N XOLD(I,J) = C(I,J) 140 CONTINUE 150 CONTINUE C ********** SET UP U AND RIGHT HAND SIDE ********** 160 CALL MULT(NM, N, A, XOLD, U) DO 180 I=1,N DO 170 J=1,N U(I,J) = U(I,J) + B(I,J) 170 CONTINUE 180 CONTINUE CALL MULT(NM, N, U, XOLD, S) DO 200 I=1,N DO 190 J=1,N S(I,J) = S(I,J) + C(I,J) 190 CONTINUE 200 CONTINUE C ********** CHECK FOR CONVERGENCE ********** XNORM = 0.0 FXNORM = 0.0 DO 220 I=1,N XNI = 0.0 FXNI = 0.0 DO 210 J=1,N XNI = XNI + CABS(XOLD(I,J)) FXNI = FXNI + CABS(S(I,J)) 210 CONTINUE IF (XNI.GT.XNORM) XNORM = XNI IF (FXNI.GT.FXNORM) FXNORM = FXNI 220 CONTINUE IF (TOL.LE.0.0) GO TO 230 IF (FXNORM.LT.TOL) GO TO 370 230 TNORM = 8.0*FLOAT(N)*ANORM*XNORM**2 + 5.0*FLOAT(N)*BNORM*XNORM * + CNORM T = 1.0 + FXNORM/TNORM IF (T.EQ.1.0) GO TO 370 IF (ITS.GE.MAXITS) GO TO 400 C ********** UPPER TRIANGULARIZATION ********** IF (ITS.NE.1) GO TO 240 MATS = 1 CALL CQZHES(NM, N, U, A, MATS, Z, S, B, C) 240 DO 260 I=1,N DO 250 J=1,N V(I,J) = A(I,J) L(I,J) = CONJG(XOLD(J,I)) EYE(I,J) = CMPLX(0.0,0.0) 250 CONTINUE EYE(I,I) = CMPLX(1.0,0.0) 260 CONTINUE IF (ITS.EQ.1) GO TO 270 MATS = 2 CALL CQZHES(NM, N, U, V, MATS, Z, S, B, C) 270 CALL CQZIT(NM, N, U, V, 0.0, MATS, Z, S, IERR) IF (IERR.NE.0) GO TO 430 C ********** LOWER TRIANGULARIZATION ********** MATS = 3 CALL CQZHES(NM, N, L, EYE, MATS, R, S, B, C) CALL CQZIT(NM, N, L, EYE, 0.0, MATS, R, S, IERR) IF (IERR.NE.0) GO TO 440 CALL CTRANS(NM, N, L) C ********** UPDATE S WITH R ********** DO 310 I=1,N DO 280 J=1,N TEMP(J) = S(I,J) S(I,J) = CMPLX(0.0,0.0) 280 CONTINUE DO 300 J=1,N DO 290 K=1,N S(I,J) = S(I,J) + TEMP(K)*R(K,J) 290 CONTINUE 300 CONTINUE 310 CONTINUE DO 330 J=1,N DO 320 I=1,N L(I,J) = L(I,J)*EYE(J,J) 320 CONTINUE EYE(J,J) = CMPLX(1.0,0.0) 330 CONTINUE C ********** BACKSOLVE THE TRANSFORMED SYSTEM ********** CALL TRISLV(NM, N, U, V, L, S, TEMP, IERR) IF (IERR.NE.0) GO TO 450 C ********** TRANSLATE BACK TO THE SOLUTION ********** CALL MULT(NM, N, Z, S, L) CALL CTRANS(NM, N, R) CALL MULT(NM, N, L, R, S) DO 350 I=1,N DO 340 J=1,N XOLD(I,J) = XOLD(I,J) - S(I,J) 340 CONTINUE 350 CONTINUE 360 CONTINUE C ********** CONVERGENCE ********** 370 IGUESS = ITS - 1 L(1,1) = CMPLX(FXNORM,0.0) DO 390 I=1,N DO 380 J=1,N S(I,J) = XOLD(I,J) 380 CONTINUE 390 CONTINUE IERR = 0 RETURN C ********** ERROR RETURNS ********** 400 IERR = 1 L(1,1) = CMPLX(FXNORM,0.0) DO 420 I=1,N DO 410 J=1,N S(I,J) = XOLD(I,J) 410 CONTINUE 420 CONTINUE RETURN 430 IERR = 2 RETURN 440 IERR = 3 RETURN 450 IERR = ITS + 10 RETURN 460 IERR = 999 RETURN END SUBROUTINE CTRANS(NM, N, A) CTR 10 C C C SUBROUTINE CTRANS FINDS THE COMPLEX CONJUGATE OF AN INPUT C MATRIX. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF MATRIX A IN THE MAIN PROGRAM. C C N IS THE ORDER OF MATRIX A. C C A IS THE INPUT MATRIX. C C C ON RETURN, C C A CONTAINS ITS CONJUGATE TRANSPOSE. C C INTEGER I, J, N, NM COMPLEX A(NM,N), TEMP COMPLEX CONJG DO 20 I=1,N DO 10 J=I,N TEMP = A(I,J) A(I,J) = CONJG(A(J,I)) A(J,I) = CONJG(TEMP) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE MULT(NM, N, A, B, C) MUL 10 C C C SUBROUTINE MULT MULTIPLIES TWO MATRICES TOGETHER. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF MATRICES A, B AND C IN THE C MAIN PROGRAM. C C N IS THE ORDER OF THE MATRICES. C C A,B ARE THE MATRICES TO BE MULTIPLIED. C C C IS A WORK MATRIX WHICH IS INTERNALLY INITIALIZED TO ZERO. C C C ON RETURN, C C A,B ARE UNALTERED. C C C CONTAINS THE MATRIX PRODUCT A*B. C C INTEGER I, J, K, N, NM COMPLEX A(NM,N), B(NM,N), C(NM,N) COMPLEX CMPLX DO 30 I=1,N DO 20 J=1,N C(I,J) = CMPLX(0.0,0.0) DO 10 K=1,N C(I,J) = C(I,J) + A(I,K)*B(K,J) 10 CONTINUE 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE CQZHES(NM, N, A, B, MATS, Z, F, G, H) CQZ 10 C C C SUBROUTINE CQZHES IS A MODIFICATION OF THE EISPACK SUBROUTINE C QZHES. ALL OPERATIONS ARE PERFORMED IN COMPLEX ARITHMETIC, C AND THE LEFT TRANSFORMATIONS MAY ALSO BE APPLIED TO AUXILIARY C MATRICES F, G AND H. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF THE MATRICES A AND B IN C THE MAIN PROGRAM. C C N IS THE ORDER OF THE MATRICES A AND B. C C A CONTAINS THE MATRIX TO BE REDUCED TO UPPER HESSENBERG C FORM. C C B CONTAINS THE MATRIX TO BE REDUCED TO UPPER TRAINGULAR C FORM. C C MATS IS AN INTEGER INPUT VARIABLE. C C IF MATS = 0, THE ACCUMULATION OF THE TRANSFORMATIONS C IS NOT DESIRED. C C IF MATS = ANY OTHER NUMBER BUT 0, THE TRANSFORMATIONS C ARE ACCUMULATED. C C IF MATS = 1, THE AUXILIARY MATRICES G AND H ARE UPDATED C WITH THE UNITARY MATRIX Q. C C IF MATS = 2, MATRIX B IS ASSUMED UPPER TRIANGULAR. C C IF MATS = 3, THE AUXILIARY MATRIX F IS NOT UPDATED C WITH THE UNITARY MATRIX Q. C C F, G AND H ARE AUXILIARY MATRICES. C C C ON RETURN, C C A IS UPPER HESSENBERG. C C B IS UPPER TRIANGULAR. C C Z CONTAINS THE HISTORY OF THE TRANSFORMATIONS, IF DESIRED. C C F, G AND H ARE UPDATED, IF DESIRED. C C INTEGER I, J, K, L, N, LB, L1, NM, NK1, NM1, NM2 COMPLEX A(NM,N), B(NM,N), Z(NM,N) COMPLEX RR, T, U1, U2, V1, V2, RHO COMPLEX F(NM,N), TF, G(NM,N), TG, H(NM,N), TH COMPLEX CMPLX, CONJG REAL R, S REAL SQRT, CABS INTEGER MATS IF (MATS.EQ.0) GO TO 30 DO 20 I=1,N DO 10 J=1,N Z(I,J) = CMPLX(0.0,0.0) 10 CONTINUE Z(I,I) = CMPLX(1.0,0.0) 20 CONTINUE C ********** REDUCE B TO UPPER TRIANGULAR FORM ********** 30 IF (N.LE.1) GO TO 260 NM1 = N - 1 IF (MATS.EQ.2) GO TO 140 DO 130 L=1,NM1 L1 = L + 1 S = 0.0 DO 40 I=L1,N S = S + CABS(B(I,L)) 40 CONTINUE IF (S.EQ.0.0) GO TO 130 S = S + CABS(B(L,L)) R = 0.0 DO 50 I=L,N B(I,L) = B(I,L)/CMPLX(S,0.0) R = R + CABS(B(I,L))**2 50 CONTINUE R = SQRT(R) RR = CMPLX(R,0.0) IF (CABS(B(L,L)).NE.0.0) RR = (B(L,L)/CABS(B(L,L)))*RR B(L,L) = B(L,L) + RR RHO = CONJG(RR)*B(L,L) DO 80 J=L1,N T = CMPLX(0.0,0.0) DO 60 I=L,N T = T + CONJG(B(I,L))*B(I,J) 60 CONTINUE T = -T/RHO DO 70 I=L,N B(I,J) = B(I,J) + T*B(I,L) 70 CONTINUE 80 CONTINUE DO 110 J=1,N T = CMPLX(0.0,0.0) TF = CMPLX(0.0,0.0) TG = CMPLX(0.0,0.0) TH = CMPLX(0.0,0.0) DO 90 I=L,N T = T + CONJG(B(I,L))*A(I,J) IF (MATS.EQ.3) GO TO 90 TF = TF + CONJG(B(I,L))*F(I,J) IF (MATS.NE.1) GO TO 90 TG = TG + CONJG(B(I,L))*G(I,J) TH = TH + CONJG(B(I,L))*H(I,J) 90 CONTINUE T = -T/RHO TF = -TF/RHO TG = -TG/RHO TH = -TH/RHO DO 100 I=L,N A(I,J) = A(I,J) + T*B(I,L) IF (MATS.EQ.3) GO TO 100 F(I,J) = F(I,J) + TF*B(I,L) IF (MATS.NE.1) GO TO 100 G(I,J) = G(I,J) + TG*B(I,L) H(I,J) = H(I,J) + TH*B(I,L) 100 CONTINUE 110 CONTINUE B(L,L) = -CMPLX(S,0.0)*RR DO 120 I=L1,N B(I,L) = CMPLX(0.0,0.0) 120 CONTINUE 130 CONTINUE C ********** REDUCE A TO UPPER HESSENBERG FORM, WHILE C KEEPING B TRIANGULAR ********** 140 IF (N.EQ.2) GO TO 260 NM2 = N - 2 DO 250 K=1,NM2 NK1 = NM1 - K DO 240 LB=1,NK1 L = N - LB L1 = L + 1 C ********** ZERO A(L+1,K) ********** S = CABS(A(L,K)) + CABS(A(L1,K)) IF (S.EQ.0.0) GO TO 240 U1 = A(L,K)/CMPLX(S,0.0) U2 = A(L1,K)/CMPLX(S,0.0) R = SQRT(CABS(U1)**2+CABS(U2)**2) RR = CMPLX(R,0.0) IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR V1 = -(U1+RR)/RR V2 = -U2/RR U2 = V2/V1 DO 150 J=K,N T = A(L,J) + CONJG(U2)*A(L1,J) A(L,J) = A(L,J) + T*V1 A(L1,J) = A(L1,J) + T*V2 150 CONTINUE A(L1,K) = CMPLX(0.0,0.0) DO 160 J=L,N T = B(L,J) + CONJG(U2)*B(L1,J) B(L,J) = B(L,J) + T*V1 B(L1,J) = B(L1,J) + T*V2 160 CONTINUE IF (MATS.EQ.3) GO TO 180 DO 170 J=1,N TF = F(L,J) + CONJG(U2)*F(L1,J) F(L,J) = F(L,J) + TF*V1 F(L1,J) = F(L1,J) + TF*V2 170 CONTINUE 180 IF (MATS.NE.1) GO TO 200 DO 190 J=1,N TG = G(L,J) + CONJG(U2) + G(L1,J) TH = H(L,J) + CONJG(U2) + H(L1,J) G(L,J) = G(L,J) + TG*V1 H(L,J) = H(L,J) + TH*V1 G(L1,J) = G(L1,J) + TG*V2 H(L1,J) = H(L1,J) + TH*V2 190 CONTINUE C ********** ZERO B(L+1,L) ********** 200 S = CABS(B(L1,L1)) + CABS(B(L1,L)) IF (S.EQ.0.0) GO TO 240 U1 = B(L1,L1)/CMPLX(S,0.0) U2 = B(L1,L)/CMPLX(S,0.0) R = SQRT(CABS(U1)**2+CABS(U2)**2) RR = CMPLX(R,0.0) IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR V1 = -(U1+RR)/RR V2 = -U2/RR U2 = V2/V1 DO 210 I=1,L1 T = B(I,L1) + CONJG(U2)*B(I,L) B(I,L1) = B(I,L1) + T*V1 B(I,L) = B(I,L) + T*V2 210 CONTINUE B(L1,L) = CMPLX(0.0,0.0) DO 220 I=1,N T = A(I,L1) + CONJG(U2)*A(I,L) A(I,L1) = A(I,L1) + T*V1 A(I,L) = A(I,L) + T*V2 220 CONTINUE IF (MATS.EQ.0) GO TO 240 DO 230 I=1,N T = Z(I,L1) + CONJG(U2)*Z(I,L) Z(I,L1) = Z(I,L1) + T*V1 Z(I,L) = Z(I,L) + T*V2 230 CONTINUE 240 CONTINUE 250 CONTINUE 260 RETURN END SUBROUTINE CQZIT(NM, N, A, B, EPS1, MATS, Z, F, IERR) CQZ 10 C C C SUBROUTINE CQZIT IS A MODIFICATION OF THE EISPACK SUBROUTINE C QZIT. ALL OPERATIONS ARE PERFORMED IN COMPLEX ARITHMETIC, C AND THE LEFT TRANSFORMATIONS MAY ALSO BE APPLIED TO AN AUXILIARY C MATRIX F. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF THE MATRICES A AND B IN C THE MAIN PROGRAM. C C N IS THE ORDER OF THE MATRICES A AND B. C C A CONTAINS AN UPPER HESSENBERG MATRIX FROM CQZHES. C C B CONTAINS AN UPPER TRIANGULAR MATRIX FROM CQZHES. C C EPS1 IS A REAL NUMBER DEFINING THE TOLERANCE USED TO DETERMINE C NEGLIGIBLE ELEMENTS OF A AND B IN THE COURSE OF THE ALG- C ORITHM. AN ELEMENT OF EITHER MATRIX WILL BE CONSIDERED C NEGLIGIBLE AND RESET TO ZERO IF IT IS NOT LARGER THAN THE C PRODUCT OF EPS1 AND THE NORM OF THE MATRIX. IF EPS1.LE.0, C RELATIVE MACHINE PRECISION WILL BE COMPUTED AND C USED INSTEAD. C C MATS IS AN INTEGER INPUT VARIABLE. IT IS SET PRIOR C TO THE CALL TO CQZHES. C C F CONTAINS AN AUXILIARY MATRIX. C C C ON RETURN, C C A IS UPPER TRIANGULAR. C C B IS UPPER TRIANGULAR. C C Z CONTAINS THE HISTORY OF THE TRANSFORMATIONS, IF DESIRED. C C F CONTAINS THE AUXILIARY MATRIX, UPDATED IF DESIRED. C C IERR IS AN INTEGER ERROR RETURN WHICH INDICATES FAILURE C OF THE QZ ALGORITHM TO REDUCE A SUBDIAGONAL ELEMENT C TO ZERO AFTER 50 ITERATIONS. C INTEGER I, J, K, L, N, EN, JJ, K1, K2, LD, LL, L1, NA, NM, ISH, * ITS, KM1, LM1 INTEGER ENM2, IERR, LOR1, ENORN COMPLEX A(NM,N), B(NM,N), Z(NM,N) COMPLEX A11, A21, A33, A34, A43, A44, B11, B22, B33, B34, B44 COMPLEX A1, A2, U1, U2, V1, V2, T, RR, SH, SS COMPLEX F(NM,N), TF COMPLEX CMPLX, CONJG, CSQRT REAL EPS1, EPSA, EPSB, ANORM, BNORM, ANI, BNI, SRELPR, R, S REAL SQRT, CABS INTEGER MAX0, MIN0, MATS IERR = 0 C ********** COMPUTE EPSA, EPSB ********** ANORM = 0.0 BNORM = 0.0 DO 20 I=1,N ANI = 0.0 IF (I.NE.1) ANI = CABS(A(I,I-1)) BNI = 0.0 DO 10 J=I,N ANI = ANI + CABS(A(I,J)) BNI = BNI + CABS(B(I,J)) 10 CONTINUE IF (ANI.GT.ANORM) ANORM = ANI IF (BNI.GT.BNORM) BNORM = BNI 20 CONTINUE IF (ANORM.EQ.0.0) ANORM = 1.0 IF (BNORM.EQ.0.0) BNORM = 1.0 SRELPR = EPS1 IF (SRELPR.GT.0.0) GO TO 40 C ********** COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO ********** SRELPR = 1.0 30 SRELPR = SRELPR/2.0 IF (1.0+SRELPR.GT.1.0) GO TO 30 SRELPR = SRELPR + SRELPR 40 EPSA = SRELPR*ANORM EPSB = SRELPR*BNORM C ********** REDUCE A TO TRIANGULAR FORM, WHILE C KEEPING B TRIANGULAR ********** LOR1 = 1 ENORN = N EN = N C ********** BEGIN QZ STEP ********** 50 IF (EN.LE.1) GO TO 220 IF (MATS.EQ.0) ENORN = EN ITS = 0 NA = EN - 1 ENM2 = NA - 1 60 ISH = 1 C ********** CHECK FOR CONVERGENCE OR REDUCIBILITY ********** DO 70 LL=1,EN LM1 = EN - LL L = LM1 + 1 IF (L.EQ.1) GO TO 90 IF (CABS(A(L,LM1)).LE.EPSA) GO TO 80 70 CONTINUE 80 A(L,LM1) = CMPLX(0.0,0.0) IF (L.LT.NA) GO TO 90 C ********** 1-BY-1 BLOCK ISOLATED ********** EN = LM1 GO TO 50 C ********** CHECK FOR SMALL TOP OF B ********** 90 LD = L L1 = L + 1 B11 = B(L,L) IF (CABS(B11).GT.EPSB) GO TO 120 B(L,L) = CMPLX(0.0,0.0) S = CABS(A(L,L)) + CABS(A(L1,L)) U1 = A(L,L)/CMPLX(S,0.0) U2 = A(L1,L)/CMPLX(S,0.0) R = SQRT(CABS(U1)**2+CABS(U2)**2) RR = CMPLX(R,0.0) IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR V1 = -(U1+RR)/RR V2 = -U2/RR U2 = V2/V1 DO 110 J=L,ENORN T = A(L,J) + CONJG(U2)*A(L1,J) A(L,J) = A(L,J) + T*V1 A(L1,J) = A(L1,J) + T*V2 T = B(L,J) + CONJG(U2)*B(L1,J) B(L,J) = B(L,J) + T*V1 B(L1,J) = B(L1,J) + T*V2 IF (MATS.EQ.3) GO TO 110 DO 100 JJ=1,N TF = F(L,JJ) + CONJG(U2)*F(L1,JJ) F(L,JJ) = F(L,JJ) + TF*V1 F(L1,JJ) = F(L1,JJ) + TF*V2 100 CONTINUE 110 CONTINUE IF (L.NE.1) A(L,LM1) = -A(L,LM1) LM1 = L L = L1 GO TO 80 120 A11 = A(L,L)/B11 A21 = A(L1,L)/B11 C ********** ITERATION STRATEGY ********** IF (ITS.EQ.50) GO TO 210 C ********** DETERMINE SHIFT ********** B22 = B(L1,L1) IF (CABS(B22).LT.EPSB) B22 = CMPLX(EPSB,0.0) B33 = B(NA,NA) IF (CABS(B33).LT.EPSB) B33 = CMPLX(EPSB,0.0) B44 = B(EN,EN) IF (CABS(B44).LT.EPSB) B44 = CMPLX(EPSB,0.0) A33 = A(NA,NA)/B33 A34 = A(NA,EN)/B44 A43 = A(EN,NA)/B33 A44 = A(EN,EN)/B44 B34 = B(NA,EN)/B44 T = CMPLX(0.5,0.0)*(A43*B34-A33-A44) RR = T*T + A34*A43 - A33*A44 C ********** DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ********** RR = CSQRT(RR) SH = -T + RR SS = -T - RR IF (CABS(SS-A44).LT.CABS(SH-A44)) SH = SS A1 = A11 - SH A2 = A21 IF (L.NE.LD) A(L,LM1) = -A(L,LM1) IF (ITS.NE.10) GO TO 130 A1 = CMPLX(1.0,0.0) A2 = CMPLX(1.1605,0.0) 130 ITS = ITS + 1 IF (MATS.EQ.0) LOR1 = LD C ********** MAIN LOOP ********** DO 200 K=L,NA K1 = K + 1 K2 = K + 2 KM1 = MAX0(K-1,L) LL = MIN0(EN,K1+ISH) C ********** ZERO A(K+1,K-1) ********** IF (K.EQ.L) GO TO 140 A1 = A(K,KM1) A2 = A(K1,KM1) 140 S = CABS(A1) + CABS(A2) IF (S.EQ.0.0) GO TO 60 U1 = A1/CMPLX(S,0.0) U2 = A2/CMPLX(S,0.0) R = SQRT(CABS(U1)**2+CABS(U2)**2) RR = CMPLX(R,0.0) IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR V1 = -(U1+RR)/RR V2 = -U2/RR U2 = V2/V1 DO 150 J=KM1,ENORN T = A(K,J) + CONJG(U2)*A(K1,J) A(K,J) = A(K,J) + T*V1 A(K1,J) = A(K1,J) + T*V2 T = B(K,J) + CONJG(U2)*B(K1,J) B(K,J) = B(K,J) + T*V1 B(K1,J) = B(K1,J) + T*V2 150 CONTINUE IF (K.NE.L) A(K1,KM1) = CMPLX(0.0,0.0) IF (MATS.EQ.3) GO TO 170 DO 160 J=1,N TF = F(K,J) + CONJG(U2)*F(K1,J) F(K,J) = F(K,J) + TF*V1 F(K1,J) = F(K1,J) + TF*V2 160 CONTINUE C ********** ZERO B(K+1,K) ********** 170 S = CABS(B(K1,K1)) + CABS(B(K1,K)) IF (S.EQ.0.0) GO TO 200 U1 = B(K1,K1)/CMPLX(S,0.0) U2 = B(K1,K)/CMPLX(S,0.0) R = SQRT(CABS(U1)**2+CABS(U2)**2) RR = CMPLX(R,0.0) IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR V1 = -(U1+RR)/RR V2 = -U2/RR U2 = V2/V1 DO 180 I=LOR1,LL T = A(I,K1) + CONJG(U2)*A(I,K) A(I,K1) = A(I,K1) + T*V1 A(I,K) = A(I,K) + T*V2 T = B(I,K1) + CONJG(U2)*B(I,K) B(I,K1) = B(I,K1) + T*V1 B(I,K) = B(I,K) + T*V2 180 CONTINUE B(K1,K) = CMPLX(0.0,0.0) IF (MATS.EQ.0) GO TO 200 DO 190 I=1,N T = Z(I,K1) + CONJG(U2)*Z(I,K) Z(I,K1) = Z(I,K1) + T*V1 Z(I,K) = Z(I,K) + T*V2 190 CONTINUE 200 CONTINUE C ********** END QZ STEP ********** GO TO 60 C ********** SET ERROR -- NEITHER BOTTOM SUBDIAGONAL ELEMENT C HAS BECOME NEGLIGIBLE AFTER 50 ITERATIONS ********** 210 IERR = EN C ********** SAVE EPSB FOR USE BY CQZVEC ********** 220 IF (N.GT.1) B(N,1) = CMPLX(EPSB,0.0) RETURN END SUBROUTINE TRISLV(NM, N, U, V, L, F, TEMP, IERR) TRI 10 C C C SUBROUTINE TRISLV BACKSOLVES A SYSTEM OF THE FORM C UY + VYL = F, WHERE U AND V ARE UPPER TRIANGULAR, AND C L IS LOWER TRIANGULAR. C C C ON ENTRY, C C NM IS THE LEADING DIMENSION OF THE MATRICES U, V, L AND F IN C THE MAIN PROGRAM. C C N IS THE ORDER OF THE MATRICES U, V, L AND F. C C U CONTAINS AN UPPER TRIANGULAR MATRIX. IT IS THE LEFT C COEFFICIENT OF Y IN THE FIRST TERM IN UY + VYL. C C V CONTAINS AN UPPER TRIANGULAR MATRIX. IT IS THE LEFT C COEFFICIENT OF Y IN THE SECOND TERM OF UY + VYL. C C L CONTAINS A LOWER TRIANGULAR MATRIX. IT IS THE RIGHT C COEFFICIENT OF Y IN THE SECOND TERM OF UY + VYL. C C F CONTAINS THE RIGHT HAND SIDE OF UY + VYL = F. C C TEMP CONTAINS A WORK VECTOR OF LENGTH AT LEAST N. C C ON RETURN, C C F CONTAINS THE SOLUTION Y. C C IERR IS AN ERROR RETURN DESIGNATING INCONSISTENCY OF THE C ORIGINAL SYSTEM. C IERR.EQ.0 FOR A NORMAL RETURN. C IERR.EQ.1 IF THE TRIANGULAR SYSTEM IS INCONSISTENT. C C INTEGER I, IERR, J, JP1, K, KK, KM1, M, N, NM, NM1 COMPLEX U(NM,N), V(NM,N), L(NM,N), F(NM,N) COMPLEX TEMP(N) COMPLEX CMPLX, DENOM, SUM REAL CABS, S, T IERR = 0 NM1 = N - 1 DO 120 KK=1,N C ********** BACKSUBSTITUTE FOR ROW K. ********** K = N - KK + 1 IF (CABS(F(K,N)).NE.0.0) GO TO 10 F(K,N) = CMPLX(0.0,0.0) GO TO 30 10 DENOM = U(K,K) + V(K,K)*L(N,N) S = CABS(DENOM) T = 1.0 + S/CABS(F(K,N)) IF (T.GT.1.0) GO TO 20 IERR = 1 RETURN 20 F(K,N) = F(K,N)/DENOM IF (N.EQ.1) RETURN 30 DO 70 I=1,NM1 J = N - I JP1 = J + 1 SUM = CMPLX(0.0,0.0) DO 40 M=JP1,N SUM = SUM + F(K,M)*L(M,J) 40 CONTINUE SUM = F(K,J) - V(K,K)*SUM IF (CABS(SUM).NE.0.0) GO TO 50 F(K,J) = CMPLX(0.0,0.0) GO TO 70 50 DENOM = U(K,K) + V(K,K)*L(J,J) S = CABS(DENOM) T = 1.0 + S/CABS(SUM) IF (T.GT.1.0) GO TO 60 IERR = 1 RETURN 60 F(K,J) = SUM/DENOM 70 CONTINUE C ********** FORM TEMP = YK-TRANS*L. ********** IF (K.EQ.1) RETURN KM1 = K - 1 DO 90 I=1,N TEMP(I) = CMPLX(0.0,0.0) DO 80 J=1,N TEMP(I) = TEMP(I) + F(K,J)*L(J,I) 80 CONTINUE 90 CONTINUE C ********** PREPARE F' WHICH IS (K-1) BY N. ********** DO 110 I=1,KM1 DO 100 J=1,N F(I,J) = F(I,J) - U(I,K)*F(K,J) - V(I,K)*TEMP(J) 100 CONTINUE 110 CONTINUE 120 CONTINUE RETURN END