SUBROUTINE JCG (NN,IA,JA,A,RHS,U,IWKSP,NW,WKSP,IPARM,RPARM,IERR)!       !     ITPACK 2C MAIN SUBROUTINE  JCG  (JACOBI CONJUGATE GRADIENT)   !     EACH OF THE MAIN SUBROUTINES:   !           JCG, JSI, SOR, SSORCG, SSORSI, RSCG, RSSI     !     CAN BE USED INDEPENDENTLY OF THE OTHERS   !       ! ... FUNCTION:   !       !          THIS SUBROUTINE, JCG, DRIVES THE JACOBI CONJUGATE!          GRADIENT ALGORITHM.!       ! ... PARAMETER LIST:       !       !          N      INPUT INTEGER.  DIMENSION OF THE MATRIX. (= NN)   !          IA,JA  INPUT INTEGER VECTORS.  THE TWO INTEGER ARRAYS OF !                 THE SPARSE MATRIX REPRESENTATION.       !          A      INPUT D.P. VECTOR.  THE D.P. ARRAY OF THE SPARSE  !                 MATRIX REPRESENTATION.!          RHS    INPUT D.P. VECTOR.  CONTAINS THE RIGHT HAND SIDE  !                 OF THE MATRIX PROBLEM.!          U      INPUT/OUTPUT D.P. VECTOR.  ON INPUT, U CONTAINS THE !                 INITIAL GUESS TO THE SOLUTION. ON OUTPUT, IT CONTAINS !                 THE LATEST ESTIMATE TO THE SOLUTION.    !          IWKSP  INTEGER VECTOR WORKSPACE OF LENGTH 3*N  !          NW     INPUT INTEGER.  LENGTH OF AVAILABLE WKSP.  ON OUTPUT, !                 IPARM(8) IS AMOUNT USED.      !          WKSP   D.P. VECTOR USED FOR WORKING SPACE.  JACOBI CONJUGATE !                 GRADIENT NEEDS THIS TO BE IN LENGTH AT LEAST      !                 4*N + 2*ITMAX,  IF ISYM = 0  (SYMMETRIC STORAGE)  !                 4*N + 4*ITMAX,  IF ISYM = 1  (NONSYMMETRIC STORAGE) !                 HERE ITMAX = IPARM(1) AND ISYM = IPARM(5) !                 (ITMAX IS THE MAXIMUM ALLOWABLE NUMBER OF ITERATIONS) !          IPARM  INTEGER VECTOR OF LENGTH 12.  ALLOWS USER TO SPECIFY!                 SOME INTEGER PARAMETERS WHICH AFFECT THE METHOD.  !          RPARM  D.P. VECTOR OF LENGTH 12. ALLOWS USER TO SPECIFY SOME !                 D.P. PARAMETERS WHICH AFFECT THE METHOD.!          IER    OUTPUT INTEGER.  ERROR FLAG. (= IERR)   !       ! ... JCG SUBPROGRAM REFERENCES:      !       !          FROM ITPACK    BISRCH, CHGCON, DETERM, DFAULT, ECHALL,   !                         ECHOUT, EIGVNS, EIGVSS, EQRT1S, ITERM, TIMER, !                         ITJCG, IVFILL, PARCON, PERMAT,  !                         PERROR5, PERVEC, PJAC, PMULT, PRBNDX,      !                         PSTOP, QSORT, DAXPY, SBELM, SCAL, DCOPY,  !                         DDOT, SUM3, UNSCAL, VEVMW, VFILL, VOUT,   !                         WEVMW, ZBRENT !          SYSTEM         DABS, DLOG10, DBLE(AMAX0), DMAX1, MOD, DSQRT!       !     VERSION:  ITPACK 2C (MARCH 1982)!       !     CODE WRITTEN BY:  DAVID KINCAID, ROGER GRIMES, JOHN RESPESS   !                       CENTER FOR NUMERICAL ANALYSIS     !                       UNIVERSITY OF TEXAS     !                       AUSTIN, TX  78712       !                       (512) 471-1242!       !     FOR ADDITIONAL DETAILS ON THE   !          (A) SUBROUTINE SEE TOMS ARTICLE 1982 !          (B) ALGORITHM  SEE CNA REPORT 150    !       !     BASED ON THEORY BY:  DAVID YOUNG, DAVID KINCAID, LOU HAGEMAN  !       !     REFERENCE THE BOOK:  APPLIED ITERATIVE METHODS      !                          L. HAGEMAN, D. YOUNG !                          ACADEMIC PRESS, 1981 !       !     **************************************************  !     *               IMPORTANT NOTE                   *  !     *                                                *  !     *      WHEN INSTALLING ITPACK ROUTINES ON A      *  !     *  DIFFERENT COMPUTER, RESET SOME OF THE VALUES  *  !     *  IN  SUBROUTNE DFAULT.   MOST IMPORTANT ARE    *  !     *                                                *  !     *   DRELPR      MACHINE RELATIVE PRECISION       *  !     *   RPARM(1)    STOPPING CRITERION               *  !     *                                                *  !     *   ALSO CHANGE SYSTEM-DEPENDENT ROUTINE         *  !     *   SECOND USED IN TIMER                         *  !     *                                                *  !     **************************************************  !       !     SPECIFICATIONS FOR ARGUMENTS    !             INTEGER IA(1),JA(1),IWKSP(1),IPARM(12),NN,NW,IERR         DOUBLE PRECISION A(1),RHS(NN),U(NN),WKSP(NW),RPARM(12)!       !     SPECIFICATIONS FOR LOCAL VARIABLES!             INTEGER IB1,IB2,IB3,IB4,IB5,IDGTS,IER,IERPER,ITMAX1,LOOP,N,NB,N3      DOUBLE PRECISION DIGIT1,DIGIT2,TEMP,TIME1,TIME2,TOL !       ! **** BEGIN: ITPACK COMMON !             INTEGER IN,IS,ISYM,ITMAX,LEVEL,NOUT             COMMON /ITCOM1/ IN,IS,ISYM,ITMAX,LEVEL,NOUT !             LOGICAL ADAPT,BETADT,CASEII,HALT,PARTAD         COMMON /ITCOM2/ ADAPT,BETADT,CASEII,HALT,PARTAD     !             DOUBLE PRECISION BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA,QA,&     &   QT,RHO,RRR,SIGE,SME,SPECR,SPR,DRELPR,STPTST,UDNM,ZETA            COMMON /ITCOM3/ BDELNM,BETAB,CME,DELNNM,DELSNM,FF,GAMMA,OMEGA,QA, &     &   QT,RHO,RRR,SIGE,SME,SPECR,SPR,DRELPR,STPTST,UDNM,ZETA      !       ! **** END  : ITPACK COMMON !       ! ... VARIABLES IN COMMON BLOCK - ITCOM1!       !     IN     - ITERATION NUMBER       !     IS     - ITERATION NUMBER WHEN PARAMETERS LAST CHANGED!     ISYM   - SYMMETRIC/NONSYMMETRIC STORAGE FORMAT SWITCH !     ITMAX  - MAXIMUM NUMBER OF ITERATIONS ALLOWED       !     LEVEL  - LEVEL OF OUTPUT CONTROL SWITCH   !     NOUT   - OUTPUT UNIT NUMBER     !       ! ... VARIABLES IN COMMON BLOCK - ITCOM2!       !     ADAPT  - FULLY ADAPTIVE PROCEDURE SWITCH  !     BETADT - SWITCH FOR ADAPTIVE DETERMINATION OF BETA  !     CASEII - ADAPTIVE PROCEDURE CASE SWITCH   !     HALT   - STOPPING TEST SWITCH   !     PARTAD - PARTIALLY ADAPTIVE PROCEDURE SWITCH!       ! ... VARIABLES IN COMMON BLOCK - ITCOM3!       !     BDELNM - TWO NORM OF B TIMES DELTA-SUPER-N!     BETAB  - ESTIMATE FOR THE SPECTRAL RADIUS OF LU MATRIX!     CME    - ESTIMATE OF LARGEST EIGENVALUE   !     DELNNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION N      !     DELSNM - INNER PRODUCT OF PSEUDO-RESIDUAL AT ITERATION S      !     FF     - ADAPTIVE PROCEDURE DAMPING FACTOR!     GAMMA  - ACCELERATION PARAMETER !     OMEGA  - OVERRELAXATION PARAMETER FOR SOR AND SSOR  !     QA     - PSEUDO-RESIDUAL RATIO  !     QT     - VIRTUAL SPECTRAL RADIUS!     RHO    - ACCELERATION PARAMETER !     RRR    - ADAPTIVE PARAMETER     !     SIGE   - PARAMETER SIGMA-SUB-E  !     SME    - ESTIMATE OF SMALLEST EIGENVALUE  !     SPECR  - SPECTRAL RADIUS ESTIMATE FOR SSOR!     DRELPR - MACHINE RELATIVE PRECISION       !     STPTST - STOPPING PARAMETER     !     UDNM   - TWO NORM OF U!     ZETA   - STOPPING CRITERION     !       ! ... INITIALIZE COMMON BLOCKS!             LEVEL = IPARM(2)            NOUT = IPARM(4)             IF (LEVEL.GE.1) WRITE (NOUT,10)    10 FORMAT ('0'///1X,'BEGINNING OF ITPACK SOLUTION MODULE  JCG')        IER = 0           IF (IPARM(1).LE.0) RETURN             N = NN            IF (IPARM(11).EQ.0) TIMJ1 = TIMER(DUMMY)        IF (LEVEL.GE.3) GO TO 20      CALL ECHOUT (IPARM,RPARM,1)           GO TO 30       20 CALL ECHALL (N,IA,JA,A,RHS,IPARM,RPARM,1)    30 TEMP = 5.0D2*DRELPR         IF (ZETA.GE.TEMP) GO TO 50            IF (LEVEL.GE.1) WRITE (NOUT,40) ZETA,DRELPR,TEMP       40 FORMAT ('0','*** W A R N I N G ************'/'0',   &     &   '    IN ITPACK ROUTINE JCG'/' ','    RPARM(1) =',D10.3,    &     &   ' (ZETA)'/' ','    A VALUE THIS SMALL MAY HINDER CONVERGENCE '/&     &   ' ','    SINCE MACHINE PRECISION DRELPR =',D10.3/' ',      &     &   '    ZETA RESET TO ',D10.3)        ZETA = TEMP    50 CONTINUE          TIME1 = RPARM(9)            TIME2 = RPARM(10)           DIGIT1 = RPARM(11)          DIGIT2 = RPARM(12)    !       ! ... VERIFY N    !             IF (N.GT.0) GO TO 70        IER = 11          IF (LEVEL.GE.0) WRITE (NOUT,60) N    60 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &     &   '    CALLED FROM ITPACK ROUTINE JCG '/' ',       &     &   '    INVALID MATRIX DIMENSION, N =',I8)      GO TO 370      70 CONTINUE    !       ! ... REMOVE ROWS AND COLUMNS IF REQUESTED      !             IF (IPARM(10).EQ.0) GO TO 90          TOL = RPARM(8)      CALL IVFILL (N,IWKSP,0)       CALL VFILL (N,WKSP,0.0D0)             CALL SBELM (N,IA,JA,A,RHS,IWKSP,WKSP,TOL,ISYM,LEVEL,NOUT,IER)       IF (IER.EQ.0) GO TO 90      IF (LEVEL.GE.0) WRITE (NOUT,80) IER,TOL      80 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &     &   '    CALLED FROM ITPACK ROUTINE JCG '/' ',       &     &   '    ERROR DETECTED IN SUBROUTINE  SBELM '/' ',  &     &   '    WHICH REMOVES ROWS AND COLUMNS OF SYSTEM '/' ',       &     &   '    WHEN DIAGONAL ENTRY TOO LARGE  '/' ','    IER = ',I5,5X,&     &   ' RPARM(8) = ',D10.3,' (TOL)')       GO TO 370   !       ! ... INITIALIZE WKSP BASE ADDRESSES. !          90 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 110         IER = 12          IF (LEVEL.GE.0) WRITE (NOUT,100) NW,IPARM(8)  100 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &     &   '    CALLED FROM ITPACK ROUTINE JCG '/' ',       &     &   '    NOT ENOUGH WORKSPACE AT ',I10/' ','    SET IPARM(8) =',I10&     &   ,' (NW)')      GO TO 370   !       ! ... PERMUTE TO  RED-BLACK SYSTEM IF REQUESTED !         110 NB = IPARM(9)       IF (NB.LT.0) GO TO 170      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 130       IF (LEVEL.GE.0) WRITE (NOUT,120) IER,NB     120 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &     &   '    CALLED FROM ITPACK ROUTINE JCG '/' ',       &     &   '    ERROR DETECTED IN SUBROUTINE  PRBNDX'/' ',  &     &   '    WHICH COMPUTES THE RED-BLACK INDEXING'/' ','    IER = ',I5&     &   ,' IPARM(9) = ',I5,' (NB)')        GO TO 370   !       ! ... PERMUTE MATRIX AND RHS!         130 IF (LEVEL.GE.2) WRITE (NOUT,140) NB         140 FORMAT (/10X,'ORDER OF BLACK SUBSYSTEM = ',I5,' (NB)')      CALL PERMAT (N,IA,JA,A,IWKSP,IWKSP(IB3),ISYM,LEVEL,NOUT,IER)        IF (IER.EQ.0) GO TO 160       IF (LEVEL.GE.0) WRITE (NOUT,150) IER        150 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &     &   '    CALLED FROM ITPACK ROUTINE JCG '/' ',       &     &   '    ERROR DETECTED IN SUBROUTINE  PERMAT'/' ',  &     &   '    WHICH DOES THE RED-BLACK PERMUTATION'/' ','    IER = ',I5)      GO TO 370     160 CALL PERVEC (N,RHS,IWKSP)             CALL PERVEC (N,U,IWKSP) !       ! ... SCALE LINEAR SYSTEM, U, AND RHS BY THE SQUARE ROOT OF THE     ! ... DIAGONAL ELEMENTS.    !         170 CONTINUE          CALL VFILL (IPARM(8),WKSP,0.0D0)      CALL SCAL (N,IA,JA,A,RHS,U,WKSP,LEVEL,NOUT,IER)           IF (IER.EQ.0) GO TO 190       IF (LEVEL.GE.0) WRITE (NOUT,180) IER        180 FORMAT ('0','*** F A T A L     E R R O R ************'/'0',   &     &   '    CALLED FROM ITPACK ROUTINE JCG '/' ',       &     &   '    ERROR DETECTED IN SUBROUTINE  SCAL  '/' ',  &     &   '    WHICH SCALES THE SYSTEM   '/' ','    IER = ',I5)            GO TO 370     190 IF (LEVEL.LE.2) GO TO 220             WRITE (NOUT,200)        200 FORMAT (///1X,'IN THE FOLLOWING, RHO AND GAMMA ARE',&     &   ' ACCELERATION PARAMETERS')        IF (ADAPT) WRITE (NOUT,210)       210 FORMAT (1X,'CME IS THE ESTIMATE OF THE LARGEST EIGENVALUE OF',&     &   ' THE JACOBI MATRIX')  220 IF (IPARM(11).NE.0) GO TO 230         TIMI1 = TIMER(DUMMY)  !       ! ... COMPUTE INITIAL PSEUDO-RESIDUAL !         230 CONTINUE          CALL DCOPY (N,RHS,1,WKSP(IB2),1)      CALL PJAC (N,IA,JA,A,U,WKSP(IB2))       CALL VEVMW (N,WKSP(IB2),U)      !       ! ... ITERATION SEQUENCE    !             ITMAX1 = ITMAX+1            DO 250 LOOP = 1,ITMAX1         IN = LOOP-1         IF (MOD(IN,2).EQ.1) GO TO 240!       ! ... CODE FOR THE EVEN ITERATIONS.   !       !     U           = U(IN)             WKSP(IB2) = DEL(IN) !     WKSP(IB1)   = U(IN-1)           WKSP(IB3) = DEL(IN-1) !                CALL ITJCG (N,IA,JA,A,U,WKSP(IB1),WKSP(IB2),WKSP(IB3),WKSP(IB4)&     &      ,WKSP(IB5))     !                IF (HALT) GO TO 280         GO TO 250!       ! ... CODE FOR THE ODD ITERATIONS.    !       !     U           = U(IN-1)           WKSP(IB2) = DEL(IN-1) !     WKSP(IB1)   = U(IN)             WKSP(IB3) = DEL(IN) !         240    CALL ITJCG (N,IA,JA,A,WKSP(IB1),U,WKSP(IB3),WKSP(IB2),WKSP(IB4)&     &      ,WKSP(IB5))     !                IF (HALT) GO TO 280  250 CONTINUE    !       ! ... ITMAX HAS BEEN REACHED!             IF (IPARM(11).NE.0) GO TO 260         TIMI2 = TIMER(DUMMY)        TIME1 = DBLE(TIMI2-TIMI1)         260 IER = 13          IF (LEVEL.GE.1) WRITE (NOUT,270) ITMAX      270 FORMAT ('0','*** W A R N I N G ************'/'0',   &     &   '    IN ITPACK ROUTINE JCG'/' ','    FAILURE TO CONVERGE IN',I5&     &   ,' ITERATIONS')          IF (IPARM(3).EQ.0) RPARM(1) = STPTST            GO TO 310   !       ! ... METHOD HAS CONVERGED  !         280 IF (IPARM(11).NE.0) GO TO 290         TIMI2 = TIMER(DUMMY)        TIME1 = DBLE(TIMI2-TIMI1)