C ALGORITHM 586, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.8, NO. 3, C SEP., 1982, P. 302. SUBROUTINE JCG(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, JCG 0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE JCG (JACOBI CONJUGATE GRADIENT) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C ... FUNCTION: C C THIS SUBROUTINE, JCG, DRIVES THE JACOBI CONJUGATE C GRADIENT ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. JACOBI CONJUGATE C GRADIENT NEEDS THIS TO BE IN LENGTH AT LEAST C 4*N + 2*ITMAX, IF ISYM = 0 (SYMMETRIC STORAGE) C 4*N + 4*ITMAX, IF ISYM = 1 (NONSYMMETRIC STORAGE) C HERE ITMAX = IPARM(1) AND ISYM = IPARM(5) C (ITMAX IS THE MAXIMUM ALLOWABLE NUMBER OF ITERATIONS) C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... JCG SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, CHGCON, DETERM, DFAULT, ECHALL, C ECHOUT, EIGVNS, EIGVSS, EQRT1S, ITERM, ITICK, C ITJCG, ITOCK, IVFILL, PARCON, PERMAT, C PERROR, PERVEC, PJAC, PMULT, PRBNDX, C PSTOP, QSORT, SAXPY, SBELM, SCAL, SCOPY, C SDOT, SUM3, UNSCAL, VEVMW, VFILL, VOUT, C WEVMW, ZBRENT C SYSTEM ABS, ALOG10, AMAX0, AMAX1, MOD, SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, IB4, IB5, IDGTS, IER, IERPER, ITIM1, ITIM2, A ITMAX1, JTIM1, JTIM2, LOOP, N, NB, NDUMMY, N3 REAL DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,40HBEGINNING OF ITPACK SOLUTION MODULE JCG ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,1) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 25H IN ITPACK ROUTINE JCG B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF(N .GT. 0 ) GO TO 15 IER = 11 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 25 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 25 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) GO TO 160 C C ... INITIALIZE WKSP BASE ADDRESSES. C 25 IB1 = 1 IB2 = IB1 + N IB3 = IB2 + N IB4 = IB3 + N IB5 = IB4 + N IPARM(8) = 4*N + 2*ITMAX IF(ISYM .NE. 0) IPARM(8) = IPARM(8) + 2*ITMAX IF (NW .GE. IPARM(8)) GO TO 40 IER = 12 IF (LEVEL .GE. 0) WRITE(NOUT,30) NW,IPARM(8) 30 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) GO TO 160 C C ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED C 40 NB = IPARM(9) IF(NB .LT. 0) GO TO 50 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 43 IF (LEVEL .GE. 0) WRITE(NOUT,42) IER,NB 42 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 C C ... PERMUTE MATRIX AND RHS C 43 IF (LEVEL .GE. 2) WRITE(NOUT,44) NB 44 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT(N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 46 IF (LEVEL .GE. 0) WRITE (NOUT,45) IER 45 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 46 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 50 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 64 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 64 IF (LEVEL .LE. 2) GO TO 70 WRITE(NOUT,66) 66 FORMAT(///1X,35HIN THE FOLLOWING, RHO AND GAMMA ARE A ,24H ACCELERATION PARAMETERS ) IF (ADAPT) WRITE(NOUT,67) 67 FORMAT(1X,48HCME IS THE ESTIMATE OF THE LARGEST EIGENVALUE OF A ,18H THE JACOBI MATRIX ) 70 IF(IPARM(11).NE.0)GO TO 75 ITIM1 = ITICK(NDUMMY) C C ... COMPUTE INITIAL PSEUDO-RESIDUAL C 75 CONTINUE CALL SCOPY(N,RHS,1,WKSP(IB2),1) CALL PJAC(N,IA,JA,A,U,WKSP(IB2)) CALL VEVMW(N,WKSP(IB2),U) C C ... ITERATION SEQUENCE C ITMAX1 = ITMAX + 1 DO 95 LOOP=1,ITMAX1 IN = LOOP -1 IF (MOD(IN,2) .EQ. 1) GO TO 85 C C ... CODE FOR THE EVEN ITERATIONS. C C U = U(IN) WKSP(IB2) = DEL(IN) C WKSP(IB1) = U(IN-1) WKSP(IB3) = DEL(IN-1) C CALL ITJCG(N,IA,JA,A,U,WKSP(IB1),WKSP(IB2),WKSP(IB3), A WKSP(IB4),WKSP(IB5)) C IF (HALT) GO TO 132 GO TO 95 C C ... CODE FOR THE ODD ITERATIONS. C C U = U(IN-1) WKSP(IB2) = DEL(IN-1) C WKSP(IB1) = U(IN) WKSP(IB3) = DEL(IN) C 85 CALL ITJCG(N,IA,JA,A,WKSP(IB1),U,WKSP(IB3),WKSP(IB2), A WKSP(IB4),WKSP(IB5)) C IF (HALT) GO TO 132 95 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF(IPARM(11).NE.0)GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IER = 13 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 25H IN ITPACK ROUTINE JCG B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF(IPARM(11).NE.0) GO TO 134 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 134 IF (LEVEL .GE. 1) WRITE(NOUT,136) IN 136 FORMAT(/1X,22HJCG HAS CONVERGED IN ,I5,11H ITERATIONS ) C C ... PUT SOLUTION INTO U IF NOT ALREADY THERE. C 140 CONTINUE IF (MOD(IN,2) .EQ. 1) CALL SCOPY(N,WKSP(IB1),1,U,1) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .LT. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(IB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,142) IERPER 142 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IPARM(8) = IPARM(8) - 2 * (ITMAX - IN) IF(IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IF(ISYM .NE. 0) IPARM(8) = IPARM(8) - 2*(ITMAX - IN) IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE JSI(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, JSI 0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE JSI (JACOBI SEMI-ITERATIVE) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C ... FUNCTION: C C THIS SUBROUTINE, JSI, DRIVES THE JACOBI SEMI- C ITERATION ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. JACOBI SI C GRADIENT NEEDS THIS TO BE IN LENGTH AT LEAST C 2*N C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... JSI SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, CHEBY, CHGSI, CHGSME, DFAULT, ECHALL, C ECHOUT, ITERM, ITICK, ITJSI, ITOCK, IVFILL, PAR C PERMAT, PERROR, PERVEC, PJAC, PMULT, PRBNDX, C PSTOP, PVTBV, QSORT, SAXPY, SBELM, SCAL, C SCOPY, SDOT, SUM3, TSTCHG, UNSCAL, VEVMW, C VFILL, VOUT, WEVMW C SYSTEM ABS, ALOG10, AMAX0, AMAX1, FLOAT, MOD, SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, ICNT, IDGTS, IER, IERPER, ITIM1, ITIM2, A ITMAX1, JTIM1, JTIM2, LOOP, N, NB, NDUMMY, N3 REAL DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,40HBEGINNING OF ITPACK SOLUTION MODULE JSI ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,2) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 25H IN ITPACK ROUTINE JSI B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF(N .GT. 0 ) GO TO 15 IER = 21 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 25 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 25 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) GO TO 160 C C ... INITIALIZE WKSP BASE ADDRESSES. C 25 IB1 = 1 IB2 = IB1 + N IB3 = IB2 + N IPARM(8) = 2*N IF (NW .GE. IPARM(8)) GO TO 40 IER = 22 IF (LEVEL .GE. 0) WRITE(NOUT,27) NW,IPARM(8) 27 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) GO TO 160 C C ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED C 40 NB = IPARM(9) IF(NB .LT. 0) GO TO 50 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 43 IF (LEVEL .GE. 0) WRITE(NOUT,42) IER,NB 42 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 C C ... PERMUTE MATRIX AND RHS C 43 IF (LEVEL .GE. 2) WRITE(NOUT,44) NB 44 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 46 IF (LEVEL .GE. 0) WRITE (NOUT,45) IER 45 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 46 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 50 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 61 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 61 IF (LEVEL .LE. 2) GO TO 64 WRITE (NOUT,63) 63 FORMAT(///1X,35HIN THE FOLLOWING, RHO AND GAMMA ARE A ,24H ACCELERATION PARAMETERS ) 64 IF(IPARM(11).NE.0)GO TO 68 ITIM1 = ITICK(NDUMMY) C C ... ITERATION SEQUENCE C 68 ITMAX1 = ITMAX + 1 DO 95 LOOP=1,ITMAX1 IN = LOOP - 1 IF (MOD(IN,2) .EQ. 1) GO TO 85 C C ... CODE FOR THE EVEN ITERATIONS. C C U = U(IN) C WKSP(IB1) = U(IN-1) C CALL ITJSI(N,IA,JA,A,RHS,U,WKSP(IB1),WKSP(IB2),ICNT) C IF (HALT) GO TO 132 GO TO 95 C C ... CODE FOR THE ODD ITERATIONS. C C U = U(IN-1) C WKSP(IB1) = U(IN) C 85 CALL ITJSI(N,IA,JA,A,RHS,WKSP(IB1),U,WKSP(IB2),ICNT) C IF (HALT) GO TO 132 95 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF(IPARM(11).NE.0)GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IER = 23 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 25H IN ITPACK ROUTINE JSI B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF(IPARM(11).NE.0)GO TO 136 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 136 IF (LEVEL .GE. 1) WRITE(NOUT,138) IN 138 FORMAT(/1X,22HJSI HAS CONVERGED IN ,I5,11H ITERATIONS ) C C ... PUT SOLUTION INTO U IF NOT ALREADY THERE. C 140 CONTINUE IF (MOD(IN,2) .EQ. 1) CALL SCOPY(N,WKSP(IB1),1,U,1) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .LT. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(IB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,142) IERPER 142 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE JSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IF(IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE SOR(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, SOR 0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE SOR (SUCCESSIVE OVERRELATION) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C ... FUNCTION: C C THIS SUBROUTINE, SOR, DRIVES THE SUCCESSIVE C OVERRELAXATION ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. SOR NEEDS THIS C TO BE IN LENGTH AT LEAST N C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... SOR SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, DFAULT, ECHALL, ECHOUT, IPSTR, ITERM, C ITICK, ITOCK, ITSOR, IVFILL, PERMAT, PERROR, C PERVEC, PFSOR1, PMULT, PRBNDX, PSTOP, QSORT, C SBELM, SCAL, SCOPY, SDOT, TAU, UNSCAL, VFILL, C VOUT, WEVMW C SYSTEM ABS, ALOG10, AMAX0, AMAX1, FLOAT, IABS, MOD, C SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, IDGTS, IER, IERPER, ITIM1, ITIM2, A ITMAX1, JTIM1, JTIM2, LOOP, N, NB, NDUMMY, N3 REAL DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,40HBEGINNING OF ITPACK SOLUTION MODULE SOR ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,3) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 25H IN ITPACK ROUTINE SOR B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF(N .GT. 0 ) GO TO 15 IER = 31 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 25 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 25 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) GO TO 160 C C ... INITIALIZE WKSP BASE ADDRESSES. C 25 IB1 = 1 IB2 = IB1 + N IB3 = IB2 + N IPARM(8) = N IF (NW .GE. IPARM(8)) GO TO 40 IER = 32 IF (LEVEL .GE. 0) WRITE(NOUT,30) NW,IPARM(8) 30 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) GO TO 160 C C ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED C 40 NB = IPARM(9) IF(NB .LT. 0) GO TO 50 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 43 IF (LEVEL .GE. 0) WRITE(NOUT,42) IER,NB 42 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 C C ... PERMUTE MATRIX AND RHS C 43 IF (LEVEL .GE. 2) WRITE(NOUT,44) NB 44 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 46 IF (LEVEL .GE. 0) WRITE (NOUT,45) IER 45 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 46 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 50 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 61 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 61 IF (LEVEL .LE. 2) GO TO 65 IF (ADAPT) WRITE (NOUT,62) 62 FORMAT(///1X,48HCME IS THE ESTIMATE OF THE LARGEST EIGENVALUE OF A ,18H THE JACOBI MATRIX ) WRITE(NOUT,64) 64 FORMAT(1X,30HOMEGA IS THE RELAXATION FACTOR ) 65 IF (IPARM(11).NE.0)GO TO 68 ITIM1 = ITICK(NDUMMY) C C ... ITERATION SEQUENCE C 68 ITMAX1 = ITMAX + 1 DO 95 LOOP=1,ITMAX1 IN = LOOP - 1 C C ... CODE FOR ONE ITERATION. C C U = U(IN) C CALL ITSOR(N,IA,JA,A,RHS,U,WKSP(IB1)) C IF (HALT) GO TO 132 95 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF(IPARM(11).NE.0)GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 25H IN ITPACK ROUTINE SOR B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IER = 33 IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF(IPARM(11).NE.0)GO TO 134 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 134 IF (LEVEL .GE. 1) WRITE(NOUT,136) IN 136 FORMAT(/1X,22HSOR HAS CONVERGED IN ,I5,11H ITERATIONS ) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C 140 CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .LT. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(IB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,142) IERPER 142 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 35H CALLED FROM ITPACK ROUTINE SOR B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IF(IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(5) = OMEGA RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE SSORCG(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, SSOR0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE SSORCG (SYMMETRIC SUCCESSIVE OVER- C RELAXATION CONJUGATE GRADIENT) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C ... FUNCTION: C C THIS SUBROUTINE, SSORCG, DRIVES THE SYMMETRIC SOR-CG C ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. SSOR-CG C NEEDS TO BE IN LENGTH AT LEAST C 6*N + 2*ITMAX, IF IPARM(5)=0 (SYMMETRIC STORAGE) C 6*N + 4*ITMAX, IF IPARM(5)=1 (NONSYMMETRIC STORAGE) C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. IF C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... SSORCG SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, CHGCON, DETERM, DFAULT, ECHALL, C ECHOUT, EIGVNS, EIGVSS, EQRT1S, ITERM, ITICK, C ITOCK, ITSRCG, IVFILL, OMEG, OMGCHG, OMGSTR, C PARCON, PBETA, PBSOR, PERMAT, PERROR, C PERVEC, PFSOR, PMULT, PRBNDX, PSTOP, PVTBV, C QSORT, SBELM, SCAL, SCOPY, SDOT, SUM3, C UNSCAL, VEVMW, VEVPW, VFILL, VOUT, WEVMW, C ZBRENT C SYSTEM ABS, ALOG, ALOG10, AMAX0, AMAX1, AMIN1, C MOD, SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, IB4, IB5, IB6, IB7, IDGTS, IER, IERPER, A ITIM1, ITIM2, ITMAX1, JTIM1, JTIM2, LOOP, N, A NB, NDUMMY, N3 REAL BETNEW, DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,43HBEGINNING OF ITPACK SOLUTION MODULE SSORCG ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,4) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 28H IN ITPACK ROUTINE SSORCG B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF(N .GT. 0 ) GO TO 15 IER = 41 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 20 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 20 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) GO TO 160 C C ... INITIALIZE WKSP BASE ADDRESSES. C 20 IB1 = 1 IB2 = IB1 + N IB3 = IB2 + N IB4 = IB3 + N IB5 = IB4 + N IB6 = IB5 + N IB7 = IB6 + N IPARM(8) = 6*N + 2*ITMAX IF(ISYM .NE. 0) IPARM(8) = IPARM(8) + 2*ITMAX IF (NW .GE. IPARM(8)) GO TO 40 IER = 42 IF (LEVEL .GE. 0) WRITE(NOUT,30) NW,IPARM(8) 30 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) GO TO 160 C C ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED C 40 NB = IPARM(9) IF(NB .LT. 0) GO TO 50 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 43 IF (LEVEL .GE. 0) WRITE(NOUT,42) IER,NB 42 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 C C ... PERMUTE MATRIX AND RHS C 43 IF (LEVEL .GE. 2) WRITE(NOUT,44) NB 44 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 46 IF (LEVEL .GE. 0) WRITE (NOUT,45) IER 45 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 46 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 50 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 61 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 61 IF (LEVEL .LE. 2) GO TO 65 WRITE(NOUT,62) 62 FORMAT(///1X,35HIN THE FOLLOWING, RHO AND GAMMA ARE A ,24H ACCELERATION PARAMETERS ) WRITE(NOUT,63) 63 FORMAT(1X,42HS-PRIME IS AN INITIAL ESTIMATE FOR NEW CME ) 65 IF (IPARM(11).NE.0)GO TO 70 ITIM1 = ITICK(NDUMMY) C C ... SPECIAL PROCEDURE FOR FULLY ADAPTIVE CASE. C 70 CONTINUE IF (.NOT.ADAPT) GO TO 80 IF (.NOT.BETADT) GO TO 75 CALL VFILL(N,WKSP(IB1),1.E0) BETNEW = PBETA(N,IA,JA,A,WKSP(IB1),WKSP(IB2),WKSP(IB3) ) A / FLOAT(N) BETAB = AMAX1 ( BETAB, .25E0, BETNEW) 75 CALL OMEG(0.E0,1) IS = 0 C C ... INITIALIZE FORWARD PSEUDO-RESIDUAL C 80 CALL SCOPY(N,RHS,1,WKSP(IB1),1) CALL SCOPY(N,U,1,WKSP(IB2),1) CALL PFSOR(N,IA,JA,A,WKSP(IB2),WKSP(IB1)) CALL VEVMW(N,WKSP(IB2),U) C C ... ITERATION SEQUENCE C ITMAX1 = ITMAX + 1 DO 95 LOOP=1,ITMAX1 IN = LOOP - 1 IF (MOD(IN,2) .EQ. 1) GO TO 85 C C ... CODE FOR THE EVEN ITERATIONS. C C U = U(IN) WKSP(IB2) = C(IN) C WKSP(IB1) = U(IN-1) WKSP(IB3) = C(IN-1) C CALL ITSRCG(N,IA,JA,A,RHS,U,WKSP(IB1),WKSP(IB2),WKSP(IB3), A WKSP(IB4),WKSP(IB5),WKSP(IB6),WKSP(IB7)) C IF (HALT) GO TO 132 GO TO 95 C C ... CODE FOR THE ODD ITERATIONS. C C U = U(IN-1) WKSP(IB2) = C(IN-1) C WKSP(IB1) = U(IN) WKSP(IB3) =C(IN) C 85 CALL ITSRCG(N,IA,JA,A,RHS,WKSP(IB1),U,WKSP(IB3),WKSP(IB2), A WKSP(IB4),WKSP(IB5),WKSP(IB6),WKSP(IB7)) C IF (HALT) GO TO 132 95 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF(IPARM(11).NE.0)GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 28H IN ITPACK ROUTINE SSORCG B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IER = 43 IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF(IPARM(11).NE.0)GO TO 134 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 134 IF (LEVEL .GE. 1) WRITE(NOUT,136) IN 136 FORMAT(/1X,25HSSORCG HAS CONVERGED IN ,I5, A 11H ITERATIONS ) C C ... PUT SOLUTION INTO U IF NOT ALREADY THERE. C 140 CONTINUE IF (MOD(IN,2) .EQ. 1) CALL SCOPY(N,WKSP(IB1),1,U,1) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .LT. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(IB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,142) IERPER 142 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IF(IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IPARM(8) = IPARM(8) - 2 * (ITMAX - IN) IF(ISYM .NE. 0) IPARM(8) = IPARM(8) - 2*(ITMAX - IN) IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(5) = OMEGA RPARM(6) = SPECR RPARM(7) = BETAB RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE SSORSI(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, SSOR0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE SSORSI (SYMMETRIC SUCCESSIVE RELAX- C ATION SEMI-ITERATION) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C ... FUNCTION: C C THIS SUBROUTINE, SSORSI, DRIVES THE SYMMETRIC SOR-SI C ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. SSORSI C NEEDS THIS TO BE IN LENGTH AT LEAST 5*N C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. IF C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... SSORSI SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, CHEBY, CHGSI, DFAULT, ECHALL, ECHOUT, C ITERM, ITICK, ITOCK, ITSRSI, IVFILL, OMEG, C OMGSTR, PARSI, PBETA, PERMAT, PERROR, C PERVEC, PFSOR, PMULT, PRBNDX, PSSOR1, C PSTOP, PVTBV, QSORT, SBELM, SCAL, SCOPY, C SDOT, SUM3, TSTCHG, UNSCAL, VEVPW, VFILL, C VOUT, WEVMW C SYSTEM ABS, ALOG, ALOG10, AMAX0, AMAX1, FLOAT, C MOD, SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, IB4, IB5, IDGTS, IER, IERPER, ITIM1, A ITIM2, ITMAX1, JTIM1, JTIM2, LOOP, N, NB, NDUMMY, N3 REAL BETNEW, DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,43HBEGINNING OF ITPACK SOLUTION MODULE SSORSI ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,5) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 28H IN ITPACK ROUTINE SSORSI B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF(N .GT. 0 ) GO TO 15 IER = 51 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 20 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 20 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) GO TO 160 C C ... INITIALIZE WKSP BASE ADDRESSES. C 20 IB1 = 1 IB2 = IB1 + N IB3 = IB2 + N IB4 = IB3 + N IB5 = IB4 + N IPARM(8) = 5*N IF (NW .GE. IPARM(8)) GO TO 40 IER = 52 IF (LEVEL .GE. 0) WRITE(NOUT,30) NW,IPARM(8) 30 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) C C ... PERMUTE TO RED-BLACK SYSTEM IF REQUESTED C 40 NB = IPARM(9) IF(NB .LT. 0) GO TO 50 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 43 IF (LEVEL .GE. 0) WRITE(NOUT,42) IER,NB 42 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 C C ... PERMUTE MATRIX AND RHS C 43 IF (LEVEL .GE. 2) WRITE(NOUT,44) NB 44 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 46 IF (LEVEL .GE. 0) WRITE (NOUT,45) IER 45 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 46 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 50 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 61 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 61 IF (LEVEL .LE. 2) GO TO 65 WRITE(NOUT,62) 62 FORMAT(///1X,35HIN THE FOLLOWING, RHO AND GAMMA ARE A ,24H ACCELERATION PARAMETERS ) 65 IF (IPARM(11).NE.0)GO TO 70 ITIM1 = ITICK(NDUMMY) C C ... SPECIAL PROCEDURE FOR FULLY ADAPTIVE CASE. C 70 CONTINUE IF (.NOT.ADAPT) GO TO 80 IF (.NOT.BETADT) GO TO 75 CALL VFILL(N,WKSP(IB1),1.E0) BETNEW = PBETA(N,IA,JA,A,WKSP(IB1),WKSP(IB2),WKSP(IB3) ) A / FLOAT(N) BETAB = AMAX1 ( BETAB, .25E0, BETNEW) 75 CALL OMEG(0.E0,1) IS = 0 C C ... ITERATION SEQUENCE C 80 ITMAX1 = ITMAX + 1 DO 95 LOOP=1,ITMAX1 IN = LOOP - 1 IF (MOD(IN,2) .EQ. 1) GO TO 85 C C ... CODE FOR THE EVEN ITERATIONS. C C U = U(IN) C WKSP(IB1) = U(IN-1) C CALL ITSRSI(N,IA,JA,A,RHS,U,WKSP(IB1),WKSP(IB2),WKSP(IB3), A WKSP(IB4),WKSP(IB5)) C IF (HALT) GO TO 132 GO TO 95 C C ... CODE FOR THE ODD ITERATIONS. C C U = U(IN-1) C WKSP(IB1) = U(IN) C 85 CALL ITSRSI(N,IA,JA,A,RHS,WKSP(IB1),U,WKSP(IB2),WKSP(IB3), A WKSP(IB4),WKSP(IB5)) C IF (HALT) GO TO 132 95 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF(IPARM(11).NE.0)GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 28H IN ITPACK ROUTINE SSORSI B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IER = 53 IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF(IPARM(11) .NE. 0)GO TO 134 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 134 IF (LEVEL .GE. 1) WRITE(NOUT,136) IN 136 FORMAT(/1X,25HSSORSI HAS CONVERGED IN ,I5,11H ITERATIONS ) C C ... PUT SOLUTION INTO U IF NOT ALREADY THERE. C 140 CONTINUE IF (MOD(IN,2) .EQ. 1) CALL SCOPY(N,WKSP(IB1),1,U,1) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .LT. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(IB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,142) IERPER 142 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 38H CALLED FROM ITPACK ROUTINE SSORSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IF (IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(5) = OMEGA RPARM(6) = SPECR RPARM(7) = BETAB RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE RSCG(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, RSCG0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE RSCG (REDUCED SYSTEM CONJUGATE C GRADIENT) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C C ... FUNCTION: C C THIS SUBROUTINE, RSCG, DRIVES THE REDUCED SYSTEM CG C ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IN THE RED-BLACK MATRIX. C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. RSCG NEEDS C THIS TO BE IN LENGTH AT LEAST C N+3*NB+2*ITMAX, IF IPARM(5)=0 (SYMMETRIC STORAGE) C N+3*NB+4*ITMAX, IF IPARM(5)=1 (NONSYMMETRIC STORAGE) C HERE NB IS THE ORDER OF THE BLACK SUBSYSTEM C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. IF C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... RSCG SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, CHGCON, DETERM, DFAULT, ECHALL, C ECHOUT, EIGVNS, EIGVSS, EQRT1S, ITERM, ITICK C ITOCK, ITRSCG, IVFILL, PARCON, PERMAT, C PERROR, PERVEC, PMULT, PRBNDX, PRSBLK, C PRSRED, PSTOP, QSORT, SBELM, SCAL, SCOPY, C SDOT, SUM3, UNSCAL, VFILL, VOUT, WEVMW, C ZBRENT C SYSTEM ABS, ALOG10, AMAX0, AMAX1, MOD, SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, IB4, IB5, IDGTS, IER, IERPER, ITIM1, ITIM2, A ITMAX1, JB3, JTIM1, JTIM2, LOOP, N, NB, NDUMMY, A NR, NRP1, N3 REAL DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,41HBEGINNING OF ITPACK SOLUTION MODULE RSCG ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,6) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 26H IN ITPACK ROUTINE RSCG B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF(N .GT. 0 ) GO TO 15 IER = 61 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 25 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 25 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) GO TO 160 C C ... INITIALIZE WKSP BASE ADDRESSES. C 25 IB1 = 1 IB2 = IB1 + N JB3 = IB2 + N C C ... PERMUTE TO RED-BLACK SYSTEM IF POSSIBLE C NB = IPARM(9) IF(NB .GE. 0) GO TO 35 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 35 IF (LEVEL .GE. 0) WRITE(NOUT,30) IER,NB 30 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 35 IF (NB .GE. 0 .AND. NB .LE. N) GO TO 43 IER = 64 IF (LEVEL .GE. 1) WRITE(NOUT,37) IER,NB 37 FORMAT(/10X,43HERROR DETECTED IN RED-BLACK SUBSYSTEM INDEX A /10X,5HIER =,I5,11H IPARM(9) =,I5,5H (NB) ) GO TO 160 43 IF (NB.NE.0 .AND. NB.NE.N) GO TO 45 NB = N/2 IF (LEVEL.GE.2 .AND. IPARM(9).GE.0) WRITE(NOUT,44) A IPARM(9),NB 44 FORMAT(/10X,12H IPARM(9) = ,I5,27H IMPLIES MATRIX IS DIAGONAL A /10X,13H NB RESET TO ,I5) C C ... PERMUTE MATRIX AND RHS C 45 IF(IPARM(9) .GE. 0) GO TO 51 IF (LEVEL .GE. 2) WRITE(NOUT,46) NB 46 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(JB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 48 IF (LEVEL .GE. 0) WRITE (NOUT,47) IER 47 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 48 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... FINISH WKSP BASE ADDRESSES C 51 IB3 = IB2 + NB IB4 = IB3 + NB IB5 = IB4 + NB NR = N - NB NRP1 = NR + 1 IPARM(8) = N + 3*NB + 2*ITMAX IF(ISYM .NE. 0) IPARM(8) = IPARM(8) + 2*ITMAX IF (NW .GE. IPARM(8)) GO TO 55 IER = 62 IF (LEVEL .GE. 0) WRITE(NOUT,53) NW,IPARM(8) 53 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) GO TO 160 C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 55 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 61 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 61 IF (LEVEL .LE. 2) GO TO 65 WRITE(NOUT,62) 62 FORMAT(///1X,35HIN THE FOLLOWING, RHO AND GAMMA ARE A ,24H ACCELERATION PARAMETERS ) IF (ADAPT) WRITE(NOUT,63) 63 FORMAT(1X,48HCME IS THE ESTIMATE OF THE LARGEST EIGENVALUE OF A ,18H THE JACOBI MATRIX ) 65 IF (IPARM(11).NE.0)GO TO 75 ITIM1 = ITICK(NDUMMY) C C ... INITIALIZE FORWARD PSEUDO-RESIDUAL C 75 CONTINUE IF (N .GT. 1) GO TO 76 U(1) = RHS(1) GO TO 132 76 CALL SCOPY(NR,RHS,1,WKSP(IB1),1) CALL PRSRED(NB,NR,IA,JA,A,U(NRP1),WKSP(IB1)) CALL SCOPY(NB,RHS(NRP1),1,WKSP(IB2),1) CALL PRSBLK(NB,NR,IA,JA,A,WKSP(IB1),WKSP(IB2)) CALL VEVMW(NB,WKSP(IB2),U(NRP1)) C C ... ITERATION SEQUENCE C ITMAX1 = ITMAX + 1 DO 95 LOOP=1,ITMAX1 IN = LOOP - 1 IF (MOD(IN,2) .EQ. 1) GO TO 85 C C ... CODE FOR THE EVEN ITERATIONS. C C U = U(IN) WKSP(IB2) = D(IN) C WKSP(IB1) = U(IN-1) WKSP(IB3) = D(IN-1) C CALL ITRSCG(N,NB,IA,JA,A,U,WKSP(IB1),WKSP(IB2),WKSP(IB3), A WKSP(IB4),WKSP(IB5) ) C IF (HALT) GO TO 132 GO TO 95 C C ... CODE FOR THE ODD ITERATIONS. C C U = U(IN-1) WKSP(IB2) = D(IN-1) C WKSP(IB1) = U(IN) WKSP(IB3) = D(IN) C 85 CALL ITRSCG(N,NB,IA,JA,A,WKSP(IB1),U,WKSP(IB3),WKSP(IB2), A WKSP(IB4),WKSP(IB5) ) C IF (HALT) GO TO 132 95 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF (IPARM(11) .NE. 0) GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 26H IN ITPACK ROUTINE RSCG B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IER = 63 IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF (IPARM(11) .NE. 0) GO TO 134 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 134 IF (LEVEL .GE. 1) WRITE(NOUT,136) IN 136 FORMAT(/1X,23HRSCG HAS CONVERGED IN ,I5,11H ITERATIONS ) C C ... PUT SOLUTION INTO U IF NOT ALREADY THERE. C 140 CONTINUE IF (N .EQ. 1) GO TO 142 IF (MOD(IN,2) .EQ. 1) CALL SCOPY(N,WKSP(IB1),1,U,1) CALL SCOPY(NR,RHS,1,U,1) CALL PRSRED(NB,NR,IA,JA,A,U(NRP1),U) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C 142 CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .GE. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(JB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,143) IERPER 143 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSCG B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IF(IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IPARM(8) = IPARM(8) - 2 * (ITMAX - IN) IF(ISYM .NE. 0) IPARM(8) = IPARM(8) - 2*(ITMAX - IN) IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE RSSI(NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP, RSSI0010 A IPARM,RPARM,IERR) C C ITPACK 2C MAIN SUBROUTINE RSSI (REDUCED SYSTEM SEMI-ITERATIVE) C EACH OF THE MAIN SUBROUTINES: C JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI C CAN BE USED INDEPENDENTLY OF THE OTHERS C C ... FUNCTION: C C THIS SUBROUTINE, RSSI, DRIVES THE REDUCED SYSTEM SI C ALGORITHM. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS C THE LATEST ESTIMATE TO THE SOLUTION. C IWKSP INTEGER VECTOR WORKSPACE OF LENGTH 3*N C NW INPUT INTEGER. LENGTH OF AVAILABLE WKSP. ON OUTPUT, C IPARM(8) IS AMOUNT USED. C WKSP REAL VECTOR USED FOR WORKING SPACE. RSSI C NEEDS THIS TO BE IN LENGTH AT LEAST N + NB C HERE NB IS THE ORDER OF THE BLACK SUBSYSTEM C IPARM INTEGER VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY C SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD. IF C RPARM REAL VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME C REAL PARAMETERS WHICH AFFECT THE METHOD. C IER OUTPUT INTEGER. ERROR FLAG. (= IERR) C C ... RSSI SUBPROGRAM REFERENCES: C C FROM ITPACK BISRCH, CHEBY, CHGSI, DFAULT, ECHALL, C ECHOUT, ITERM, ITICK, ITOCK, ITRSSI, IVFILL, C PARSI, PERMAT, PERROR, PERVEC, PMULT, C PRBNDX, PRSBLK, PRSRED, PSTOP, QSORT, C SAXPY, SBELM, SCAL, SCOPY, SDOT, SUM3, C TSTCHG, UNSCAL, VEVMW, VFILL, VOUT, C WEVMW C SYSTEM ABS, ALOG10, AMAX0, AMAX1, FLOAT, MOD C SQRT C C VERSION: ITPACK 2C (MARCH 1982) C C CODE WRITTEN BY: DAVID KINCAID, ROGER GRIMES, JOHN RESPESS C CENTER FOR NUMERICAL ANALYSIS C UNIVERSITY OF TEXAS C AUSTIN, TX 78712 C (512) 471-1242 C C FOR ADDITIONAL DETAILS ON THE C (A) SUBROUTINE SEE TOMS ARTICLE 1982 C (B) ALGORITHM SEE CNA REPORT 150 C C BASED ON THEORY BY: DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN C C REFERENCE THE BOOK: APPLIED ITERATIVE METHODS C L. HAGEMAN, D. YOUNG C ACADEMIC PRESS, 1981 C C ************************************************** C * IMPORTANT NOTE * C * * C * WHEN INSTALLING ITPACK ROUTINES ON A * C * DIFFERENT COMPUTER, RESET SOME OF THE VALUES * C * IN SUBROUTNE DFAULT. MOST IMPORTANT ARE * C * * C * SRELPR MACHINE RELATIVE PRECISION * C * RPARM(1) STOPPING CRITERION * C * * C * ALSO CHANGE SYSTEM-DEPENDENT ROUTINE * C * CPTIME USED IN ITICK AND ITOCK * C * * C ************************************************** C C C SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), IWKSP(1), IPARM(12), A NN, NW, IERR REAL A(1), RHS(NN), U(NN), WKSP(NW), RPARM(12) C C SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IB1, IB2, IB3, IB4, IB5, IDGTS, IER, IERPER, ITIM1, ITIM2, A ITMAX1, JB3, JTIM1, JTIM2, LOOP, N, NB, NDUMMY, A NR, NRP1, N3 REAL DIGIT1, DIGIT2, TEMP, TIME1, TIME2, TOL C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C ... VARIABLES IN COMMON BLOCK - ITCOM1 C C IN - ITERATION NUMBER C IS - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED C ISYM - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH C ITMAX - MAXIMUM NUMBER OF ITERATIONS ALLOWED C LEVEL - LEVEL OF OUTPUT CONTROL SWITCH C NOUT - OUTPUT UNIT NUMBER C C ... VARIABLES IN COMMON BLOCK - ITCOM2 C C ADAPT - FULLY ADAPTIVE PROCEDURE SWITCH C BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA C CASEII - ADAPTIVE PROCEDURE CASE SWITCH C HALT - STOPPING TEST SWITCH C PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH C C ... VARIABLES IN COMMON BLOCK - ITCOM3 C C BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N C BETAB - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX C CME - ESTIMATE OF LARGEST EIGENVALUE C DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N C DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S C FF - ADAPTIVE PROCEDURE DAMPING FACTOR C GAMMA - ACCELERATION PARAMETER C OMEGA - OVERRELAXATION PARAMETER FOR SOR AND SSOR C QA - PSEUDO-RESIDUAL RATIO C QT - VIRTUAL SPECTRAL RADIUS C RHO - ACCELERATION PARAMETER C RRR - ADAPTIVE PARAMETER C SIGE - PARAMETER SIGMA-SUB-E C SME - ESTIMATE OF SMALLEST EIGENVALUE C SPECR - SPECTRAL RADIUS ESTIMATE FOR SSOR C SRELPR - MACHINE RELATIVE PRECISION C STPTST - STOPPING PARAMETER C UDNM - TWO NORM OF U C ZETA - STOPPING CRITERION C C C ... INITIALIZE COMMON BLOCKS C LEVEL = IPARM(2) NOUT = IPARM(4) IF (LEVEL .GE. 1) WRITE(NOUT,1) 1 FORMAT(1H0 ///1X,41HBEGINNING OF ITPACK SOLUTION MODULE RSSI ) IER = 0 IF(IPARM(1) .LE. 0) RETURN N = NN IF(IPARM(11) .EQ. 0) JTIM1 = ITICK(NDUMMY) IF (LEVEL .GE. 3) GO TO 2 CALL ECHOUT(IPARM,RPARM,7) GO TO 3 2 CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,1) 3 TEMP = 5.0E2*SRELPR IF (ZETA .GE. TEMP) GO TO 5 IF (LEVEL .GE. 1) WRITE (NOUT,4) ZETA,SRELPR,TEMP 4 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 26H IN ITPACK ROUTINE RSSI B / 1H , 14H RPARM(1) =,E10.3,7H (ZETA) C / 1H , 46H A VALUE THIS SMALL MAY HINDER CONVERGENCE D / 1H , 36H SINCE MACHINE PRECISION SRELPR = ,E10.3 E / 1H , 18H ZETA RESET TO ,E10.3) ZETA = TEMP 5 CONTINUE TIME1 = RPARM(9) TIME2 = RPARM(10) DIGIT1 = RPARM(11) DIGIT2 = RPARM(12) C C ... VERIFY N C IF (N.GT.0) GO TO 15 IER = 71 IF (LEVEL .GE. 0) WRITE(NOUT,14) N 14 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 33H INVALID MATRIX DIMENSION, N =,I8) GO TO 160 15 CONTINUE C C ... REMOVE ROWS AND COLUMNS IF REQUESTED C IF(IPARM(10) .EQ. 0) GO TO 25 TOL = RPARM(8) CALL IVFILL(N,IWKSP,0) CALL VFILL(N,WKSP,0.0E0) CALL SBELM(N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER) IF(IER .EQ. 0) GO TO 25 IF (LEVEL .GE. 0) WRITE(NOUT,23) IER,TOL 23 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE SBELM C / 1H , 45H WHICH REMOVES ROWS AND COLUMNS OF SYSTEM D / 1H , 35H WHEN DIAGONAL ENTRY TOO LARGE E / 1H , 10H IER = ,I5,5X,12H RPARM(8) = ,E10.3,6H (TOL) ) C C ... INITIALIZE WKSP BASE ADDRESSES. 25 IB1 = 1 IB2 = IB1 + N JB3 = IB2 + N C C ... PERMUTE TO RED-BLACK SYSTEM IF POSSIBLE C NB = IPARM(9) IF(NB .GE. 0) GO TO 35 N3 = 3*N CALL IVFILL(N3,IWKSP,0) CALL PRBNDX(N,NB,IA,JA,IWKSP,IWKSP(IB2),LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 35 IF (LEVEL .GE. 0) WRITE(NOUT,30) IER,NB 30 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PRBNDX C / 1H , 41H WHICH COMPUTES THE RED-BLACK INDEXING D / 1H , 10H IER = ,I5,12H IPARM(9) = ,I5,5H (NB) ) GO TO 160 35 IF (NB .GE. 0 .AND. NB .LE. N) GO TO 43 IER = 74 IF (LEVEL .GE. 1) WRITE(NOUT,37) IER,NB 37 FORMAT(/10X,43HERROR DETECTED IN RED-BLACK SUBSYSTEM INDEX A /10X,5HIER =,I5,11H IPARM(9) =,I5,5H (NB) ) GO TO 160 43 IF (NB.NE.0 .AND. NB.NE.N) GO TO 45 NB = N/2 IF (LEVEL.GE.2 .AND. IPARM(9).GE.0) WRITE(NOUT,44) A IPARM(9),NB 44 FORMAT(/10X,12H IPARM(9) = ,I5,27H IMPLIES MATRIX IS DIAGONAL A /10X,13H NB RESET TO ,I5) C C ... PERMUTE MATRIX AND RHS C 45 IF(IPARM(9) .GE. 0) GO TO 52 IF (LEVEL .GE. 2) WRITE(NOUT,46) NB 46 FORMAT (/10X,27HORDER OF BLACK SUBSYSTEM = ,I5,5H (NB) ) CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(JB3),ISYM,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 51 IF (LEVEL .GE. 0) WRITE (NOUT,48) IER 48 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 40H WHICH DOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) GO TO 160 51 CALL PERVEC(N,RHS,IWKSP) CALL PERVEC(N,U,IWKSP) C C ... INITIALIZE WKSP BASE ADDRESSES C 52 NR = N - NB C NRP1 = NR + 1 IPARM(8) = N + NB IF (NW .GE. IPARM(8)) GO TO 55 IER = 72 IF (LEVEL .GE. 0) WRITE(NOUT,53) NW,IPARM(8) 53 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 28H NOT ENOUGH WORKSPACE AT ,I10 C / 1H , 18H SET IPARM(8) =,I10,5H (NW) ) GO TO 160 C C ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE C ... DIAGONAL ELEMENTS. C 55 CONTINUE CALL VFILL(IPARM(8),WKSP,0.0E0) CALL SCAL(N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER) IF (IER .EQ. 0) GO TO 61 IF (LEVEL .GE. 0) WRITE(NOUT,60) IER 60 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE SCAL C / 1H , 30H WHICH SCALES THE SYSTEM D / 1H , 10H IER = ,I5) GO TO 160 61 IF (LEVEL .LE. 2) GO TO 65 WRITE(NOUT,62) 62 FORMAT(///1X,35HIN THE FOLLOWING, RHO AND GAMMA ARE A ,24H ACCELERATION PARAMETERS ) 65 IF (IPARM(11).NE.0)GO TO 76 ITIM1 = ITICK(NDUMMY) C C ... ITERATION SEQUENCE C 76 IF (N .GT. 1) GO TO 77 U(1) = RHS(1) GO TO 132 77 ITMAX1 = ITMAX + 1 DO 90 LOOP=1,ITMAX1 IN = LOOP - 1 IF (MOD(IN,2) .EQ. 1) GO TO 80 C C ... CODE FOR THE EVEN ITERATIONS. C C U = U(IN) C WKSP(IB1) = U(IN-1) C CALL ITRSSI(N,NB,IA,JA,A,RHS,U,WKSP(IB1),WKSP(IB2) ) C IF (HALT) GO TO 132 GO TO 90 C C ... CODE FOR THE ODD ITERATIONS. C C U = U(IN-1) C WKSP(IB1) = U(IN) C 80 CALL ITRSSI(N,NB,IA,JA,A,RHS,WKSP(IB1),U,WKSP(IB2) ) C IF (HALT) GO TO 132 90 CONTINUE C C ... ITMAX HAS BEEN REACHED C IF (IPARM(11) .NE. 0) GO TO 125 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 125 IF (LEVEL .GE. 1) WRITE(NOUT,130) ITMAX 130 FORMAT(1H0, 30H*** W A R N I N G ************ A / 1H0, 26H IN ITPACK ROUTINE RSSI B / 1H , 26H FAILURE TO CONVERGE IN ,I5, 11H ITERATIONS ) IER = 73 IF (IPARM(3) .EQ. 0) RPARM(1) = STPTST GO TO 140 C C ... METHOD HAS CONVERGED C 132 IF (IPARM(11) .NE. 0) GO TO 134 ITIM2 = ITOCK(NDUMMY) TIME1 = AMAX0(0,ITIM2-ITIM1)*1.E-3 134 IF (LEVEL .GE. 1) WRITE(NOUT,136) IN 136 FORMAT(/1X,23HRSSI HAS CONVERGED IN ,I5,11H ITERATIONS ) C C ... PUT SOLUTION INTO U IF NOT ALREADY THERE. C 140 CONTINUE IF (N .EQ. 1) GO TO 142 IF (MOD(IN,2) .EQ. 1) CALL SCOPY(N,WKSP(IB1),1,U,1) CALL SCOPY(NR,RHS,1,U,1) CALL PRSRED(NB,NR,IA,JA,A,U(NRP1),U) C C ... UNSCALE THE MATRIX, SOLUTION, AND RHS VECTORS. C 142 CALL UNSCAL(N,IA,JA,A,RHS,U,WKSP) C C ... UN-PERMUTE MATRIX,RHS, AND SOLUTION C IF(IPARM(9) .GE. 0) GO TO 145 CALL PERMAT (N,IA,JA,A,IWKSP(IB2),IWKSP(JB3),ISYM, A LEVEL,NOUT,IERPER) IF (IERPER .EQ. 0) GO TO 144 IF (LEVEL .GE. 0) WRITE(NOUT,143) IERPER 143 FORMAT(1H0, 40H*** F A T A L E R R O R ************ A / 1H0, 36H CALLED FROM ITPACK ROUTINE RSSI B / 1H , 40H ERROR DETECTED IN SUBROUTINE PERMAT C / 1H , 45H WHICH UNDOES THE RED-BLACK PERMUTATION D / 1H , 10H IER = ,I5) IF(IER .EQ. 0) IER = IERPER GO TO 160 144 CALL PERVEC(N,RHS,IWKSP(IB2)) CALL PERVEC(N,U,IWKSP(IB2)) C C ... OPTIONAL ERROR ANALYSIS C 145 IDGTS = IPARM(12) IF(IDGTS .LT. 0) GO TO 150 IF(IPARM(2) .LE. 0) IDGTS = 0 CALL PERROR(N,IA,JA,A,RHS,U,WKSP,DIGIT1,DIGIT2,IDGTS) C C ... SET RETURN PARAMETERS IN IPARM AND RPARM C 150 IF(IPARM(11) .NE. 0) GO TO 155 JTIM2 = ITOCK(NDUMMY) TIME2 = AMAX0(0,JTIM2 - JTIM1)*1.E-3 155 IF (IPARM(3) .NE. 0) GO TO 160 IPARM(1) = IN IPARM(9) = NB RPARM(2) = CME RPARM(3) = SME RPARM(9) = TIME1 RPARM(10) = TIME2 RPARM(11) = DIGIT1 RPARM(12) = DIGIT2 C 160 CONTINUE IERR = IER IF (LEVEL .GE. 3) CALL ECHALL(N,IA,JA,A,RHS,IPARM,RPARM,2) C RETURN END SUBROUTINE ITJCG(NN,IA,JA,A,U,U1,D,D1,DTWD,TRI) ITJC0010 C C ... FUNCTION: C C THIS SUBROUTINE, ITJCG, PERFORMS ONE ITERATION OF THE C JACOBI CONJUGATE GRADIENT ALGORITHM. IT IS CALLED BY JCG. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. CONTAINS INFORMATION DEFINING C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. CONTAINS THE NONZERO VALUES OF THE C LINEAR SYSTEM. C U INPUT REAL VECTOR. CONTAINS THE VALUE OF THE C SOLUTION VECTOR AT THE END OF IN ITERATIONS. C U1 INPUT/OUTPUT REAL VECTOR. ON INPUT, IT CONTAINS C THE VALUE OF THE SOLUTION AT THE END OF THE IN-1 C ITERATION. ON OUTPUT, IT WILL CONTAIN THE NEWEST C ESTIMATE FOR THE SOLUTION VECTOR. C D INPUT REAL VECTOR. CONTAINS THE PSEUDO-RESIDUAL C VECTOR AFTER IN ITERATIONS. C D1 INPUT/OUTPUT REAL VECTOR. ON INPUT, D1 CONTAINS C THE PSEUDO-RESIDUAL VECTOR AFTER IN-1 ITERATIONS. ON C OUTPUT, IT WILL CONTAIN THE NEWEST PSEUDO-RESIDUAL C VECTOR. C DTWD REAL ARRAY. USED IN THE COMPUTATIONS OF THE C ACCELERATION PARAMETER GAMMA AND THE NEW PSEUDO- C RESIDUAL. C TRI REAL ARRAY. STORES THE TRIDIAGONAL MATRIX ASSOCIATED C WITH THE EIGENVALUES OF THE CONJUGATE GRADIENT C POLYNOMIAL. C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), A NN REAL A(1), U(NN), U1(NN), D(NN), D1(NN), DTWD(NN), TRI(2,1) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER N REAL CON, C1, C2, C3, C4, DNRM, DTNRM, GAMOLD, RHOOLD, RHOTMP LOGICAL Q1 C C ... SPECIFICATIONS FOR FUNCTION SUBPROGRAMS C REAL SDOT C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE JCG C C ... COMPUTE NEW ESTIMATE FOR CME IF ADAPT = .TRUE. C IF (ADAPT) CALL CHGCON(TRI,GAMOLD,RHOOLD,1) C C ... TEST FOR STOPPING C N = NN DELNNM = SDOT(N,D,1,D,1) DNRM = DELNNM CON = CME CALL PSTOP(N,U,DNRM,CON,1,Q1) IF (HALT) GO TO 16 C C ... COMPUTE RHO AND GAMMA - ACCELERATION PARAMETERS C CALL VFILL(N,DTWD,0.E0) CALL PJAC(N,IA,JA,A,D,DTWD) DTNRM = SDOT(N,D,1,DTWD,1) IF(ISYM .EQ. 0) GO TO 12 RHOTMP = SDOT(N,DTWD,1,D1,1) CALL PARCON(DTNRM,C1,C2,C3,C4,GAMOLD,RHOTMP,1) RHOOLD = RHOTMP GO TO 14 12 CALL PARCON(DTNRM,C1,C2,C3,C4,GAMOLD,RHOOLD,1) C C ... COMPUTE U(IN+1) AND D(IN+1) C 14 CALL SUM3(N,C1,D,C2,U,C3,U1) CALL SUM3(N,C1,DTWD,C4,D,C3,D1) C C ... OUTPUT INTERMEDIATE INFORMATION C 16 CALL ITERM(N,A,U,DTWD,1) C C RETURN END SUBROUTINE ITJSI(NN,IA,JA,A,RHS,U,U1,D,ICNT) ITJS0010 C C ... FUNCTION: C C THIS SUBROUTINE, ITJSI, PERFORMS ONE ITERATION OF THE C JACOBI SEMI-ITERATIVE ALGORITHM. IT IS CALLED BY JSI. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT REAL VECTOR. CONTAINS THE ESTIMATE FOR THE C SOLUTION VECTOR AFTER IN ITERATIONS. C U1 INPUT/OUTPUT REAL VECTOR. ON INPUT, U1 CONTAINS THE C SOLUTION VECTOR AFTER IN-1 ITERATIONS. ON OUTPUT, C IT WILL CONTAIN THE NEWEST ESTIMATE FOR THE SOLUTION C VECTOR. C D REAL ARRAY. D IS USED FOR THE COMPUTATION OF THE C PSEUDO-RESIDUAL ARRAY FOR THE CURRENT ITERATION. C ICNT NUMBER OF ITERATIONS SINCE LAST CHANGE OF SME C C ... SPECIFICATIONS OF ARGUMENTS C INTEGER IA(1), JA(1), A NN, ICNT REAL A(1), RHS(NN), U(NN), U1(NN), D(NN) C C ... SPECIFICATIONS OF LOCAL VARIABLES C INTEGER N REAL CON, C1, C2, C3, DNRM, DTNRM, OLDNRM LOGICAL Q1 C C ... SPECIFICATIONS OF FUNCTION SUBPROGRAMS C REAL SDOT LOGICAL TSTCHG, CHGSME C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE JSI C N = NN IF (IN .EQ. 0) ICNT = 0 C C ... COMPUTE PSEUDO-RESIDUALS C CALL SCOPY(N,RHS,1,D,1) CALL PJAC(N,IA,JA,A,U,D) CALL VEVMW(N,D,U) C C ... STOPPING AND ADAPTIVE CHANGE TESTS C OLDNRM = DELNNM DELNNM = SDOT(N,D,1,D,1) DNRM = DELNNM CON = CME CALL PSTOP(N,U,DNRM,CON,1,Q1) IF (HALT) GO TO 40 IF (.NOT.ADAPT) GO TO 30 IF (.NOT.TSTCHG(1)) GO TO 10 C C ... CHANGE ITERATIVE PARAMETERS (CME) C DTNRM = PVTBV(N,IA,JA,A,D) CALL CHGSI(DTNRM,1) IF (.NOT.ADAPT) GO TO 30 GO TO 20 C C ... TEST IF SME NEEDS TO BE CHANGED AND CHANGE IF NECESSARY. C 10 CONTINUE IF (CASEII) GO TO 30 IF (.NOT. CHGSME(OLDNRM,ICNT)) GO TO 30 ICNT = 0 C C ... COMPUTE U(IN+1) AFTER CHANGE OF PARAMETERS C 20 CALL SCOPY(N,U,1,U1,1) CALL SAXPY(N,GAMMA,D,1,U1,1) GO TO 40 C C ... COMPUTE U(IN+1) WITHOUT CHANGE OF PARAMETERS C 30 CALL PARSI(C1,C2,C3,1) CALL SUM3(N,C1,D,C2,U,C3,U1) C C ... OUTPUT INTERMEDIATE INFORMATION C 40 CALL ITERM(N,A,U,D,2) C RETURN END SUBROUTINE ITSOR (NN,IA,JA,A,RHS,U,WK) ITSO0010 C C ... FUNCTION: C C THIS SUBROUTINE, ITSOR, PERFORMS ONE ITERATION OF THE C SUCCESSIVE OVERRELAXATION ALGORITHM. IT IS CALLED BY SOR. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT/OUTPUT REAL VECTOR. ON INPUT, U CONTAINS THE C SOLUTION VECTOR AFTER IN ITERATIONS. ON OUTPUT, C IT WILL CONTAIN THE NEWEST ESTIMATE FOR THE SOLUTION C VECTOR. C WK REAL ARRAY. WORK VECTOR OF LENGTH N. C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), A NN REAL A(1), RHS(NN), U(NN), WK(NN) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER IP, IPHAT, IPSTAR, ISS, N REAL DNRM, H, OMEGAP, SPCRM1 LOGICAL CHANGE, Q1 C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE SOR C C ... SET INITIAL PARAMETERS NOT ALREADY SET C N = NN IF (IN .NE. 0) GO TO 10 CALL PSTOP(N,U,0.E0,0.E0,0,Q1) IF (ADAPT) GO TO 5 CHANGE = .FALSE. IP = 0 IPHAT = 2 ISS = 0 GO TO 20 C 5 CHANGE = .TRUE. IP = 0 OMEGAP = OMEGA OMEGA = 1.E0 ISS = 0 IPHAT = 2 IPSTAR = 4 IF (OMEGAP .LE. 1.E0) CHANGE = .FALSE. C C ... RESET OMEGA, IPHAT, AND IPSTAR (CIRCLE A IN FLOWCHART) C 10 IF (.NOT. CHANGE) GO TO 20 CHANGE = .FALSE. IS = IS + 1 IP = 0 ISS = 0 OMEGA = AMIN1 (OMEGAP,TAU(IS)) IPHAT = MAX0 ( 3 , IFIX( (OMEGA-1.E0)/(2.E0-OMEGA) ) ) IPSTAR = IPSTR (OMEGA) C C ... COMPUTE U (IN + 1) AND NORM OF DEL(S,P) - CIRCLE B IN FLOW CHART C 20 CONTINUE DELSNM = DELNNM SPCRM1 = SPECR CALL SCOPY(N,RHS,1,WK,1) CALL PFSOR1(N,IA,JA,A,U,WK) IF (DELNNM .EQ. 0.E0) GO TO 30 IF (IN .NE. 0) SPECR = DELNNM / DELSNM IF (IP .LT. IPHAT) GO TO 60 C C ... STOPPING TEST, SET H C IF (SPECR .GE. 1.E0) GO TO 60 IF (.NOT.(SPECR .GT. (OMEGA - 1.E0))) GO TO 30 H = SPECR GO TO 40 30 ISS = ISS + 1 H = OMEGA - 1.E0 C C ... PERFORM STOPPING TEST. C 40 CONTINUE DNRM = DELNNM ** 2 CALL PSTOP(N,U,DNRM,H,1,Q1) IF (HALT) GO TO 60 C C ... METHOD HAS NOT CONVERGED YET, TEST FOR CHANGING OMEGA C IF (.NOT. ADAPT) GO TO 60 IF (IP .LT. IPSTAR) GO TO 60 IF (OMEGA .GT. 1.E0) GO TO 50 CME = SQRT(ABS(SPECR)) OMEGAP = 2.E0 / (1.E0 + SQRT(ABS(1.E0 - SPECR))) CHANGE = .TRUE. GO TO 60 50 IF (ISS .NE. 0) GO TO 60 IF (SPECR .LE. (OMEGA - 1.E0)**FF) GO TO 60 IF ((SPECR + 5.E-5) .LE. SPCRM1) GO TO 60 C C ... CHANGE PARAMETERS C CME = (SPECR + OMEGA - 1.E0) / (SQRT(ABS(SPECR))*OMEGA) OMEGAP = 2.E0 / (1.E0 + SQRT(ABS(1.E0 - CME*CME))) CHANGE = .TRUE. C C ... OUTPUT INTERMEDIATE INFORMATION C 60 CALL ITERM(N,A,U,WK,3) IP = IP + 1 C C RETURN END SUBROUTINE ITSRCG(NN,IA,JA,A,RHS,U,U1,C,C1,D,DL,WK,TRI) ITSR0010 C C ... FUNCTION: C C THIS SUBROUTINE, ITSRCG, PERFORMS ONE ITERATION OF THE C SYMMETRIC SOR CONJUGATE GRADIENT ALGORITHM. IT IS CALLED BY C SSORCG. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT REAL VECTOR. CONTAINS THE ESTIMATE OF THE C SOLUTION VECTOR AFTER IN ITERATIONS. C U1 INPUT/OUTPUT REAL VECTOR. ON INPUT, U1 CONTAINS THE C THE ESTIMATE FOR THE SOLUTION AFTER IN-1 ITERATIONS. C ON OUTPUT, U1 CONTAINS THE UPDATED ESTIMATE. C C INPUT REAL VECTOR. CONTAINS THE FORWARD RESIDUAL C AFTER IN ITERATIONS. C C1 INPUT/OUTPUT REAL VECTOR. ON INPUT, C1 CONTAINS C THE FORWARD RESIDUAL AFTER IN-1 ITERATIONS. ON C OUTPUT, C1 CONTAINS THE UPDATED FORWARD RESIDUAL. C D REAL VECTOR. IS USED TO COMPUTE THE BACKWARD PSEUDO- C RESIDUAL VECTOR FOR THE CURRENT ITERATION. C DL REAL VECTOR. IS USED IN THE COMPUTATIONS OF THE C ACCELERATION PARAMETERS. C WK REAL VECTOR. WORKING SPACE OF LENGTH N. C TRI REAL VECTOR. STORES THE TRIDIAGONAL MATRIX ASSOCIATED C WITH THE CONJUGATE GRADIENT ACCELERATION. C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), A NN REAL A(1), RHS(NN), U(NN), U1(NN), C(NN), C1(NN), D(NN), A DL(NN), WK(NN), TRI(2,1) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER N REAL BETNEW, CON, DNRM, GAMOLD, RHOOLD, RHOTMP, T1, A T2, T3, T4 LOGICAL Q1 C C ... SPECIFICATIONS FOR FUNCTION SUBPROGRAMS C REAL SDOT LOGICAL OMGCHG, OMGSTR C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE SSORCG C C ... CALCULATE S-PRIME FOR ADAPTIVE PROCEDURE. C N = NN IF (ADAPT.OR.PARTAD) CALL CHGCON(TRI,GAMOLD,RHOOLD,3) C C ... COMPUTE BACKWARD RESIDUAL C CALL SCOPY(N,RHS,1,WK,1) CALL SCOPY(N,C,1,D,1) CALL VEVPW(N,D,U) CALL PBSOR(N,IA,JA,A,D,WK) CALL VEVMW(N,D,U) C C ... COMPUTE ACCELERATION PARAMETERS AND THEN U(IN+1) (IN U1) C CALL SCOPY(N,D,1,DL,1) CALL VFILL(N,WK,0.E0) CALL PFSOR(N,IA,JA,A,DL,WK) CALL WEVMW(N,D,DL) DELNNM = SDOT(N,C,1,C,1) IF (DELNNM .EQ. 0.E0) GO TO 5 DNRM = SDOT(N,C,1,DL,1) IF(DNRM .EQ. 0.E0) GO TO 5 IF(ISYM .EQ. 0) GO TO 2 RHOTMP = SDOT(N,C,1,C1,1) - SDOT(N,DL,1,C1,1) CALL PARCON(DNRM,T1,T2,T3,T4,GAMOLD,RHOTMP,3) RHOOLD = RHOTMP GO TO 4 2 CALL PARCON(DNRM,T1,T2,T3,T4,GAMOLD,RHOOLD,3) 4 CALL SUM3(N,T1,D,T2,U,T3,U1) C C ... TEST FOR STOPPING C 5 BDELNM = SDOT(N,D,1,D,1) DNRM = BDELNM CON = SPECR CALL PSTOP(N,U,DNRM,CON,1,Q1) IF (HALT) GO TO 50 C C ... IF NON- OR PARTIALLY-ADAPTIVE, COMPUTE C(IN+1) AND EXIT. C IF (ADAPT) GO TO 10 CALL SUM3(N,-T1,DL,T2,C,T3,C1) GO TO 50 C C ... FULLY ADAPTIVE PROCEDURE C 10 CONTINUE IF (OMGSTR(1)) GO TO 40 IF (OMGCHG(1)) GO TO 20 C C ... PARAMETERS HAVE BEEN UNCHANGED. COMPUTE C(IN+1) AND EXIT. C CALL SUM3(N,-T1,DL,T2,C,T3,C1) GO TO 50 C C ... IT HAS BEEN DECIDED TO CHANGE PARAMETERS C (1) COMPUTE NEW BETAB IF BETADT = .TRUE. C 20 CONTINUE IF (.NOT.BETADT) GO TO 30 BETNEW = PBETA(N,IA,JA,A,D,WK,C1) / BDELNM BETAB = AMAX1(BETAB,.25E0,BETNEW) C C ... (2) COMPUTE NEW CME, OMEGA, AND SPECR C 30 CONTINUE IF(CASEII) GO TO 35 DNRM = PVTBV(N,IA,JA,A,D) GO TO 37 35 CALL VFILL(N,WK,0.E0) CALL PJAC(N,IA,JA,A,D,WK) DNRM = SDOT(N,WK,1,WK,1) 37 CALL OMEG(DNRM,3) C C ... (3) COMPUTE NEW FORWARD RESIDUAL SINCE OMEGA HAS BEEN CHANGED. C 40 CONTINUE CALL SCOPY(N,RHS,1,WK,1) CALL SCOPY(N,U1,1,C1,1) CALL PFSOR(N,IA,JA,A,C1,WK) CALL VEVMW(N,C1,U1) C C ... OUTPUT INTERMEDIATE RESULTS. C 50 CALL ITERM(N,A,U,WK,4) C RETURN END SUBROUTINE ITSRSI(NN,IA,JA,A,RHS,U,U1,C,D,CTWD,WK) ITSR0010 C ... FUNCTION: C C THIS SUBROUTINE, ITSRSI, PERFORMS ONE ITERATION OF THE C SYMMETRIC SOR SEMI-ITERATION ALGORITHM. IT IS CALLED BY C SSORSI. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. (= NN) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C U INPUT REAL VECTOR. CONTAINS THE ESTIMATE OF THE C SOLUTION VECTOR AFTER IN ITERATIONS. C U1 INPUT/OUTPUT REAL VECTOR. ON INPUT, U1 CONTAINS THE C THE ESTIMATE FOR THE SOLUTION AFTER IN-1 ITERATIONS. C ON OUTPUT, U1 CONTAINS THE UPDATED ESTIMATE. C C REAL VECTOR. IS USED TO COMPUTE THE FORWARD PSEUDO- C RESIDUAL VECTOR FOR THE CURRENT ITERATION. C D REAL VECTOR. IS USED TO COMPUTE THE BACKWARD PSEUDO- C RESIDUAL VECTOR FOR THE CURRENT ITERATION. C CTWD REAL VECTOR. IS USED IN THE COMPUTATIONS OF THE C ACCELERATION PARAMETERS. C WK REAL VECTOR. WORKING SPACE OF LENGTH N. C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), A NN REAL A(1), RHS(NN), U(NN), U1(NN), C(NN), D(NN), CTWD(NN), A WK(NN) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER N REAL BETNEW, CON, C1, C2, C3, DNRM LOGICAL Q1 C C ... SPECIFICATIONS FOR FUNCTION SUBPROGRAMS C REAL SDOT LOGICAL OMGSTR, TSTCHG C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE SSORSI C C ... COMPUTE PSEUDO-RESIDUALS (FORWARD AND BACKWARD) C N = NN CALL SCOPY(N,RHS,1,WK,1) CALL SCOPY(N,U,1,CTWD,1) CALL PSSOR1(N,IA,JA,A,CTWD,WK,C,D) C C ... COMPUTE U(IN+1) -- CONTAINED IN THE VECTOR U1. C CALL PARSI(C1,C2,C3,3) CALL SUM3(N,C1,D,C2,U,C3,U1) C C ... TEST FOR STOPPING C BDELNM = SDOT(N,D,1,D,1) DNRM = BDELNM CON = SPECR CALL PSTOP(N,U,DNRM,CON,1,Q1) IF (HALT .OR. .NOT.(ADAPT.OR.PARTAD)) GO TO 20 C C ... ADAPTIVE PROCEDURE C IF (OMGSTR(1)) GO TO 20 DELNNM = SDOT(N,C,1,C,1) IF (IN .EQ. IS) DELSNM = DELNNM IF (IN.EQ.0 .OR. .NOT.TSTCHG(1)) GO TO 20 C C ... IT HAS BEEN DECIDED TO CHANGE PARAMETERS. C ... (1) COMPUTE CTWD C CALL SCOPY(N,D,1,CTWD,1) CALL VFILL(N,WK,0.E0) CALL PFSOR(N,IA,JA,A,CTWD,WK) CALL VEVPW(N,CTWD,C) CALL VEVMW(N,CTWD,D) C C ... (2) COMPUTE NEW SPECTRAL RADIUS FOR CURRENT OMEGA. C DNRM = SDOT(N,C,1,CTWD,1) CALL CHGSI(DNRM,3) IF (.NOT. ADAPT) GO TO 20 C C ... (3) COMPUTE NEW BETAB IF BETADT = .TRUE. C IF (.NOT.BETADT) GO TO 10 BETNEW = PBETA(N,IA,JA,A,D,WK,CTWD) / BDELNM BETAB = AMAX1(BETAB, .25E0, BETNEW) C C ... (4) COMPUTE NEW CME, OMEGA, AND SPECR. C 10 CONTINUE IF(CASEII) GO TO 15 DNRM = PVTBV(N,IA,JA,A,D) GO TO 17 15 CALL VFILL(N,WK,0.E0) CALL PJAC(N,IA,JA,A,D,WK) DNRM = SDOT(N,WK,1,WK,1) 17 CALL OMEG(DNRM,3) C C ... OUTPUT INTERMEDIATE INFORMATION C 20 CALL ITERM(N,A,U,WK,5) C RETURN END SUBROUTINE ITRSCG(N,NNB,IA,JA,A,UB,UB1,DB,DB1,WB,TRI) ITRS0010 C C ... FUNCTION: C C THIS SUBROUTINE, ITRSCG, PERFORMS ONE ITERATION OF THE C REDUCED SYSTEM CONJUGATE GRADIENT ALGORITHM. IT IS C CALLED BY RSCG. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. C NB INPUT INTEGER. CONTAINS THE NUMBER OF BLACK POINTS C IN THE RED-BLACK MATRIX. (= NNB) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C UB INPUT REAL VECTOR. CONTAINS THE ESTIMATE FOR THE C SOLUTION ON THE BLACK POINTS AFTER IN ITERATIONS. C UB1 INPUT/OUTPUT REAL VECTOR. ON INPUT, UB1 CONTAINS THE C SOLUTION VECTOR AFTER IN-1 ITERATIONS. ON OUTPUT, C IT WILL CONTAIN THE NEWEST ESTIMATE FOR THE SOLUTION C VECTOR. THIS IS ONLY FOR THE BLACK POINTS. C DB INPUT REAL ARRAY. DB CONTAINS THE VALUE OF THE C CURRENT PSEUDO-RESIDUAL ON THE BLACK POINTS. C DB1 INPUT/OUTPUT REAL ARRAY. DB1 CONTAINS THE PSEUDO- C RESIDUAL ON THE BLACK POINTS FOR THE IN-1 ITERATION C ON INPUT. ON OUTPUT, IT IS FOR THE IN+1 ITERATION. C WB REAL ARRAY. WB IS USED FOR COMPUTATIONS INVOLVING C BLACK VECTORS. C TRI REAL ARRAY. STORES THE TRIDIAGONAL MATRIX ASSOCIATED C WITH CONJUGATE GRADIENT ACCELERATION. C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), A N, NNB REAL A(1), UB(N), UB1(N), DB(NNB), DB1(N), WB(NNB), TRI(2,1) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER NB, NR, NRP1 REAL CON, C1, C2, C3, C4, DNRM, GAMOLD, RHOOLD, RHOTMP LOGICAL Q1 C C ... SPECIFICATIONS FOR FUNCTION SUBPROGRAMS C REAL SDOT C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE RSCG C C ... COMPUTE NEW ESTIMATE FOR CME IF ADAPT = .TRUE. C NB = NNB NR = N - NB NRP1 = NR + 1 IF (ADAPT) CALL CHGCON(TRI,GAMOLD,RHOOLD,2) C C ... TEST FOR STOPPING C DELNNM = SDOT(NB,DB,1,DB,1) DNRM = DELNNM CON = CME CALL PSTOP(NB,UB(NRP1),DNRM,CON,2,Q1) IF (HALT) GO TO 6 C C ... COMPUTE ACCELERATION PARAMETERS C CALL VFILL(NR,UB1,0.E0) CALL PRSRED(NB,NR,IA,JA,A,DB,UB1) CALL VFILL(NB,WB,0.E0) CALL PRSBLK(NB,NR,IA,JA,A,UB1,WB) DNRM = SDOT(NB,DB,1,WB,1) IF(ISYM .EQ. 0) GO TO 2 RHOTMP = SDOT(NB,WB,1,DB1,1) CALL PARCON(DNRM,C1,C2,C3,C4,GAMOLD,RHOTMP,2) RHOOLD = RHOTMP GO TO 4 2 CALL PARCON(DNRM,C1,C2,C3,C4,GAMOLD,RHOOLD,2) C C ... COMPUTE UB(IN+1) AND DB(IN+1) C 4 CALL SUM3(NB,C1,DB,C2,UB(NRP1),C3,UB1(NRP1)) CALL SUM3(NB,C1,WB,C4,DB,C3,DB1) C C ... OUTPUT INTERMEDIATE INFORMATION C 6 CALL ITERM(NB,A(NRP1), UB(NRP1),WB,6) C RETURN END SUBROUTINE ITRSSI(N,NNB,IA,JA,A,RHS,UB,UB1,DB) ITRS0010 C C ... FUNCTION: C C THIS SUBROUTINE, ITRSSI, PERFORMS ONE ITERATION OF THE C REDUCED SYSTEM SEMI-ITERATION ALGORITHM. IT IS C CALLED BY RSSI. C C ... PARAMETER LIST: C C N INPUT INTEGER. DIMENSION OF THE MATRIX. C NB INPUT INTEGER. CONTAINS THE NUMBER OF BLACK POINTS C IN THE RED-BLACK MATRIX. (= NNB) C IA,JA INPUT INTEGER VECTORS. THE TWO INTEGER ARRAYS OF C THE SPARSE MATRIX REPRESENTATION. C A INPUT REAL VECTOR. THE REAL ARRAY OF THE SPARSE C MATRIX REPRESENTATION. C RHS INPUT REAL VECTOR. CONTAINS THE RIGHT HAND SIDE C OF THE MATRIX PROBLEM. C UB INPUT REAL VECTOR. CONTAINS THE ESTIMATE FOR THE C SOLUTION ON THE BLACK POINTS AFTER IN ITERATIONS. C UB1 INPUT/OUTPUT REAL VECTOR. ON INPUT, UB1 CONTAINS THE C SOLUTION VECTOR AFTER IN-1 ITERATIONS. ON OUTPUT, C IT WILL CONTAIN THE NEWEST ESTIMATE FOR THE SOLUTION C VECTOR. THIS IS ONLY FOR THE BLACK POINTS. C DB INPUT REAL ARRAY. DB CONTAINS THE VALUE OF THE C CURRENT PSEUDO-RESIDUAL ON THE BLACK POINTS. C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IA(1), JA(1), A N, NNB REAL A(1), RHS(N), UB(N), UB1(N), DB(NNB) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER NB, NR, NRP1 REAL CONST, C1, C2, C3, DNRM LOGICAL Q1 C C ... SPECIFICATIONS FOR FUNCTION SUBPROGRAMS C REAL SDOT LOGICAL TSTCHG C C *** BEGIN: ITPACK COMMON C INTEGER IN, IS, ISYM, ITMAX, LEVEL, NOUT COMMON /ITCOM1/ IN, IS, ISYM, ITMAX, LEVEL, NOUT C LOGICAL ADAPT, BETADT, CASEII, HALT, PARTAD COMMON /ITCOM2/ ADAPT, BETADT, CASEII, HALT, PARTAD C REAL BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA, A QA,QT,RHO,RRR,SIGE,SME,SPECR,SPR,SRELPR,STPTST,UDNM,ZETA C C *** END : ITPACK COMMON C C DESCRIPTION OF VARIABLES IN COMMON BLOCKS IN SUBROUTINE RSSI C C ... COMPUTE UR(IN) INTO UB C NB = NNB NR = N - NB NRP1 = NR + 1 CALL SCOPY(NR,RHS,1,UB,1) CALL PRSRED(NB,NR,IA,JA,A,UB(NRP1),UB) C C ... COMPUTE PSEUDO-RESIDUAL, DB(IN) C CALL SCOPY(NB,RHS(NRP1),1,DB,1) CALL PRSBLK(NB,NR,IA,JA,A,UB,DB) CALL VEVMW(NB,DB,UB(NRP1)) C C ... TEST FOR STOPPING C DELNNM = SDOT(NB,DB,1,DB,1) DNRM = DELNNM CONST = CME CALL PSTOP(NB,UB(NRP1),DNRM,CONST,2,Q1) IF (HALT) GO TO 20 IF (.NOT.ADAPT) GO TO 10 C C ... TEST TO CHANGE PARAMETERS C IF (.NOT. TSTCHG(2)) GO TO 10 C C ... CHANGE PARAMETERS C CALL VFILL(NR,UB1,0.E0) CALL PRSRED(NB,NR,IA,JA,A,DB,UB1) DNRM = SDOT(NR,UB1,1,UB1,1) CALL CHGSI(DNRM,2) IF (.NOT.ADAPT) GO TO 10 C C ... COMPUTE UB(N+1) AFTER CHANGING PARAMETERS C CALL SCOPY(NB,UB(NRP1),1,UB1(NRP1),1) CALL SAXPY(NB,GAMMA,DB,1,UB1(NRP1),1) GO TO 20 C C ... COMPUTE UB(N+1) WITHOUT CHANGE OF PARAMETERS C 10 CALL PARSI(C1,C2,C3,2) CALL SUM3(NB,C1,DB,C2,UB(NRP1),C3,UB1(NRP1)) C C ... OUTPUT INTERMEDIATE INFORMATION C 20 CALL ITERM(NB,A(NRP1), UB(NRP1),DB,7) C RETURN END INTEGER FUNCTION BISRCH ( N, K, L ) BISR0010 C C ... BISRCH IS AN INTEGER FUNCTION WHICH USES A BISECTION SEARCH C TO FIND THE ENTRY J IN THE ARRAY K SUCH THAT THE VALUE L IS C GREATER THAN OR EQUAL TO K(J) AND STRICTLY LESS THAN K(J+1). C C ... PARAMETER LIST: C C N INTEGER LENGTH OF VECTOR K C K INTEGER VECTOR C L INTEGER CONSTANT SUCH THAT K(J) .GE. L .LT. K(J+1) C WITH J RETURNED AS VALUE OF INTEGER FUNCTION BISRCH C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER N, L, K(N) C C ... SPECIFICATIONS FOR LOCAL VARIABLES C INTEGER JLEFT, JMID, JRIGHT C JLEFT = 1 JRIGHT = N IF ( N .EQ. 2 ) GO TO 40 JMID = ( N + 1 ) / 2 C 10 IF ( L .GE. K(JMID) ) GO TO 20 C C ...... L .GE. K(LEFT) AND L .LT. K(JMID) C JRIGHT = JMID GO TO 30 C C ...... L .GE. K(JMID) AND L .LT. K(JRIGHT) C 20 JLEFT = JMID C C ...... TEST FOR CONVERGENCE C 30 IF ( JRIGHT - JLEFT .EQ. 1 ) GO TO 40 JMID = JLEFT + ( JRIGHT - JLEFT + 1 ) / 2 GO TO 10 C C ...... BISECTION SEARCH FINISHED C 40 BISRCH = JLEFT C RETURN END REAL FUNCTION CHEBY(QA,QT,RRR,IP,CME,SME) CHEB0010 C C COMPUTES THE SOLUTION TO THE CHEBYSHEV EQUATION C C ... PARAMETER LIST: C C QA RATIO OF PSEUDO-RESIDUALS C QT VIRTUAL SPECTRAL RADIUS C RRR ADAPTIVE PARAMETER C IP NUMBER OF ITERATIONS SINCE LAST CHANGE OF C PARAMETERS C CME, ESTIMATES FOR THE LARGEST AND SMALLEST EIGEN- C SME VALUES OF THE ITERATION MATRIX C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IP REAL CME, QA, QT, RRR, SME C C ... SPECIFICATIONS FOR LOCAL VARIABLES C REAL X, Y, Z C Z = .5E0*(QA + SQRT(ABS(QA**2 - QT**2)))*(1.E0 + RRR**IP) X = Z **(1.E0/FLOAT(IP)) Y = (X + RRR/X)/(1.E0 + RRR) C CHEBY = .5E0*(CME + SME + Y * (2.E0 - CME - SME) ) C RETURN END SUBROUTINE CHGCON(TRI,GAMOLD,RHOOLD,IBMTH) CHGC0010 C C COMPUTES THE NEW ESTIMATE FOR THE LARGEST EIGENVALUE FOR C CONJUGATE GRADIENT ACCELERATION. C C ... PARAMETER LIST: C C TRI TRIDIAGONAL MATRIX ASSOCIATED WITH THE EIGENVALUES C OF THE CONJUGATE GRADIENT POLYNOMIAL C GAMOLD C AND C RHOOLD PREVIOUS VALUES OF ACCELERATION PARAMETERS C IBMTH INDICATOR OF BASIC METHOD BEING ACCELERATED BY CG C IBMTH = 1, JACOBI C = 2, REDUCED SYSTEM C = 3, SSOR C C ... SPECIFICATIONS FOR ARGUMENTS C INTEGER IBMTH REAL TRI(2,1), A GAMOLD, RHOOLD C C ... SPECI