*DECK SLPDOC SUBROUTINE SLPDOC C***BEGIN PROLOGUE SLPDOC C***PURPOSE Sparse Linear Algebra Package Version 2.0.2 Documentation. C Routines to solve large sparse symmetric and nonsymmetric C positive definite linear systems, Ax = b, using precondi- C tioned iterative methods. C***LIBRARY SLATEC (SLAP) C***CATEGORY D2A4, D2B4, Z C***TYPE SINGLE PRECISION (SLPDOC-S, DLPDOC-D) C***KEYWORDS BICONJUGATE GRADIENT SQUARED, DOCUMENTATION, C GENERALIZED MINIMUM RESIDUAL, ITERATIVE IMPROVEMENT, C NORMAL EQUATIONS, ORTHOMIN, C PRECONDITIONED CONJUGATE GRADIENT, SLAP, C SPARSE ITERATIVE METHODS C***AUTHOR Seager, Mark. K., (LLNL) C User Systems Division C Lawrence Livermore National Laboratory C PO BOX 808, L-60 C Livermore, CA 94550 C (FTS) 543-3141, (510) 423-3141 C seager@llnl.gov C***DESCRIPTION C The C Sparse Linear Algebra Package C C @@@@@@@ @ @@@ @@@@@@@@ C @ @ @ @ @ @ @ C @ @ @ @ @ @ C @@@@@@@ @ @ @ @@@@@@@@ C @ @ @@@@@@@@@ @ C @ @ @ @ @ @ C @@@@@@@ @@@@@@@@@ @ @ @ C C @ @ @@@@@@@ @@@@@ C @ @ @ @ @ @@ C @ @ @@@@@@@ @ @@ @ @ @ @ C @ @ @ @ @@ @ @@@@@@ @ @ @ C @ @ @@@@@@@@@ @ @ @ @ @ C @ @ @ @ @ @@@ @@ @ C @@@ @@@@@@@ @ @@@@@@@@@ @@@ @@@@@ C C C ================================================================= C ========================== Introduction ========================= C ================================================================= C This package was originally derived from a set of iterative C routines written by Anne Greenbaum, as announced in "Routines C for Solving Large Sparse Linear Systems", Tentacle, Lawrence C Livermore National Laboratory, Livermore Computing Center C (January 1986), pp 15-21. C C This document contains the specifications for the SLAP Version C 2.0 package, a Fortran 77 package for the solution of large C sparse linear systems, Ax = b, via preconditioned iterative C methods. Included in this package are "core" routines to do C Iterative Refinement (Jacobi's method), Conjugate Gradient, C Conjugate Gradient on the normal equations, AA'y = b, (where x = C A'y and A' denotes the transpose of A), BiConjugate Gradient, C BiConjugate Gradient Squared, Orthomin and Generalized Minimum C Residual Iteration. These "core" routines do not require a C "fixed" data structure for storing the matrix A and the C preconditioning matrix M. The user is free to choose any C structure that facilitates efficient solution of the problem at C hand. The drawback to this approach is that the user must also C supply at least two routines (MATVEC and MSOLVE, say). MATVEC C must calculate, y = Ax, given x and the user's data structure for C A. MSOLVE must solve, r = Mz, for z (*NOT* r) given r and the C user's data structure for M (or its inverse). The user should C choose M so that inv(M)*A is approximately the identity and the C solution step r = Mz is "easy" to solve. For some of the "core" C routines (Orthomin, BiConjugate Gradient and Conjugate Gradient C on the normal equations) the user must also supply a matrix C transpose times vector routine (MTTVEC, say) and (possibly, C depending on the "core" method) a routine that solves the C transpose of the preconditioning step (MTSOLV, say). C Specifically, MTTVEC is a routine which calculates y = A'x, given C x and the user's data structure for A (A' is the transpose of A). C MTSOLV is a routine which solves the system r = M'z for z given r C and the user's data structure for M. C C This process of writing the matrix vector operations can be time C consuming and error prone. To alleviate these problems we have C written drivers for the "core" methods that assume the user C supplies one of two specific data structures (SLAP Triad and SLAP C Column format), see below. Utilizing these data structures we C have augmented each "core" method with two preconditioners: C Diagonal Scaling and Incomplete Factorization. Diagonal scaling C is easy to implement, vectorizes very well and for problems that C are not too ill-conditioned reduces the number of iterations C enough to warrant its use. On the other hand, an Incomplete C factorization (Incomplete Cholesky for symmetric systems and C Incomplete LU for nonsymmetric systems) may take much longer to C calculate, but it reduces the iteration count (for most problems) C significantly. Our implementations of IC and ILU vectorize for C machines with hardware gather scatter, but the vector lengths can C be quite short if the number of non-zeros in a column is not C large. C C ================================================================= C ==================== Supplied Data Structures =================== C ================================================================= C The following describes the data structures supplied with the C package: SLAP Triad and Column formats. C C ====================== S L A P Triad format ===================== C C In the SLAP Triad format only the non-zeros are stored. They may C appear in *ANY* order. The user supplies three arrays of length C NELT, where NELT is the number of non-zeros in the matrix: C (IA(NELT), JA(NELT), A(NELT)). If the matrix is symmetric then C one need only store the lower triangle (including the diagonal) C and NELT would be the corresponding number of non-zeros stored. C For 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 location C of the A array. This is an extremely easy data structure to C generate. On the other hand, it is not very efficient on vector C computers for the iterative solution of linear systems. Hence, C SLAP changes this input data structure to the SLAP Column format C for the iteration (but does not change it back). C C Here is an example of the SLAP Triad storage format for a C nonsymmetric 5x5 Matrix. NELT=11. Recall that the entries may C 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 In the SLAP Column format the non-zeros are stored counting down C columns (except for the diagonal entry, which must appear first C in each "column") and are stored in the real array A. In other C words, for each column in the matrix first put the diagonal entry C in A. Then put in the other non-zero elements going down the C column (except the diagonal) in order. The IA array holds the C row index for each non-zero. The JA array holds the offsets into C the IA, A arrays for the beginning of each column. That is, C IA(JA(ICOL)), A(JA(ICOL)) are the first elements of the ICOL-th C column in IA and A. IA(JA(ICOL+1)-1), A(JA(ICOL+1)-1) are the C last elements 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 matrix C and NELT is the number of non-zeros in the matrix. If the matrix C is symmetric one need only store the lower triangle (including C the diagonal) and NELT would be the corresponding number of C non-zeros stored. C C Here is an example of the SLAP Column storage format for a C nonsymmetric 5x5 Matrix (in the A and IA arrays '|' denotes the C end of a 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 ================================================================= C ====================== Which Method To Use ====================== C ================================================================= C C BACKGROUND C In solving a large sparse linear system Ax = b using an iterative C method, it is not necessary to actually store the matrix A. C Rather, what is needed is a procedure for multiplying the matrix C A times a given vector y to obtain the matrix-vector product, Ay. C SLAP has been written to take advantage of this fact. The higher C level routines in the package require storage only of the non-zero C elements of A (and their positions), and even this can be C avoided, if the user writes his own subroutine for multiplying C the matrix times a vector and calls the lower-level iterative C routines in the package. C C If the matrix A is ill-conditioned, then most iterative methods C will be slow to converge (if they converge at all!). To improve C the convergence rate, one may use a "matrix splitting," or, C "preconditioning matrix," say, M. It is then necessary to solve, C at each iteration, a linear system with coefficient matrix M. A C good preconditioner M should have two properties: (1) M should C "approximate" A, in the sense that the matrix inv(M)*A (or some C variant thereof) is better conditioned than the original matrix C A; and (2) linear systems with coefficient matrix M should be C much easier to solve than the original system with coefficient C matrix A. Preconditioning routines in the SLAP package are C separate from the iterative routines, so that any of the C preconditioners provided in the package, or one that the user C codes himself, can be used with any of the iterative routines. C C CHOICE OF PRECONDITIONER C If you willing to live with either the SLAP Triad or Column C matrix data structure you can then choose one of two types of C preconditioners to use: diagonal scaling or incomplete C factorization. To choose between these two methods requires C knowing something about the computer you're going to run these C codes on and how well incomplete factorization approximates the C inverse of your matrix. C C Let us suppose you have a scalar machine. Then, unless the C incomplete factorization is very, very poor this is *GENERALLY* C the method to choose. It will reduce the number of iterations C significantly and is not all that expensive to compute. So if C you have just one linear system to solve and "just want to get C the job done" then try incomplete factorization first. If you C are thinking of integrating some SLAP iterative method into your C favorite "production code" then try incomplete factorization C first, but also check to see that diagonal scaling is indeed C slower for a large sample of test problems. C C Let us now suppose you have a vector computer with hardware C gather/scatter support (Cray X-MP, Y-MP, SCS-40 or Cyber 205, ETA C 10, ETA Piper, Convex C-1, etc.). Then it is much harder to C choose between the two methods. The versions of incomplete C factorization in SLAP do in fact vectorize, but have short vector C lengths and the factorization step is relatively more expensive. C Hence, for most problems (i.e., unless your problem is ill C conditioned, sic!) diagonal scaling is faster, with its very C fast set up time and vectorized (with long vectors) C preconditioning step (even though it may take more iterations). C If you have several systems (or right hand sides) to solve that C can utilize the same preconditioner then the cost of the C incomplete factorization can be amortized over these several C solutions. This situation gives more advantage to the incomplete C factorization methods. If you have a vector machine without C hardware gather/scatter (Cray 1, Cray 2 & Cray 3) then the C advantages for incomplete factorization are even less. C C If you're trying to shoehorn SLAP into your favorite "production C code" and can not easily generate either the SLAP Triad or Column C format then you are left to your own devices in terms of C preconditioning. Also, you may find that the preconditioners C supplied with SLAP are not sufficient for your problem. In this C situation we would recommend that you talk with a numerical C analyst versed in iterative methods about writing other C preconditioning subroutines (e.g., polynomial preconditioning, C shifted incomplete factorization, SOR or SSOR iteration). You C can always "roll your own" by using the "core" iterative methods C and supplying your own MSOLVE and MATVEC (and possibly MTSOLV and C MTTVEC) routines. C C SYMMETRIC SYSTEMS C If your matrix is symmetric then you would want to use one of the C symmetric system solvers. If your system is also positive C definite, (Ax,x) (Ax dot product with x) is positive for all C non-zero vectors x, then use Conjugate Gradient (SCG, SSDCG, C SSICSG). If you're not sure it's SPD (symmetric and Positive C Definite) then try SCG anyway and if it works, fine. If you're C sure your matrix is not positive definite then you may want to C try the iterative refinement methods (SIR) or the GMRES code C (SGMRES) if SIR converges too slowly. C C NONSYMMETRIC SYSTEMS C This is currently an area of active research in numerical C analysis and there are new strategies being developed. C Consequently take the following advice with a grain of salt. If C you matrix is positive definite, (Ax,x) (Ax dot product with x C is positive for all non-zero vectors x), then you can use any of C the methods for nonsymmetric systems (Orthomin, GMRES, C BiConjugate Gradient, BiConjugate Gradient Squared and Conjugate C Gradient applied to the normal equations). If your system is not C too ill conditioned then try BiConjugate Gradient Squared (BCGS) C or GMRES (SGMRES). Both of these methods converge very quickly C and do not require A' or M' (' denotes transpose) information. C SGMRES does require some additional storage, though. If the C system is very ill conditioned or nearly positive indefinite C ((Ax,x) is positive, but may be very small), then GMRES should C be the first choice, but try the other methods if you have to C fine tune the solution process for a "production code". If you C have a great preconditioner for the normal equations (i.e., M is C an approximation to the inverse of AA' rather than just A) then C this is not a bad route to travel. Old wisdom would say that the C normal equations are a disaster (since it squares the condition C number of the system and SCG convergence is linked to this number C of infamy), but some preconditioners (like incomplete C factorization) can reduce the condition number back below that of C the original system. C C ================================================================= C ======================= Naming Conventions ====================== C ================================================================= C SLAP iterative methods, matrix vector and preconditioner C calculation routines follow a naming convention which, when C understood, allows one to determine the iterative method and data C structure(s) used. The subroutine naming convention takes the C following form: C P[S][M]DESC C where C P stands for the precision (or data type) of the routine and C is required in all names, C S denotes whether or not the routine requires the SLAP Triad C or Column format (it does if the second letter of the name C is S and does not otherwise), C M stands for the type of preconditioner used (only appears C in drivers for "core" routines), and C DESC is some number of letters describing the method or purpose C of the routine. The following is a list of the "DESC" C fields for iterative methods and their meaning: C BCG,BC: BiConjugate Gradient C CG: Conjugate Gradient C CGN,CN: Conjugate Gradient on the Normal equations C CGS,CS: biConjugate Gradient Squared C GMRES,GMR,GM: Generalized Minimum RESidual C IR,R: Iterative Refinement C JAC: JACobi's method C GS: Gauss-Seidel C OMN,OM: OrthoMiN C C In the single precision version of SLAP, all routine names start C with an S. The brackets around the S and M designate that these C fields are optional. C C Here are some examples of the routines: C 1) SBCG: Single precision BiConjugate Gradient "core" routine. C One can deduce that this is a "core" routine, because the S and C M fields are missing and BiConjugate Gradient is an iterative C method. C 2) SSDBCG: Single precision, SLAP data structure BCG with Diagonal C scaling. C 3) SSLUBC: Single precision, SLAP data structure BCG with incom- C plete LU factorization as the preconditioning. C 4) SCG: Single precision Conjugate Gradient "core" routine. C 5) SSDCG: Single precision, SLAP data structure Conjugate Gradient C with Diagonal scaling. C 6) SSICCG: Single precision, SLAP data structure Conjugate Gra- C dient with Incomplete Cholesky factorization preconditioning. C C C ================================================================= C ===================== USER CALLABLE ROUTINES ==================== C ================================================================= C The following is a list of the "user callable" SLAP routines and C their one line descriptions. The headers denote the file names C where the routines can be found, as distributed for UNIX systems. C C Note: Each core routine, SXXX, has a corresponding stop routine, C ISSXXX. If the stop routine does not have the specific stop C test the user requires (e.g., weighted infinity norm), then C the user should modify the source for ISSXXX accordingly. C C ============================= sir.f ============================= C SIR: Preconditioned Iterative Refinement Sparse Ax = b Solver. C SSJAC: Jacobi's Method Iterative Sparse Ax = b Solver. C SSGS: Gauss-Seidel Method Iterative Sparse Ax = b Solver. C SSILUR: Incomplete LU Iterative Refinement Sparse Ax = b Solver. C C ============================= scg.f ============================= C SCG: Preconditioned Conjugate Gradient Sparse Ax=b Solver. C SSDCG: Diagonally Scaled Conjugate Gradient Sparse Ax=b Solver. C SSICCG: Incomplete Cholesky Conjugate Gradient Sparse Ax=b Solver. C C ============================= scgn.f ============================ C SCGN: Preconditioned CG Sparse Ax=b Solver for Normal Equations. C SSDCGN: Diagonally Scaled CG Sparse Ax=b Solver for Normal Eqn's. C SSLUCN: Incomplete LU CG Sparse Ax=b Solver for Normal Equations. C C ============================= sbcg.f ============================ C SBCG: Preconditioned BiConjugate Gradient Sparse Ax = b Solver. C SSDBCG: Diagonally Scaled BiConjugate Gradient Sparse Ax=b Solver. C SSLUBC: Incomplete LU BiConjugate Gradient Sparse Ax=b Solver. C C ============================= scgs.f ============================ C SCGS: Preconditioned BiConjugate Gradient Squared Ax=b Solver. C SSDCGS: Diagonally Scaled CGS Sparse Ax=b Solver. C SSLUCS: Incomplete LU BiConjugate Gradient Squared Ax=b Solver. C C ============================= somn.f ============================ C SOMN: Preconditioned Orthomin Sparse Iterative Ax=b Solver. C SSDOMN: Diagonally Scaled Orthomin Sparse Iterative Ax=b Solver. C SSLUOM: Incomplete LU Orthomin Sparse Iterative Ax=b Solver. C C ============================ sgmres.f =========================== C SGMRES: Preconditioned GMRES Iterative Sparse Ax=b Solver. C SSDGMR: Diagonally Scaled GMRES Iterative Sparse Ax=b Solver. C SSLUGM: Incomplete LU GMRES Iterative Sparse Ax=b Solver. C C ============================ smset.f ============================ C The following routines are used to set up preconditioners. C C SSDS: Diagonal Scaling Preconditioner SLAP Set Up. C SSDSCL: Diagonally Scales/Unscales a SLAP Column Matrix. C SSD2S: Diagonal Scaling Preconditioner SLAP Normal Eqns Set Up. C SS2LT: Lower Triangle Preconditioner SLAP Set Up. C SSICS: Incomplete Cholesky Decomp. Preconditioner SLAP Set Up. C SSILUS: Incomplete LU Decomposition Preconditioner SLAP Set Up. C C ============================ smvops.f =========================== C Most of the incomplete factorization (LL' and LDU) solvers C in this file require an intermediate routine to translate C from the SLAP MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, C IWORK) calling convention to the calling sequence required C by the solve routine. This generally is accomplished by C fishing out pointers to the preconditioner (stored in RWORK) C from the IWORK array and then making a call to the routine C that actually does the backsolve. C C SSMV: SLAP Column Format Sparse Matrix Vector Product. C SSMTV: SLAP Column Format Sparse Matrix (transpose) Vector Prod. C SSDI: Diagonal Matrix Vector Multiply. C SSLI: SLAP MSOLVE for Lower Triangle Matrix (set up for SSLI2). C SSLI2: Lower Triangle Matrix Backsolve. C SSLLTI: SLAP MSOLVE for LDL' (IC) Fact. (set up for SLLTI2). C SLLTI2: Backsolve routine for LDL' Factorization. C SSLUI: SLAP MSOLVE for LDU Factorization (set up for SSLUI2). C SSLUI2: SLAP Backsolve for LDU Factorization. C SSLUTI: SLAP MTSOLV for LDU Factorization (set up for SSLUI4). C SSLUI4: SLAP Backsolve for LDU Factorization. C SSMMTI: SLAP MSOLVE for LDU Fact of Normal Eq (set up for SSMMI2). C SSMMI2: SLAP Backsolve for LDU Factorization of Normal Equations. C C =========================== slaputil.f ========================== C The following utility routines are useful additions to SLAP. C C SBHIN: Read Sparse Linear System in the Boeing/Harwell Format. C SCHKW: SLAP WORK/IWORK Array Bounds Checker. C SCPPLT: Printer Plot of SLAP Column Format Matrix. C SS2Y: SLAP Triad to SLAP Column Format Converter. C QS2I1R: Quick Sort Integer array, moving integer and real arrays. C (Used by SS2Y.) C STIN: Read in SLAP Triad Format Linear System. C STOUT: Write out SLAP Triad Format Linear System. C C C***REFERENCES 1. Mark K. Seager, A SLAP for the Masses, in C G. F. Carey, Ed., Parallel Supercomputing: Methods, C Algorithms and Applications, Wiley, 1989, pp.135-155. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 880715 DATE WRITTEN C 890404 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 -----( This produced Version 2.0.1. )----- C 891003 Rearranged list of user callable routines to agree with C order in source deck. (FNF) C 891004 Updated reference. C 910411 Prologue converted to Version 4.0 format. (BAB) C -----( This produced Version 2.0.2. )----- C 910506 Minor improvements to prologue. (FNF) C 920511 Added complete declaration section. (WRB) C 920929 Corrected format of reference. (FNF) C 921019 Improved one-line descriptions, reordering some. (FNF) C***END PROLOGUE SLPDOC C***FIRST EXECUTABLE STATEMENT SLPDOC C C This is a *DUMMY* subroutine and should never be called. C RETURN C------------- LAST LINE OF SLPDOC FOLLOWS ----------------------------- END