C     TEST OF ALGORITHM 24
      INTEGER N,IMULT,NIN,NOUT
      REAL B(100),C(100),G(100),T(100),V(100),TOL 
      EXTERNAL FRANK
C  USING FRANK MATRIX AND UNIT SOLN
C  I/O CHANNELS
      NIN=15
      NOUT=16
   1  READ(NIN,900)N,IMULT,TOL
 900  FORMAT(2I5,F10.5)
      WRITE(NOUT,960)N,IMULT,TOL
 960  FORMAT('0ORDER=',I5,' IMULT=',I5,' TOLERANCE=',1PE16.8)
      IF(N.LE.0)STOP
      CALL SETUP(N,B,C,FRANK) 
      CALL A24CG(N,B,C,TOL,G,NOUT,FRANK,V,T,IMULT)
      CALL TESCG(N,B,C,FRANK,V,NOUT)
      GOTO 1
      END 
      SUBROUTINE SETUP(N,B,C,APR)
C  PUTS RHS = APR*B WHERE B = VECTOR OF ONES
C  PUTS B   = NULL VECTOR
C  APR  = NAME OF MATRIX*VECTOR MULTIPLY
C  N=ORDER OF PROBLEM
      INTEGER N,I
      REAL B(N),C(N)
      DO 10 I=1,N
        B(I)=1.0
  10  CONTINUE
      CALL APR(N,B,C)
      DO 20 I=1,N
        B(I)=0.0
  20  CONTINUE
      RETURN
      END 
      SUBROUTINE TESCG(N,B,C,APR,V,NOUT)
C  TESTS RESULTS OF A24CG
C  B= PROPOSED SOLUTION
C  N= ORDER OF PROBLEM
C  C= GIVEN RHS OF EQUATIONS
C  APR = NAME OF MATRIX*VECTOR ROUTINE DESCRIBING COEFFICIENT MATRIX
C  V= WORKING VECTOR OF LENGTH N
C  NOUT= PRINTER CHANNEL
      INTEGER N,I,NOUT
      REAL B(N),C(N),V(N),S,T 
      CALL APR(N,B,V)
      S=0.0
      T=0.0
      DO 10 I=1,N
        V(I)=V(I)-C(I)
        T=T+ABS(B(I)-1.0)
        S=S+ABS(V(I))
        WRITE(NOUT,950)I,B(I),V(I)
 950    FORMAT('  B(',I3,')=',1PE16.8,'  RESIDUAL=',0PE16.8)
  10  CONTINUE
      T=T/N
      S=S/N
      WRITE(NOUT,951)S,T
 951  FORMAT('0MEAN ABSOLUTE -- RESIDUAL=',1PE16.8,' -- ERROR=',E16.8)
      RETURN
      END 
      SUBROUTINE A24CG(N,B,C,TOL,G,IPR,APR,V,T,IMULT)
C  ALGORITHM 24 CONJUGATE GRADIENT SOLUTION OF LINEAR EQUATIONS
C    HAVING POSITIVE DEFINITE COEFFICIENT MATRIX
C   J.C.NASH  FEBRUARY 1980
C  N     = ORDER OF PROBLEM AND LENGTH OF VECTORS B,C,G,V,T 
C  B     = INITIAL GUESS FOR SOLUTION (INPUT) MAY BE NULL
C        = SOLUTION (OUTPUT)
C  C     = RIGHT HAND SIDE OF EQUATIONS 
C  TOL   = CONVERGENCE TOLERANCE ON SIZE OF RESIDUAL SUM OF SQUARES
C  G     = VECTOR OF APPROXIMATE RESIDUALS (OUTPUT)
C  IPR   = PRINTER CHANNEL. IF IPR.GT.0 PRINT TO CHANNEL IPR
C  APR   = NAME OF SUBROUTINE TO PUT A*T IN V
C         CALLING SEQUENCE:   CALL APR(N,T,V)
C  V,T   = WORK ARRAYS OF N ELEMENTS EACH
C  IMULT = NO. OF MATRIX MULTIPLICATIONS ALLOWED (INPUT) OR USED (OUT)
C  STEP 0 
      INTEGER N,IPR,IMULT,LIMIT,I,ITN,COUNT
      REAL TOL,G2,G2L,K,T2,B(N),C(N),G(N),T(N),V(N)
      LIMIT=IMULT
      IMULT=0
C  STEP 1 
   5  IMULT=IMULT+1 
      IF(IMULT.GT.LIMIT)RETURN
      CALL APR(N,B,G)
      DO 10 I=1,N
        G(I)=G(I)-C(I)
  10  CONTINUE
C  STEP 2 
  20  G2=0.0
      DO 25 I=1,N
        G2=G2+G(I)**2
        T(I)=-G(I)
  25  CONTINUE
      IF(IPR.GT.0)WRITE(IPR,950)IMULT,G2
 950  FORMAT(13H AFTER STEP 2,I4,15H PRODUCTS, RSS=,1PE16.8)
C  STEP 3 
      IF(G2.LT.TOL)RETURN
C  STEP 4 
      DO 110 ITN=1,N
C  STEP 5 
        IMULT=IMULT+1
        IF(IMULT.GT.LIMIT)RETURN
      CALL APR(N,T,V)
C  STEP 6 
        T2=0.0
        DO 60 I=1,N 
          T2=T2+T(I)*V(I)
  60    CONTINUE
C  STEP 7 
        K=G2/T2
        G2L=G2
C  STEP 8 
        G2=0.0
        COUNT=0
        DO 80 I=1,N 
          G(I)=G(I)+K*V(I)
          T2=B(I)
          B(I)=T2+K*T(I)
          IF(T2.EQ.B(I))COUNT=COUNT+1
          G2=G2+G(I)**2
  80    CONTINUE
C  STEP 9 
      IF(COUNT.EQ.N)GOTO 5
C     IF(COUNT.EQ.N)GOTO 20
C       IF(G2.LT.TOL)RETURN
C   USING ALTERNATIVE STEP 9 ROUTE
        IF(G2.LT.TOL)GOTO 5
      IF(IPR.GT.0)WRITE(IPR,951)IMULT,G2
 951  FORMAT(13H AFTER STEP 9,I4,15H PRODUCTS, RSS=,1PE16.8)
C  STEP 10
        T2=G2/G2L
        DO 100 I=1,N
          T(I)=T2*T(I)-G(I)
 100    CONTINUE
C  STEP 11
 110  CONTINUE
C  STEP 12
      GOTO 5
C  NO STEP 13 SINCE RETURN USED
      END 
      SUBROUTINE FRANK(N,T,V) 
      INTEGER N,I,J 
      REAL T(N),V(N),S
      DO 10 I=1,N
        S=0.0
        DO 5 J=1,N
          S=S+AMIN0(I,J)*T(J) 
   5    CONTINUE
        V(I)=S
  10  CONTINUE
      RETURN
      END 
   10   50 0.000001 
   10   10 0.000001 
   10    5 0.000001 
   50   51 0.000001 
   50   10 0.000001 
   50  100 0.001