SUBROUTINE PCGQS(NN,AA,LENAA,MAXA,PP,YP,RHO,START,WORK,IFLAG)
C
C THIS SUBROUTINE SOLVES A SYSTEM OF EQUATION USING THE METHOD OF
C CONJUGATE GRADIENTS.  THE SYSTEM TO BE SOLVED IS IN THE FORM
C
C   (AUG)*X = B,
C
C WHERE
C
C            +--          --+        +-    -+
C            |        |     |        |      |
C            |   AA   | -PP |        |      |
C      AUG = |        |     | ,  B = | -RHO |
C            +--------------+        |      |
C            |      YP      |        |      |
C            +--          --+        +-    -+
C
C
C
C THE SYSTEM IS SOLVED BY SPLITTING  AUG  INTO TWO MATRICES
C AUG = M + L,  WHERE
C
C       +-        -+                                +-     -+
C       |      |   |                                |       |
C       |  AA  | C |                                | -PP-C |
C   M = |      |   |  ,  L = U * [E(NN+1)**T],  U = |       | ,
C       +------+---+                                +-------+
C       |  C   | D |                                |   0   |
C       +-        -+                                +-     -+
C
C E(NN+1) IS THE (NN+1) X 1 VECTOR CONSISTING OF ALL ZEROS EXCEPT FOR
C A '1' IN ITS LAST POSITION.
C
C THE FINAL SOLUTION VECTOR,  X,  IS GIVEN BY
C
C            +-                                    -+
C            |           [SOL(U)]*[E(NN+1)**T]      |
C       X =  | I  -  -----------------------------  | * SOL(B)
C            |        {[(SOL(U))**T]*E(NN+1)}+1.0   |
C            +-                                    -+
C
C WHERE SOL(A)=[M**(-1)]*A.  THE TWO SYSTEMS  (MZ=U, MZ=B) ARE SOLVED
C BY A PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
C
C
C ON INPUT:
C
C NN  = THE DIMENSION OF THE MATRIX PACKED IN AA.
C
C AA(1:LENAA)  CONTAINS THE MATRIX  AA,  STORED IN PACKED SKYLINE
C    FORMAT.  LENAA  AND  MAXA  DESCRIBE THE DATA STRUCTURE.
C
C LENAA  = THE LENGTH OF THE ONE-DIMENSIONAL ARRAY  AA.
C
C MAXA(1:NN+2)  IN ITS FIRST  N+1  COMPONENTS CONTAINS THE INDICES OF
C    THE DIAGONAL ELEMENTS OF THE MATRIX STORED IN  AA.
C    MAXA(NN+2)  IS ASSIGNED THE VALUE  LENNAA + NN + 2.
C
C    AS AN EXAMPLE OF THE PACKED SKYLINE STORAGE FORMAT, CONSIDER THE
C    SYMMETRIC MATRIX
C
C
C                +--             --+
C                |  1  3  0  0  0  |
C                |  3  2  0  7  0  |
C                |  0  0  4  6  0  |
C                |  0  7  6  5  9  |
C                |  0  0  0  9  8  |
C                +--             --+
C
C    THIS WOULD RESULT IN  NN=5, LENAA=9, MAXA=(1,2,4,5,8,10,16),
C    AND  AA=(1,2,3,4,5,6,7,8,9).
C
C PP(1:NN)  = THE NEGATIVE OF THE LAST COLUMN OF  AUG.
C
C YP(1:NN+1)  = THE LAST ROW OF  AUG.
C
C RHO(1:NN+1)  = THE NEGATIVE OF THE RIGHT HAND SIDE OF THE EQUATION TO
C    BE SOLVED.
C
C WORK(1:6*(NN+1)+LENAA)  IS A WORK ARRAY DIVIDED UP AS FOLLOWS:
C
C    WORK(1:NN+1) = TEMPORARY WORKING STORAGE.
C
C    WORK(NN+2:2*NN+2) = INTERMEDIATE SOLUTION VECTOR Z FOR  MZ=U.
C       THE INPUT VALUE IS USED AS THE INITIAL ESTIMATE FOR Z.
C
C    WORK(2*NN+3:3*NN+3) = INTERMEDIATE SOLUTION VECTOR Z FOR  MZ=B.
C
C    WORK(3*NN+4:4*NN+4) = STORAGE FOR THE RESIDUAL VECTORS.
C
C    WORK(4*NN+5:5*NN+5) = STORAGE FOR THE DIRECTION VECTORS.
C
C    WORK(5*NN+6:6*NN+6+LENAA) = STORAGE FOR THE PRECONDITIONING
C       MATRIX Q.
C
C
C ON OUTPUT:
C
C NN, AA, LENAA, MAXA, PP, YP, AND RHO  ARE UNCHANGED.
C
C START(1:N+1)  CONTAINS THE SOLUTION VECTOR  X.
C
C IFLAG  IS UNCHANGED UNLESS THE CONJUGATE GRADIENT ITERATION
C    FAILS TO CONVERGE IN 10*(NN+1) ITERATIONS (MOST LIKELY DUE
C    TO A SINGULAR JACOBIAN MATRIX).  IN THIS CASE,  PCGQS  RETURNS
C    IFLAG = 4, AND DOES NOT COMPUTE  X.
C
C
C CALLS  D1MACH, DAXPY, DCOPY, DDOT, DNRM2, DSCAL, GMFADS, MULTDS,
C    SOLVDS.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION D1MACH, DDOT, DNRM2
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION AB, AU, BB, BU, DZNRM, PBNPRD, PUNPRD,
     $     RBNPRD, RBTOL, RNPRD, RUNPRD, RUTOL, TEMP, ZLEN, ZTOL
        INTEGER DIR, IMAX, J, LENQ, NP1, Q, RES, ZB, ZU
        LOGICAL STILLU, STILLB
C
C     SCALAR ARGUMENTS
C
        INTEGER NN, LENAA, IFLAG
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION AA(LENAA), PP(NN), YP(NN+1), RHO(NN+1),
     $     START(NN+1), WORK(6*(NN+1)+LENAA)
        INTEGER MAXA(NN+2)
C
C ***** END OF DECLARATIONS *****
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
C SET UP BASES FOR VECTORS STORED IN WORK ARRAY.
C
      NP1=NN+1
      ZU=NN+2
      ZB=(2*NN)+3
      RES=(3*NN)+4
      DIR=(4*NN)+5
      Q=(5*NN)+6
C
C INITIALIZE  PRECONDITIONING MATRIX  Q,  SET VALUES OF  MAXA(NN+1)
C AND  MAXA(NN+2),  COMPUTE PRECONDITIONER.
C
      CALL DCOPY(LENAA,AA,1,WORK(Q),1)
      CALL DCOPY(NP1,YP,1,WORK(Q+LENAA),1)
      MAXA(NN+1)=LENAA+1
      MAXA(NN+2)=LENAA+NN+2
      LENQ = MAXA(NN+2)-1
      CALL GMFADS(NP1,WORK(Q),LENQ,MAXA)
C
C COMPUTE ALL TOLERANCES NEEDED FOR EXIT CRITERIA.
C
      CALL DCOPY(NN,PP,1,WORK,1)
      WORK(NP1)=0.0
      CALL DAXPY(NP1,1.0D0,YP,1,WORK,1)
      IMAX=10*NP1
      STILLU=.TRUE.
      STILLB=.TRUE.
      ZTOL=100.0*D1MACH(4)
      RUTOL=ZTOL*DNRM2(NP1,WORK,1)
      RBTOL=ZTOL*DNRM2(NP1,RHO,1)
C
C ***** END OF INITIALIZATION *****
C
C ***** SOLVE SYSTEM  M Z = U *****
C
C COMPUTE INITIAL RESIDUAL VECTOR FOR THE SYSTEM M Z = U .
C     RES = (Q**(-1))*(U - M*Z.)
*
      CALL MULTDS(WORK(RES),AA,WORK(ZU),MAXA,NN,LENAA)
      WORK(RES+NN)= DDOT(NN,YP,1,WORK(ZU),1)
      CALL DAXPY(NP1,WORK(ZU+NN),YP,1,WORK(RES),1)
      CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
      CALL DAXPY(NN,-1.0D0,PP,1,WORK(RES),1)
      CALL DAXPY(NN,-1.0D0,YP,1,WORK(RES),1)
      CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
C
C COMPUTE INITIAL DIRECTION VECTOR.
C     DIR = (A**T)*(Q**(-T))*RES.
C
      CALL DCOPY(NP1,WORK(RES),1,WORK,1)
      CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
      CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
      WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
      CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
C
C COMPUTE INITIAL INNER PRODUCTS.
C
      RUNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
      PUNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
C
C REPEAT UNTIL CONVERGENCE OR TOO MANY ITERATIONS.
C
      J=1
C
C     DO WHILE ((STILLU) .AND. (J .LE. IMAX))
100   IF (.NOT. ((STILLU) .AND. (J .LE. IMAX)) ) GO TO 200
C
C        IF ||RESIDUAL|| IS STILL NOT SMALL ENOUGH, CONTINUE.
C
         IF (SQRT(RUNPRD) .GT. RUTOL) THEN
C
C           IF DIRECTION VECTOR IS ZERO, THEN RE-COMPUTE RESIDUAL,
C           DIRECTION VECTOR, AND INNER PRODUCTS FROM SCRATCH
C           (RATHER THAN FROM UPDATES OF PREVIOUS VALUES).
C
            IF (PUNPRD .EQ. 0.0) THEN
C
C              COMPUTE RESIDUAL.
C
               CALL MULTDS(WORK(RES),AA,WORK(ZU),MAXA,NN,LENAA)
               WORK(RES+NN)= DDOT(NN,YP,1,WORK(ZU),1)
               CALL DAXPY(NP1,WORK(ZU+NN),YP,1,WORK(RES),1)
               CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
               CALL DAXPY(NN,-1.0D0,PP,1,WORK(RES),1)
               CALL DAXPY(NN,-1.0D0,YP,1,WORK(RES),1)
               CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
C
C              COMPUTE DIRECTION VECTOR.
C
               CALL DCOPY(NP1,WORK(RES),1,WORK,1)
               CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
               CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
               WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
               CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
C
C              COMPUTE INNER PRODUCTS
C
               RUNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
               PUNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
C
C              CHECK FOR CONVERGENCE.
C
               IF (SQRT(RUNPRD) .LE. RUTOL) THEN
                  STILLU=.FALSE.
               ENDIF
            ENDIF
            IF (STILLU) THEN
C
C              UPDATE SOLUTION VECTOR.
C                Z = Z + AU*DIR, WHERE AU= RUNPRD/PUNPRD.
C
               AU=RUNPRD/PUNPRD
               CALL DAXPY(NP1,AU,WORK(DIR),1,WORK(ZU),1)
C
C              COMPUTE RELATIVE CHANGE IN THE SOLUTION.
C
               DZNRM=AU*SQRT(PUNPRD)
               ZLEN=DNRM2(NP1,WORK(ZU),1)
C
C              IF RELATIVE CHANGE IN SOLUTIONS IS SMALL ENOUGH, EXIT.
C
               IF ( (DZNRM/ZLEN) .LT. ZTOL) STILLU=.FALSE.
            ENDIF
         ELSE
            STILLU=.FALSE.
         ENDIF
C
C        IF NO EXIT CRITERIA FOR  MZ=U  HAVE BEEN MET, UPDATE RESIDUAL,
C          DIRECTION VECTORS, AND INNER PRODUCTS FOR NEXT ITERATION.
C
         IF (STILLU) THEN
C
C           UPDATE RESIDUAL VECTOR; COMPUTE <RES,RES>.
C             RES = RES - AU*(Q**(-1))*M*DIR.
C
            CALL MULTDS(WORK,AA,WORK(DIR),MAXA,NN,LENAA)
            WORK(NP1)=DDOT(NN,YP,1,WORK(DIR),1)
            CALL DAXPY(NP1,WORK(DIR+NN),YP,1,WORK,1)
            CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
            CALL DAXPY(NP1,-AU,WORK,1,WORK(RES),1)
            RNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
C
C           UPDATE DIRECTION VECTOR; COMPUTE <DIR,DIR>.
C             DIR = (M**T)*(Q**(-T))*RES + BU*DIR,
C             WHERE  BU = RNPRD/RUNPRD. (NOTE: START IS USED AS
C             A WORK ARRAY HERE).
C
            BU=RNPRD/RUNPRD
            RUNPRD=RNPRD
            CALL DCOPY(NP1,WORK(RES),1,WORK,1)
            CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
            CALL MULTDS(START,AA,WORK,MAXA,NN,LENAA)
            START(NP1)=DDOT(NN,YP,1,WORK,1)
            CALL DAXPY(NP1,WORK(NP1),YP,1,START,1)
            CALL DAXPY(NP1,BU,WORK(DIR),1,START,1)
            CALL DCOPY(NP1,START,1,WORK(DIR),1)
            PUNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
         ENDIF
C
         J=J+1
      GO TO 100
200   CONTINUE
C     END DO
C
C SET ERROR FLAG IF THE CONJUGATE GRADIENT ITERATION DID NOT CONVERGE.
C
      IF (J .GT. IMAX) THEN
         IFLAG=4
         RETURN
      ENDIF
C
C ***** END OF M Z = U SYSTEM *****
C
C ***** SOLVE SYSTEM M Z = B *****
C
C
C COMPUTE INITIAL RESIDUAL VECTOR FOR THE SYSTEM M Z = B .
C
      CALL MULTDS(WORK(RES),AA,WORK(ZB),MAXA,NN,LENAA)
      WORK(RES+NN)=DDOT(NN,YP,1,WORK(ZB),1)
      CALL DAXPY(NP1,WORK(ZB+NN),YP,1,WORK(RES),1)
      CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
      CALL DAXPY(NP1,-1.0D0,RHO,1,WORK(RES),1)
      CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
C
C COMPUTE INITIAL DIRECTION VECTOR.
C
      CALL DCOPY(NP1,WORK(RES),1,WORK,1)
      CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
      CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
      WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
      CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
C
C COMPUTE INITIAL INNER PRODUCTS.
C
      RBNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
      PBNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
C
C REPEAT UNTIL CONVERGENCE, OR TOO MANY ITERATIONS.
C
      J=1
C
C     DO WHILE ( STILLB  .AND. (J .LE. IMAX) )
300   IF (.NOT. ( STILLB  .AND. (J .LE. IMAX) ) ) GO TO 400
C
C        IF ||RESIDUAL|| IS STILL NOT SMALL ENOUGH, CONTINUE.
C
         IF (SQRT(RBNPRD) .GT. RBTOL) THEN
C
C           IF DIRECTION VECTOR IS ZERO, RE-COMPUTE RESIDUAL,
C             DIRECTION VECTOR, AND INNER PRODUCTS FROM SCRATCH.
C
            IF (PBNPRD .EQ. 0.0) THEN
C
C              COMPUTE RESIDUAL.
C
               CALL MULTDS(WORK(RES),AA,WORK(ZB),MAXA,NN,LENAA)
               WORK(RES+NN)=DDOT(NN,YP,1,WORK(ZB),1)
               CALL DAXPY(NP1,WORK(ZB+NN),YP,1,WORK(RES),1)
               CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
               CALL DAXPY(NP1,-1.0D0,RHO,1,WORK(RES),1)
               CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
C
C              COMPUTE DIRECTION VECTOR.
C
               CALL DCOPY(NP1,WORK(RES),1,WORK,1)
               CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
               CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
               WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
               CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
C
C              COMPUTE INNER PRODUCTS.
C
               RBNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
               PBNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
C
C              CHECK FOR CONVERGENCE.
C
               IF (SQRT(RBNPRD) .LE. RBTOL) THEN
                  STILLB=.FALSE.
               ENDIF
            ENDIF
            IF (STILLB) THEN
C
C              UPDATE SOLUTION VECTOR.
C                Z = Z + AB*DIR, WHERE AB=RBNPRD/PBNPRD.
C
               AB=RBNPRD/PBNPRD
               CALL DAXPY(NP1,AB,WORK(DIR),1,WORK(ZB),1)
C
C              COMPUTE RELATIVE CHANGE IN SOLUTIONS.
C
               DZNRM=AB*SQRT(PBNPRD)
               ZLEN=DNRM2(NP1,WORK(ZB),1)
C
C              IF RELATIVE CHANGE IN SOLUTIONS IS SMALL ENOUGH, EXIT.
C
               IF ( (DZNRM/ZLEN) .LT. ZTOL) STILLB=.FALSE.
            ENDIF
         ELSE
            STILLB=.FALSE.
         ENDIF
C
C        IF NO EXIT CRITERIA FOR  MZ=B  HAVE BEEN MET, UPDATE RESIDUAL,
C          DIRECTION VECTORS, AND INNER PRODUCTS.
C
         IF (STILLB) THEN
C
C           UPDATE RESIDUAL VECTOR; COMPUTE <RES,RES>.
C             RES = RES - AB*(Q**(-1))*M*DIR.
C
            CALL MULTDS(WORK,AA,WORK(DIR),MAXA,NN,LENAA)
            WORK(NP1)=DDOT(NN,YP,1,WORK(DIR),1)
            CALL DAXPY(NP1,WORK(DIR+NN),YP,1,WORK,1)
            CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
            CALL DAXPY(NP1,-AB,WORK,1,WORK(RES),1)
            RNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
C
C           UPDATE DIRECTION VECTOR; COMPUTE <DIR,DIR>.
C             DIR = (M**T)*(Q**(-T))*RES + BB*DIR,
C             WHERE BB=RNPRD/RBNPRD.
C           (NOTE:  START IS USED AS A WORK ARRAY HERE).
C
            BB=RNPRD/RBNPRD
            RBNPRD=RNPRD
            CALL DCOPY(NP1,WORK(RES),1,WORK,1)
            CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
            CALL MULTDS(START,AA,WORK,MAXA,NN,LENAA)
            START(NP1)=DDOT(NN,YP,1,WORK,1)
            CALL DAXPY(NP1,WORK(NP1),YP,1,START,1)
            CALL DAXPY(NP1,BB,WORK(DIR),1,START,1)
            CALL DCOPY(NP1,START,1,WORK(DIR),1)
            PBNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
         ENDIF
C
         J=J+1
      GO TO 300
400   CONTINUE
C     END DO
C
C SET ERROR FLAG IF THE CONJUGATE GRADIENT ITERATION DID NOT CONVERGE.
C
      IF (J .GT. IMAX) THEN
         IFLAG=4
         RETURN
      ENDIF
C
C ***** END OF  M Z = B  SYSTEM *****
C
C COMPUTE FINAL SOLUTION VECTOR  X,  AND RETURN IT IN START.
C
      TEMP=-WORK(ZB+NN)/(1.0D0+WORK(ZU+NN))
      CALL DCOPY(NP1,WORK(ZB),1,START,1)
      CALL DAXPY(NP1,TEMP,WORK(ZU),1,START,1)
C
      RETURN
      END

.