# ssiccg.f

          SUBROUTINE SSICCG (N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
         +   ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW)
    C***BEGIN PROLOGUE  SSICCG
    C***PURPOSE  Incomplete Cholesky Conjugate Gradient Sparse Ax=b Solver.
    C            Routine to solve a symmetric positive definite linear
    C            system  Ax = b  using the incomplete Cholesky
    C            Preconditioned Conjugate Gradient method.
    C***LIBRARY   SLATEC (SLAP)
    C***CATEGORY  D2B4
    C***TYPE      SINGLE PRECISION (SSICCG-S, DSICCG-D)
    C***KEYWORDS  INCOMPLETE CHOLESKY, ITERATIVE PRECONDITION, SLAP, SPARSE,
    C             SYMMETRIC LINEAR SYSTEM
    C***AUTHOR  Greenbaum, Anne, (Courant Institute)
    C           Seager, Mark K., (LLNL)
    C             Lawrence Livermore National Laboratory
    C             PO BOX 808, L-60
    C             Livermore, CA 94550 (510) 423-3141
    C             seager@llnl.gov
    C***DESCRIPTION
    C
    C *Usage:
    C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
    C     INTEGER ITER, IERR, IUNIT, LENW, IWORK(NL+2*N+1), LENIW
    C     REAL B(N), X(N), A(NELT), TOL, ERR, RWORK(NL+5*N)
    C
    C     CALL SSICCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
    C    $     ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
    C
    C *Arguments:
    C N      :IN       Integer.
    C         Order of the Matrix.
    C B      :IN       Real B(N).
    C         Right-hand side vector.
    C X      :INOUT    Real X(N).
    C         On input X is your initial guess for solution vector.
    C         On output X is the final approximate solution.
    C NELT   :IN       Integer.
    C         Number of Non-Zeros stored in A.
    C IA     :INOUT    Integer IA(NELT).
    C JA     :INOUT    Integer JA(NELT).
    C A      :INOUT    Real A(NELT).
    C         These arrays should hold the matrix A in either the SLAP
    C         Triad format or the SLAP Column format.  See "Description",
    C         below.  If the SLAP Triad format is chosen it is changed
    C         internally to the SLAP Column format.
    C ISYM   :IN       Integer.
    C         Flag to indicate symmetric storage format.
    C         If ISYM=0, all non-zero entries of the matrix are stored.
    C         If ISYM=1, the matrix is symmetric, and only the upper
    C         or lower triangle of the matrix is stored.
    C ITOL   :IN       Integer.
    C         Flag to indicate type of convergence criterion.
    C         If ITOL=1, iteration stops when the 2-norm of the residual
    C         divided by the 2-norm of the right-hand side is less than TOL.
    C         If ITOL=2, iteration stops when the 2-norm of M-inv times the
    C         residual divided by the 2-norm of M-inv times the right hand
    C         side is less than TOL, where M-inv is the inverse of the
    C         diagonal of A.
    C         ITOL=11 is often useful for checking and comparing different
    C         routines.  For this case, the user must supply the "exact"
    C         solution or a very accurate approximation (one with an error
    C         much less than TOL) through a common block,
    C             COMMON /SSLBLK/ SOLN( )
    C         If ITOL=11, iteration stops when the 2-norm of the difference
    C         between the iterative approximation and the user-supplied
    C         solution divided by the 2-norm of the user-supplied solution
    C         is less than TOL.  Note that this requires the user to set up
    C         the "COMMON /SSLBLK/ SOLN(LENGTH)" in the calling routine.
    C         The routine with this declaration should be loaded before the
    C         stop test so that the correct length is used by the loader.
    C         This procedure is not standard Fortran and may not work
    C         correctly on your system (although it has worked on every
    C         system the authors have tried).  If ITOL is not 11 then this
    C         common block is indeed standard Fortran.
    C TOL    :INOUT    Real.
    C         Convergence criterion, as described above.  (Reset if IERR=4.)
    C ITMAX  :IN       Integer.
    C         Maximum number of iterations.
    C ITER   :OUT      Integer.
    C         Number of iterations required to reach convergence, or
    C         ITMAX+1 if convergence criterion could not be achieved in
    C         ITMAX iterations.
    C ERR    :OUT      Real.
    C         Error estimate of error in final approximate solution, as
    C         defined by ITOL.
    C IERR   :OUT      Integer.
    C         Return error flag.
    C           IERR = 0 => All went well.
    C           IERR = 1 => Insufficient space allocated for WORK or IWORK.
    C           IERR = 2 => Method failed to converge in ITMAX steps.
    C           IERR = 3 => Error in user input.
    C                       Check input values of N, ITOL.
    C           IERR = 4 => User error tolerance set too tight.
    C                       Reset to 500*R1MACH(3).  Iteration proceeded.
    C           IERR = 5 => Preconditioning matrix, M, is not positive
    C                       definite.  (r,z) < 0.
    C           IERR = 6 => Matrix A is not positive definite.  (p,Ap) < 0.
    C           IERR = 7 => Incomplete factorization broke down and was
    C                       fudged.  Resulting preconditioning may be less
    C                       than the best.
    C IUNIT  :IN       Integer.
    C         Unit number on which to write the error at each iteration,
    C         if this is desired for monitoring convergence.  If unit
    C         number is 0, no writing will occur.
    C RWORK  :WORK     Real RWORK(LENW).
    C         Real array used for workspace.
    C LENW   :IN       Integer.
    C         Length of the real workspace, RWORK.  LENW >= NL+5*N.
    C         NL is the number of non-zeros in the lower triangle of the
    C         matrix (including the diagonal).
    C IWORK  :WORK     Integer IWORK(LENIW).
    C         Integer array used for workspace.
    C         Upon return the following locations of IWORK hold information
    C         which may be of use to the user:
    C         IWORK(9)  Amount of Integer workspace actually used.
    C         IWORK(10) Amount of Real workspace actually used.
    C LENIW  :IN       Integer.
    C         Length of the integer workspace, IWORK.  LENIW >= NL+N+11.
    C         NL is the number of non-zeros in the lower triangle of the
    C         matrix (including the diagonal).
    C
    C *Description:
    C       This routine  performs  preconditioned  conjugate   gradient
    C       method on the   symmetric positive  definite  linear  system
    C       Ax=b.   The preconditioner  is  the incomplete Cholesky (IC)
    C       factorization of the matrix A.  See  SSICS for details about
    C       the incomplete   factorization algorithm.  One   should note
    C       here however, that the  IC factorization is a  slow  process
    C       and  that  one should   save  factorizations  for  reuse, if
    C       possible.  The   MSOLVE operation (handled  in  SSLLTI) does
    C       vectorize on machines  with  hardware  gather/scatter and is
    C       quite fast.
    C
    C       The Sparse Linear Algebra Package (SLAP) utilizes two matrix
    C       data structures: 1) the  SLAP Triad  format or  2)  the SLAP
    C       Column format.  The user can hand this routine either of the
    C       of these data structures and SLAP  will figure out  which on
    C       is being used and act accordingly.
    C
    C       =================== S L A P Triad format ===================
    C
    C       This routine requires that the  matrix A be   stored in  the
    C       SLAP  Triad format.  In  this format only the non-zeros  are
    C       stored.  They may appear in  *ANY* order.  The user supplies
    C       three arrays of  length NELT, where  NELT is  the number  of
    C       non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)).  For
    C       each non-zero the user puts the row and column index of that
    C       matrix element  in the IA and  JA arrays.  The  value of the
    C       non-zero   matrix  element is  placed  in  the corresponding
    C       location of the A array.   This is  an  extremely  easy data
    C       structure to generate.  On  the  other hand it   is  not too
    C       efficient on vector computers for  the iterative solution of
    C       linear systems.  Hence,   SLAP changes   this  input    data
    C       structure to the SLAP Column format  for  the iteration (but
    C       does not change it back).
    C
    C       Here is an example of the  SLAP Triad   storage format for a
    C       5x5 Matrix.  Recall that the entries may appear in any order.
    C
    C           5x5 Matrix      SLAP Triad format for 5x5 matrix on left.
    C                              1  2  3  4  5  6  7  8  9 10 11
    C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
    C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
    C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
    C       | 0  0  0 44  0|
    C       |51  0 53  0 55|
    C
    C       =================== S L A P Column format ==================
    C
    C       This routine  requires that  the matrix A  be stored in  the
    C       SLAP Column format.  In this format the non-zeros are stored
    C       counting down columns (except for  the diagonal entry, which
    C       must appear first in each  "column")  and are stored  in the
    C       real array A.  In other words, for each column in the matrix
    C       put the diagonal entry in A.  Then put in the other non-zero
    C       elements going down   the  column (except  the diagonal)  in
    C       order.  The IA array holds the row  index for each non-zero.
    C       The JA array holds the offsets into the IA, A arrays for the
    C       beginning of   each    column.    That  is,    IA(JA(ICOL)),
    C       A(JA(ICOL)) points to the beginning of the ICOL-th column in
    C       IA and  A.  IA(JA(ICOL+1)-1),  A(JA(ICOL+1)-1) points to the
    C       end  of   the ICOL-th  column.  Note   that  we  always have
    C       JA(N+1) = NELT+1, where  N  is the number of columns in  the
    C       matrix and  NELT   is the number of non-zeros in the matrix.
    C
    C       Here is an example of the  SLAP Column  storage format for a
    C       5x5 Matrix (in the A and IA arrays '|'  denotes the end of a
    C       column):
    C
    C           5x5 Matrix      SLAP Column format for 5x5 matrix on left.
    C                              1  2  3    4  5    6  7    8    9 10 11
    C       |11 12  0  0 15|   A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
    C       |21 22  0  0  0|  IA:  1  2  5 |  2  1 |  3  5 |  4 |  5  1  3
    C       | 0  0 33  0 35|  JA:  1  4  6    8  9   12
    C       | 0  0  0 44  0|
    C       |51  0 53  0 55|
    C
    C *Side Effects:
    C       The SLAP Triad format (IA, JA, A) is modified internally to be
    C       the SLAP Column format.  See above.
    C
    C *Cautions:
    C     This routine will attempt to write to the Fortran logical output
    C     unit IUNIT, if IUNIT .ne. 0.  Thus, the user must make sure that
    C     this logical unit is attached to a file or terminal before calling
    C     this routine with a non-zero value for IUNIT.  This routine does
    C     not check for the validity of a non-zero IUNIT unit number.
    C
    C***SEE ALSO  SCG, SSLLTI
    C***REFERENCES  1. Louis Hageman and David Young, Applied Iterative
    C                  Methods, Academic Press, New York, 1981.
    C               2. Concus, Golub and O'Leary, A Generalized Conjugate
    C                  Gradient Method for the Numerical Solution of
    C                  Elliptic Partial Differential Equations, in Sparse
    C                  Matrix Computations, Bunch and Rose, Eds., Academic
    C                  Press, New York, 1979.
    C***ROUTINES CALLED  SCG, SCHKW, SS2Y, SSICS, SSLLTI, SSMV, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   871119  DATE WRITTEN
    C   881213  Previous REVISION DATE
    C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
    C   890921  Removed TeX from comments.  (FNF)
    C   890922  Numerous changes to prologue to make closer to SLATEC
    C           standard.  (FNF)
    C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
    C   900805  Changed XERRWV calls to calls to XERMSG.  (RWC)
    C   910411  Prologue converted to Version 4.0 format.  (BAB)
    C   920407  COMMON BLOCK renamed SSLBLK.  (WRB)
    C   920511  Added complete declaration section.  (WRB)
    C   920929  Corrected format of references.  (FNF)
    C   921019  Corrected NEL to NL.  (FNF)
    C***END PROLOGUE  SSICCG