SUBROUTINE DBCG(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MTTVEC, & MSOLVE, MTSOLV, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, & R, Z, P, RR, ZZ, PP, DZ, RWORK, IWORK) !***BEGIN PROLOGUE DBCG !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2A4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(DBCG-D), ! Non-Symmetric Linear system, Sparse, ! Iterative Precondition, BiConjugate Gradient !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Preconditioned BiConjugate Gradient Sparse Ax=b solver. ! Routine to solve a Non-Symmetric linear system Ax = b ! using the Preconditioned BiConjugate Gradient method. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, IWORK(USER DEFINABLE) ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N), P(N) ! DOUBLE PRECISION RR(N), ZZ(N), PP(N), DZ(N) ! DOUBLE PRECISION RWORK(USER DEFINABLE) ! EXTERNAL MATVEC, MTTVEC, MSOLVE, MTSOLV ! ! CALL DBCG(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MTTVEC, ! $ MSOLVE, MTSOLV, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, ! $ R, Z, P, RR, ZZ, PP, DZ, RWORK, IWORK) ! ! *Arguments: ! N :IN Integer ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :IN Integer IA(NELT). ! JA :IN Integer JA(NELT). ! A :IN Double Precision A(NELT). ! These arrays contain the matrix data structure for A. ! It could take any form. See "Description", below for more ! late breaking details... ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! MATVEC :EXT External. ! Name of a routine which performs the matrix vector multiply ! operation Y = A*X given A and X. The name of the MATVEC ! routine must be declared external in the calling program. ! The calling sequence of MATVEC is: ! CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM ) ! Where N is the number of unknowns, Y is the product A*X upon ! return, X is an input vector. NELT, IA, JA, A and ISYM ! define the SLAP matrix data structure: see Description,below. ! MTTVEC :EXT External. ! Name of a routine which performs the matrix transpose vector ! multiply y = A'*X given A and X (where ' denotes transpose). ! The name of the MTTVEC routine must be declared external in ! the calling program. The calling sequence to MTTVEC is the ! same as that for MTTVEC, viz.: ! CALL MTTVEC( N, X, Y, NELT, IA, JA, A, ISYM ) ! Where N is the number of unknowns, Y is the product A'*X ! upon return, X is an input vector. NELT, IA, JA, A and ISYM ! define the SLAP matrix data structure: see Description,below. ! MSOLVE :EXT External. ! Name of a routine which solves a linear system MZ = R for Z ! given R with the preconditioning matrix M (M is supplied via ! RWORK and IWORK arrays). The name of the MSOLVE routine ! must be declared external in the calling program. The ! calling sequence of MSLOVE is: ! CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! Where N is the number of unknowns, R is the right-hand side ! vector, and Z is the solution upon return. NELT, IA, JA, A ! and ISYM define the SLAP matrix data structure: see ! Description, below. RWORK is a double precision array that ! can be used ! to pass necessary preconditioning information and/or ! workspace to MSOLVE. IWORK is an integer work array for the ! same purpose as RWORK. ! MTSOLV :EXT External. T ! Name of a routine which solves a linear system M ZZ = RR for ! ZZ given RR with the preconditioning matrix M (M is supplied ! via RWORK and IWORK arrays). The name of the MTSOLV routine ! must be declared external in the calling program. The call- ! ing sequence to MTSOLV is: ! CALL MTSOLV(N, RR, ZZ, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! Where N is the number of unknowns, RR is the right-hand side ! vector, and ZZ is the solution upon return. NELT, IA, JA, A ! and ISYM define the SLAP matrix data structure: see ! Description, below. RWORK is a double precision array that ! can be used ! to pass necessary preconditioning information and/or ! workspace to MTSOLV. IWORK is an integer work array for the ! same purpose as RWORK. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! R :WORK Double Precision R(N). ! Z :WORK Double Precision Z(N). ! P :WORK Double Precision P(N). ! RR :WORK Double Precision RR(N). ! ZZ :WORK Double Precision ZZ(N). ! PP :WORK Double Precision PP(N). ! DZ :WORK Double Precision DZ(N). ! RWORK :WORK Double Precision RWORK(USER DEFINED). ! Double Precision array that can be used for workspace in ! MSOLVE and MTSOLV. ! IWORK :WORK Integer IWORK(USER DEFINED). ! Integer array that can be used for workspace in MSOLVE ! and MTSOLV. ! ! *Description ! This routine does not care what matrix data structure is ! used for A and M. It simply calls the MATVEC and MSOLVE ! routines, with the arguments as described above. The user ! could write any type of structure and the appropriate MATVEC ! and MSOLVE routines. It is assumed that A is stored in the ! IA, JA, A arrays in some fashion and that M (or INV(M)) is ! stored in IWORK and RWORK in some fashion. The SLAP ! routines SDBCG and DSLUBC are examples of this procedure. ! ! Two examples of matrix data structures are the: 1) SLAP ! Triad format and 2) SLAP Column format. ! ! =================== S L A P Triad format =================== ! In this format only the non-zeros are stored. They may ! appear in *ANY* order. The user supplies three arrays of ! length NELT, where NELT is the number of non-zeros in the ! matrix: (IA(NELT), JA(NELT), A(NELT)). For each non-zero ! the user puts the row and column index of that matrix ! element in the IA and JA arrays. The value of the non-zero ! matrix element is placed in the corresponding location of ! the A array. This is an extremely easy data structure to ! generate. On the other hand it is not too efficient on ! vector computers for the iterative solution of linear ! systems. Hence, SLAP changes this input data structure to ! the SLAP Column format for the iteration (but does not ! change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *See Also: ! SDBCG, DSLUBC !***REFERENCES (NONE) !***ROUTINES CALLED MATVEC, MTTVEC, MSOLVE, MTSOLV, ISDBCG, ! DCOPY, DDOT, DAXPY, D1MACH !***END PROLOGUE DBCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX INTEGER ITER, IERR, IUNIT, IWORK(*) DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N), P(N) DOUBLE PRECISION RR(N), ZZ(N), PP(N), DZ(N), RWORK(*) EXTERNAL MATVEC, MTTVEC, MSOLVE, MTSOLV ! ! Check some of the input data. ! !***FIRST EXECUTABLE STATEMENT DBCG ITER = 0 IERR = 0 IF( N.LT.1 ) THEN IERR = 3 RETURN end if FUZZ = epsilon ( fuzz ) TOLMIN = 500.0*FUZZ FUZZ = FUZZ*FUZZ IF( TOL.LT.TOLMIN ) THEN TOL = TOLMIN IERR = 4 end if ! ! Calculate initial residual and pseudo-residual, and check ! stopping criterion. ! CALL MATVEC(N, X, R, NELT, IA, JA, A, ISYM) R(1:n) = B(1:n) - R(1:n) RR(1:n) = R(1:n) CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) CALL MTSOLV(N, RR, ZZ, NELT, IA, JA, A, ISYM, RWORK, IWORK) IF( ISDBCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, PP, & DZ, RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) & GO TO 200 IF( IERR.NE.0 ) RETURN ! ! iteration loop ! DO 100 K=1,ITMAX ITER = K ! ! Calculate coefficient BK and direction vectors P and PP. ! BKNUM = DDOT(N, Z, 1, RR, 1) IF( ABS(BKNUM).LE.FUZZ ) THEN IERR = 6 RETURN end if IF(ITER .EQ. 1) THEN CALL DCOPY(N, Z, 1, P, 1) CALL DCOPY(N, ZZ, 1, PP, 1) ELSE BK = BKNUM/BKDEN DO 20 I = 1, N P(I) = Z(I) + BK*P(I) PP(I) = ZZ(I) + BK*PP(I) 20 CONTINUE end if BKDEN = BKNUM ! ! Calculate coefficient AK, new iterate X, new resids R and RR, ! and new pseudo-resids Z and ZZ. ! CALL MATVEC(N, P, Z, NELT, IA, JA, A, ISYM) AKDEN = DDOT(N, PP, 1, Z, 1) AK = BKNUM/AKDEN IF( ABS(AKDEN).LE.FUZZ ) THEN IERR = 6 RETURN end if CALL DAXPY(N, AK, P, 1, X, 1) CALL DAXPY(N, -AK, Z, 1, R, 1) CALL MTTVEC(N, PP, ZZ, NELT, IA, JA, A, ISYM) CALL DAXPY(N, -AK, ZZ, 1, RR, 1) CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) CALL MTSOLV(N, RR, ZZ, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! ! Check stopping criterion. ! IF( ISDBCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, & PP, DZ, RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) & GO TO 200 100 CONTINUE ! ! Stopping criterion not satisfied. ! ITER = ITMAX + 1 IERR = 2 200 RETURN END SUBROUTINE DSDBCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) !***BEGIN PROLOGUE DSDBCG !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2A4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(SSDBCG-D), ! Non-Symmetric Linear system, Sparse, ! Iterative Precondition !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Diagonally Scaled BiConjugate Gradient Sparse Ax=b solver. ! Routine to solve a linear system Ax = b using the ! BiConjugate Gradient method with diagonal scaling. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, LENW, IWORK(10), LENIW ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(8*N) ! ! CALL DSDBCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) ! ! *Arguments: ! N :IN Integer ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :INOUT Integer IA(NELT). ! JA :INOUT Integer JA(NELT). ! A :INOUT Double Precision A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See "Description", ! below. If the SLAP Triad format is chosen it is changed ! internally to the SLAP Column format. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! RWORK :WORK Double Precision RWORK(LENW). ! Double Precision array used for workspace. ! LENW :IN Integer. ! Length of the double precision workspace, RWORK. ! LENW >= 8*N. ! IWORK :WORK Integer IWORK(LENIW). ! Used to hold pointers into the RWORK array. ! LENIW :IN Integer. ! Length of the integer workspace, IWORK. LENIW >= 10. ! Upon return the following locations of IWORK hold information ! which may be of use to the user: ! IWORK(9) Amount of Integer workspace actually used. ! IWORK(10) Amount of Double Precision workspace actually used. ! ! *Description: ! This routine performs preconditioned BiConjugate gradient ! method on the Non-Symmetric positive definite linear system ! Ax=b. The preconditioner is M = DIAG(A), the diagonal of the ! matrix A. This is the simplest of preconditioners and ! vectorizes very well. ! ! The Sparse Linear Algebra Package (SLAP) utilizes two matrix ! data structures: 1) the SLAP Triad format or 2) the SLAP ! Column format. The user can hand this routine either of the ! of these data structures and SLAP will figure out which on ! is being used and act accordingly. ! ! =================== S L A P Triad format =================== ! ! This routine requires that the matrix A be stored in the ! SLAP Triad format. In this format only the non-zeros are ! stored. They may appear in *ANY* order. The user supplies ! three arrays of length NELT, where NELT is the number of ! non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For ! each non-zero the user puts the row and column index of that ! matrix element in the IA and JA arrays. The value of the ! non-zero matrix element is placed in the corresponding ! location of the A array. This is an extremely easy data ! structure to generate. On the other hand it is not too ! efficient on vector computers for the iterative solution of ! linear systems. Hence, SLAP changes this input data ! structure to the SLAP Column format for the iteration (but ! does not change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *Side Effects: ! The SLAP Triad format (IA, JA, A) is modified internally to ! be the SLAP Column format. See above. ! ! *See Also: ! DBCG, DLUBCG !***REFERENCES (NONE) !***ROUTINES CALLED DS2Y, DCHKW, DSDS, DBCG !***END PROLOGUE DSDBCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER INTEGER IERR, LENW, IWORK(LENIW), LENIW DOUBLE PRECISION B(N), X(N), A(N), TOL, ERR, RWORK(LENW) EXTERNAL DSMV, DSMTV, DSDI PARAMETER (LOCRB=1, LOCIB=11) ! ! Change the SLAP input matrix IA, JA, A to SLAP-Column format. !***FIRST EXECUTABLE STATEMENT DSDBCG IERR = 0 IF( N.LT.1 .OR. NELT.LT.1 ) THEN IERR = 3 RETURN end if CALL DS2Y( N, NELT, IA, JA, A, ISYM ) ! ! Set up the workspace. Compute the inverse of the ! diagonal of the matrix. LOCIW = LOCIB LOCDIN = LOCRB LOCR = LOCDIN + N LOCZ = LOCR + N LOCP = LOCZ + N LOCRR = LOCP + N LOCZZ = LOCRR + N LOCPP = LOCZZ + N LOCDZ = LOCPP + N LOCW = LOCDZ + N ! ! Check the workspace allocations. CALL DCHKW( 'DSDBCG', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR ) IF( IERR.NE.0 ) RETURN ! IWORK(4) = LOCDIN IWORK(9) = LOCIW IWORK(10) = LOCW ! CALL DSDS(N, NELT, IA, JA, A, ISYM, RWORK(LOCDIN)) ! ! Perform the Diagonally Scaled BiConjugate gradient algorithm. CALL DBCG(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSMTV, & DSDI, DSDI, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, & RWORK(LOCR), RWORK(LOCZ), RWORK(LOCP), & RWORK(LOCRR), RWORK(LOCZZ), RWORK(LOCPP), & RWORK(LOCDZ), RWORK(1), IWORK(1)) RETURN !------------- LAST LINE OF DSDBCG FOLLOWS ---------------------------- END SUBROUTINE DSLUBC(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) !***BEGIN PROLOGUE DSLUBC !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2A4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(SSLUBC-D), ! Non-Symmetric Linear system, Sparse, ! Iterative incomplete LU Precondition !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Incomplete LU BiConjugate Gradient Sparse Ax=b solver. ! Routine to solve a linear system Ax = b using the ! BiConjugate Gradient method with Incomplete LU ! decomposition preconditioning. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, LENW, IWORK(NEL+NU+4*N+2), LENIW ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(NEL+NU+8*N) ! ! CALL DSLUBC(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW) ! ! *Arguments: ! N :IN Integer. ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :INOUT Integer IA(NELT). ! JA :INOUT Integer JA(NELT). ! A :INOUT Double Precision A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See "Description", ! below. If the SLAP Triad format is chosen it is changed ! internally to the SLAP Column format. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN( ) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IERR = 7 => Incomplete factorization broke down ! and was fudged. Resulting preconditioning may ! be less than the best. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! RWORK :WORK Double Precision RWORK(LENW). ! Double Precision array used for workspace. NEL is the ! number of non- ! zeros in the lower triangle of the matrix (including the ! diagonal). NU is the number of nonzeros in the upper ! triangle of the matrix (including the diagonal). ! LENW :IN Integer. ! Length of the double precision workspace, RWORK. ! LENW >= NEL+NU+8*N. ! IWORK :WORK Integer IWORK(LENIW). ! Integer array used for workspace. NEL is the number of non- ! zeros in the lower triangle of the matrix (including the ! diagonal). NU is the number of nonzeros in the upper ! triangle of the matrix (including the diagonal). ! Upon return the following locations of IWORK hold information ! which may be of use to the user: ! IWORK(9) Amount of Integer workspace actually used. ! IWORK(10) Amount of Double Precision workspace actually used. ! LENIW :IN Integer. ! Length of the integer workspace, IWORK. ! LENIW >= NEL+NU+4*N+12. ! ! *Description: ! This routine is simply a driver for the DBCGN routine. It ! calls the DSILUS routine to set up the preconditioning and ! then calls DBCGN with the appropriate MATVEC, MTTVEC and ! MSOLVE, MTSOLV routines. ! ! The Sparse Linear Algebra Package (SLAP) utilizes two matrix ! data structures: 1) the SLAP Triad format or 2) the SLAP ! Column format. The user can hand this routine either of the ! of these data structures and SLAP will figure out which on ! is being used and act accordingly. ! ! =================== S L A P Triad format =================== ! ! This routine requires that the matrix A be stored in the ! SLAP Triad format. In this format only the non-zeros are ! stored. They may appear in *ANY* order. The user supplies ! three arrays of length NELT, where NELT is the number of ! non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For ! each non-zero the user puts the row and column index of that ! matrix element in the IA and JA arrays. The value of the ! non-zero matrix element is placed in the corresponding ! location of the A array. This is an extremely easy data ! structure to generate. On the other hand it is not too ! efficient on vector computers for the iterative solution of ! linear systems. Hence, SLAP changes this input data ! structure to the SLAP Column format for the iteration (but ! does not change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *Side Effects: ! The SLAP Triad format (IA, JA, A) is modified internally to ! be the SLAP Column format. See above. ! ! *See Also: ! DBCG, SDBCG !***REFERENCES (NONE) !***ROUTINES CALLED DS2Y, DCHKW, DSILUS, DBCG, DSMV, DSMTV !***END PROLOGUE DSLUBC IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER INTEGER IERR, IUNIT, LENW, IWORK(LENIW), LENIW DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW) EXTERNAL DSMV, DSMTV, DSLUI, DSLUTI PARAMETER (LOCRB=1, LOCIB=11) ! ! Change the SLAP input matrix IA, JA, A to SLAP-Column format. !***FIRST EXECUTABLE STATEMENT DSLUBC IERR = 0 IF( N.LT.1 .OR. NELT.LT.1 ) THEN IERR = 3 RETURN end if CALL DS2Y( N, NELT, IA, JA, A, ISYM ) ! ! Count number of Non-Zero elements preconditioner ILU matrix. ! Then set up the work arrays. NL = 0 NU = 0 DO 20 ICOL = 1, N ! Don't count diagonal. JBGN = JA(ICOL)+1 JEND = JA(ICOL+1)-1 IF( JBGN.LE.JEND ) THEN !VD$ NOVECTOR DO J = JBGN, JEND IF( IA(J).GT.ICOL ) THEN NL = NL + 1 IF( ISYM.NE.0 ) NU = NU + 1 ELSE NU = NU + 1 end if end do end if 20 CONTINUE ! LOCIL = LOCIB LOCJL = LOCIL + N+1 LOCIU = LOCJL + NL LOCJU = LOCIU + NU LOCNR = LOCJU + N+1 LOCNC = LOCNR + N LOCIW = LOCNC + N ! LOCL = LOCRB LOCDIN = LOCL + NL LOCU = LOCDIN + N LOCR = LOCU + NU LOCZ = LOCR + N LOCP = LOCZ + N LOCRR = LOCP + N LOCZZ = LOCRR + N LOCPP = LOCZZ + N LOCDZ = LOCPP + N LOCW = LOCDZ + N ! ! Check the workspace allocations. CALL DCHKW( 'DSLUBC', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR ) IF( IERR.NE.0 ) RETURN ! IWORK(1) = LOCIL IWORK(2) = LOCJL IWORK(3) = LOCIU IWORK(4) = LOCJU IWORK(5) = LOCL IWORK(6) = LOCDIN IWORK(7) = LOCU IWORK(9) = LOCIW IWORK(10) = LOCW ! ! Compute the Incomplete LU decomposition. CALL DSILUS( N, NELT, IA, JA, A, ISYM, NL, IWORK(LOCIL), & IWORK(LOCJL), RWORK(LOCL), RWORK(LOCDIN), NU, IWORK(LOCIU), & IWORK(LOCJU), RWORK(LOCU), IWORK(LOCNR), IWORK(LOCNC) ) ! ! Perform the incomplete LU preconditioned ! BiConjugate Gradient algorithm. CALL DBCG(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSMTV, & DSLUI, DSLUTI, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, & RWORK(LOCR), RWORK(LOCZ), RWORK(LOCP), & RWORK(LOCRR), RWORK(LOCZZ), RWORK(LOCPP), & RWORK(LOCDZ), RWORK, IWORK ) RETURN !------------- LAST LINE OF DSLUBC FOLLOWS ---------------------------- END FUNCTION ISDBCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, & TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, PP, DZ, & RWORK, IWORK, AK, BK, BNRM, SOLNRM) !***BEGIN PROLOGUE ISDBCG !***REFER TO DBCG, DSDBCG, DSLUBC !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2A4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(ISDBCG-D), ! Non-Symmetric Linear system, Sparse, ! Iterative Precondition, Stop Test !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Preconditioned BiConjugate Gradient Stop Test. ! This routine calculates the stop test for the BiConjugate ! Gradient iteration scheme. It returns a nonzero if the ! error estimate (the type of which is determined by ITOL) ! is less than the user specified tolerance TOL. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER ! INTEGER IERR, IUNIT, IWORK(USER DEFINED) ! DOUBLE PRECISION B(N), X(N), A(N), TOL, ERR, R(N), Z(N), P(N) ! DOUBLE PRECISION RR(N), ZZ(N), PP(N), DZ(N) ! DOUBLE PRECISION RWORK(USER DEFINED), AK, BK, BNRM, SOLNRM ! EXTERNAL MSOLVE ! ! IF( ISDBCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, RR, ZZ, PP, DZ, ! $ RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) ! $ THEN ITERATION DONE ! ! *Arguments: ! N :IN Integer ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :IN Integer IA(NELT). ! JA :IN Integer JA(NELT). ! A :IN Double Precision A(NELT). ! These arrays contain the matrix data structure for A. ! It could take any form. See "Description", in the SLAP ! routine DBCG for more late breaking details... ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! MSOLVE :EXT External. ! Name of a routine which solves a linear system MZ = R for Z ! given R with the preconditioning matrix M (M is supplied via ! RWORK and IWORK arrays). The name of the MSOLVE routine ! must be declared external in the calling program. The ! calling sequence of MSLOVE is: ! CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! Where N is the number of unknowns, R is the right-hand side ! vector, and Z is the solution upon return. NELT, IA, JA, A ! and ISYM define the SLAP matrix data structure: see ! Description, below. RWORK is a double precision array that ! can be used ! to pass necessary preconditioning information and/or ! workspace to MSOLVE. IWORK is an integer work array for the ! same purpose as RWORK. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than tol, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than tol) through a common block, ! COMMON /SOLBLK/ SOLN( ) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than tol. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Error flag. IERR is set to 3 if ITOL is not on of the ! acceptable values, see above. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! R :IN Double Precision R(N). ! The residual r = b - Ax. ! Z :WORK Double Precision Z(N). ! P :DUMMY Double Precision P(N). ! RR :DUMMY Double Precision RR(N). ! ZZ :DUMMY Double Precision ZZ(N). ! PP :DUMMY Double Precision PP(N). ! DZ :WORK Double Precision DZ(N). ! If ITOL.eq.0 then DZ is used to hold M-inv * B on the first ! call. If ITOL.eq.11 then DZ is used to hold X-SOLN. ! RWORK :WORK Double Precision RWORK(USER DEFINED). ! Double Precision array that can be used for workspace in ! MSOLVE and MTSOLV. ! IWORK :WORK Integer IWORK(USER DEFINED). ! Integer array that can be used for workspace in MSOLVE ! and MTSOLV. ! AK :IN Double Precision. ! Current iterate BiConjugate Gradient iteration parameter. ! BK :IN Double Precision. ! Current iterate BiConjugate Gradient iteration parameter. ! BNRM :INOUT Double Precision. ! Norm of the right hand side. Type of norm depends on ITOL. ! Calculated only on the first call. ! SOLNRM :INOUT Double Precision. ! 2-Norm of the true solution, SOLN. Only computed and used ! if ITOL = 11. ! ! *Function Return Values: ! 0 : Error estimate (determined by ITOL) is *NOT* less than the ! specified tolerance, TOL. The iteration must continue. ! 1 : Error estimate (determined by ITOL) is less than the ! specified tolerance, TOL. The iteration can be considered ! complete. ! ! *Precision: Double Precision !***REFERENCES (NONE) !***ROUTINES CALLED MSOLVE, DNRM2 !***COMMON BLOCKS SOLBLK !***END PROLOGUE ISDBCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX INTEGER ITER, IERR, IUNIT, IWORK(1) DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N), P(N) DOUBLE PRECISION RR(N), ZZ(N), PP(N), DZ(N), RWORK(*) DOUBLE PRECISION AK, BK, BNRM, SOLNRM COMMON /SOLBLK/ SOLN(1) EXTERNAL MSOLVE ! !***FIRST EXECUTABLE STATEMENT ISDBCG ISDBCG = 0 ! IF( ITOL.EQ.1 ) THEN ! err = ||Residual||/||RightHandSide|| (2-Norms). IF(ITER .EQ. 0) BNRM = DNRM2(N, B, 1) ERR = DNRM2(N, R, 1)/BNRM ELSE IF( ITOL.EQ.2 ) THEN ! -1 -1 ! err = ||M Residual||/||M RightHandSide|| (2-Norms). IF(ITER .EQ. 0) THEN CALL MSOLVE(N, B, DZ, NELT, IA, JA, A, ISYM, RWORK, IWORK) BNRM = DNRM2(N, DZ, 1) end if ERR = DNRM2(N, Z, 1)/BNRM ELSE IF( ITOL.EQ.11 ) THEN ! err = ||x-TrueSolution||/||TrueSolution|| (2-Norms). IF(ITER .EQ. 0) SOLNRM = DNRM2(N, SOLN, 1) DZ(1:n) = X(1:n) - SOLN(1:n) ERR = DNRM2(N, DZ, 1)/SOLNRM ELSE ! ! If we get here ITOL is not one of the acceptable values. ERR = 1.0E10 IERR = 3 end if ! IF(IUNIT .NE. 0) THEN IF( ITER.EQ.0 ) THEN WRITE(IUNIT,1000) N, ITOL end if WRITE(IUNIT,1010) ITER, ERR, AK, BK end if IF(ERR .LE. TOL) ISDBCG = 1 ! RETURN 1000 FORMAT(' Preconditioned BiConjugate Gradient for N, ITOL = ', & I5,I5,/' ITER',' Error Estimate',' Alpha', & ' Beta') 1010 FORMAT(1X,I4,1X,E16.7,1X,E16.7,1X,E16.7) !------------- LAST LINE OF ISDBCG FOLLOWS ---------------------------- END SUBROUTINE DCG(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSOLVE, & ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, & RWORK, IWORK ) !***BEGIN PROLOGUE DCG !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2B4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(DCG-D), ! Symmetric Linear system, Sparse, Iterative Precondition !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Preconditioned Conjugate Gradient iterative Ax=b solver. ! Routine to solve a symmetric positive definite linear ! system Ax = b using the Preconditioned Conjugate ! Gradient method. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, IWORK(USER DEFINABLE) ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N) ! DOUBLE PRECISION P(N), DZ(N), RWORK(USER DEFINABLE) ! EXTERNAL MATVEC, MSOLVE ! ! CALL DCG(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MSLOVE, ! $ ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, ! $ RWORK, IWORK ) ! ! *Arguments: ! N :IN Integer. ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :IN Integer IA(NELT). ! JA :IN Integer JA(NELT). ! A :IN Integer A(NELT). ! These arrays contain the matrix data structure for A. ! It could take any form. See ``Description'', below ! for more late breaking details... ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! MATVEC :EXT External. ! Name of a routine which performs the matrix vector multiply ! Y = A*X given A and X. The name of the MATVEC routine must ! be declared external in the calling program. The calling ! sequence to MATVEC is: ! ! CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM ) ! ! Where N is the number of unknowns, Y is the product A*X ! upon return X is an input vector, NELT is the number of ! non-zeros in the SLAP IA, JA, A storage for the matrix A. ! ISYM is a flag which, if non-zero, denotest that A is ! symmetric and only the lower or upper triangle is stored. ! MSOLVE :EXT External. ! Name of a routine which solves a linear system MZ = R for ! Z given R with the preconditioning matrix M (M is supplied via ! RWORK and IWORK arrays). The name of the MSOLVE routine must ! be declared external in the calling program. The calling ! sequence to MSLOVE is: ! ! CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! ! Where N is the number of unknowns, R is the right-hand side ! vector, and Z is the solution upon return. RWORK is a double ! precision ! array that can be used to pass necessary preconditioning ! information and/or workspace to MSOLVE. IWORK is an integer ! work array for the same purpose as RWORK. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! R :WORK Double Precision R(N). ! Z :WORK Double Precision Z(N). ! P :WORK Double Precision P(N). ! DZ :WORK Double Precision DZ(N). ! RWORK :WORK Double Precision RWORK(USER DEFINABLE). ! Double Precision array that can be used by MSOLVE. ! IWORK :WORK Integer IWORK(USER DEFINABLE). ! Integer array that can be used by MSOLVE. ! ! *Description ! This routine does not care what matrix data structure is ! used for A and M. It simply calls the MATVEC and MSOLVE ! routines, with the arguments as described above. The user ! could write any type of structure and the appropriate MATVEC ! and MSOLVE routines. It is assumed that A is stored in the ! IA, JA, A arrays in some fashion and that M (or INV(M)) is ! stored in IWORK and RWORK in some fashion. The SLAP ! routines DSDCG and DSICCG are examples of this procedure. ! ! Two examples of matrix data structures are the: 1) SLAP ! Triad format and 2) SLAP Column format. ! ! =================== S L A P Triad format =================== ! ! In this format only the non-zeros are stored. They may ! appear in *ANY* order. The user supplies three arrays of ! length NELT, where NELT is the number of non-zeros in the ! matrix: (IA(NELT), JA(NELT), A(NELT)). For each non-zero ! the user puts the row and column index of that matrix ! element in the IA and JA arrays. The value of the non-zero ! matrix element is placed in the corresponding location of ! the A array. This is an extremely easy data structure to ! generate. On the other hand it is not too efficient on ! vector computers for the iterative solution of linear ! systems. Hence, SLAP changes this input data structure to ! the SLAP Column format for the iteration (but does not ! change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *See Also: ! DSDCG, DSICCG !***REFERENCES 1. Louis Hageman \& David Young, ``Applied Iterative ! Methods'', Academic Press, New York (1981) ISBN ! 0-12-313340-8. ! ! 2. Concus, Golub \& O'Leary, ``A Generalized Conjugate ! Gradient Method for the Numerical Solution of ! Elliptic Partial Differential Equations,'' in Sparse ! Matrix Computations (Bunch \& Rose, Eds.), Academic ! Press, New York (1979). !***ROUTINES CALLED MATVEC, MSOLVE, ISDCG, DCOPY, DDOT, DAXPY, D1MACH !***END PROLOGUE DCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER INTEGER IUNIT, IERR, IWORK(*) DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N), P(N) DOUBLE PRECISION DZ(N), RWORK(*) EXTERNAL MATVEC, MSOLVE ! ! Check some of the input data. !***FIRST EXECUTABLE STATEMENT DCG ITER = 0 IERR = 0 IF( N.LT.1 ) THEN IERR = 3 RETURN end if TOLMIN = 500.0 * epsilon ( tolmin ) IF( TOL.LT.TOLMIN ) THEN TOL = TOLMIN IERR = 4 end if ! ! Calculate initial residual and pseudo-residual, and check ! stopping criterion. CALL MATVEC(N, X, R, NELT, IA, JA, A, ISYM) R(1:n) = B(1:n) - R(1:n) CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! IF( ISDCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, & RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) GO TO 200 IF( IERR.NE.0 ) RETURN ! ! ***** Iteration loop ***** ! DO 100 K=1,ITMAX ITER = K ! ! Calculate coefficient bk and direction vector p. BKNUM = DDOT(N, Z, 1, R, 1) IF( BKNUM.LE.0.0D0 ) THEN IERR = 5 RETURN end if IF(ITER .EQ. 1) THEN CALL DCOPY(N, Z, 1, P, 1) ELSE BK = BKNUM/BKDEN DO 20 I = 1, N P(I) = Z(I) + BK*P(I) 20 CONTINUE end if BKDEN = BKNUM ! ! Calculate coefficient ak, new iterate x, new residual r, ! and new pseudo-residual z. CALL MATVEC(N, P, Z, NELT, IA, JA, A, ISYM) AKDEN = DDOT(N, P, 1, Z, 1) IF( AKDEN.LE.0.0D0 ) THEN IERR = 6 RETURN end if AK = BKNUM/AKDEN CALL DAXPY(N, AK, P, 1, X, 1) CALL DAXPY(N, -AK, Z, 1, R, 1) CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! ! check stopping criterion. IF( ISDCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, RWORK, & IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) GO TO 200 ! 100 CONTINUE ! ! ***** end of loop ***** ! ! stopping criterion not satisfied. ITER = ITMAX + 1 IERR = 2 ! 200 RETURN END SUBROUTINE DSDCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) !***BEGIN PROLOGUE DSDCG !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2B4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(DSDCG-D), ! Symmetric Linear system, Sparse, Iterative Precondition !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Diagonally Scaled Conjugate Gradient Sparse Ax=b Solver. ! Routine to solve a symmetric positive definite linear ! system Ax = b using the Preconditioned Conjugate ! Gradient method. The preconditioner is diagonal ! scaling. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, LENW, IWORK(10), LENIW ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(5*N) ! ! CALL DSDCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) ! ! *Arguments: ! N :IN Integer. ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :INOUT Integer IA(NELT). ! JA :INOUT Integer JA(NELT). ! A :INOUT Integer A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See ``Description'', ! below. If the SLAP Triad format is chosen it is changed ! internally to the SLAP Column format. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! RWORK :WORK Double Precision RWORK(LENW). ! Double Precision array used for workspace. ! LENW :IN Integer. ! Length of the double precision workspace, RWORK. LENW >= 5*N. ! IWORK :WORK Integer IWORK(LENIW). ! Used to hold pointers into the double precision workspace, ! RWORK. Upon return the following locations of IWORK hold ! information which may be of use to the user: ! IWORK(9) Amount of Integer workspace actually used. ! IWORK(10) Amount of Double Precision workspace actually used. ! LENIW :IN Integer. ! Length of the integer workspace, IWORK. LENIW >= 10. ! ! *Description: ! This routine performs preconditioned conjugate gradient ! method on the symmetric positive definite linear system ! Ax=b. The preconditioner is M = DIAG(A), the diagonal of ! the matrix A. This is the simplest of preconditioners and ! vectorizes very well. This routine is simply a driver for ! the DCG routine. It calls the DSDS routine to set up the ! preconditioning and then calls DCG with the appropriate ! MATVEC and MSOLVE routines. ! ! The Sparse Linear Algebra Package (SLAP) utilizes two matrix ! data structures: 1) the SLAP Triad format or 2) the SLAP ! Column format. The user can hand this routine either of the ! of these data structures and SLAP will figure out which on ! is being used and act accordingly. ! ! =================== S L A P Triad format =================== ! ! This routine requires that the matrix A be stored in the ! SLAP Triad format. In this format only the non-zeros are ! stored. They may appear in *ANY* order. The user supplies ! three arrays of length NELT, where NELT is the number of ! non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For ! each non-zero the user puts the row and column index of that ! matrix element in the IA and JA arrays. The value of the ! non-zero matrix element is placed in the corresponding ! location of the A array. This is an extremely easy data ! structure to generate. On the other hand it is not too ! efficient on vector computers for the iterative solution of ! linear systems. Hence, SLAP changes this input data ! structure to the SLAP Column format for the iteration (but ! does not change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *Side Effects: ! The SLAP Triad format (IA, JA, A) is modified internally to ! be the SLAP Column format. See above. ! ! *See Also: ! DCG, DSICCG !***REFERENCES 1. Louis Hageman \& David Young, ``Applied Iterative ! Methods'', Academic Press, New York (1981) ISBN ! 0-12-313340-8. ! 2. Concus, Golub \& O'Leary, ``A Generalized Conjugate ! Gradient Method for the Numerical Solution of ! Elliptic Partial Differential Equations,'' in Sparse ! Matrix Computations (Bunch \& Rose, Eds.), Academic ! Press, New York (1979). !***ROUTINES CALLED DS2Y, DCHKW, DSDS, DCG !***END PROLOGUE DSDCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL INTEGER ITMAX, ITER, IERR, IUNIT, LENW, IWORK(LENIW), LENIW DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW) EXTERNAL DSMV, DSDI PARAMETER (LOCRB=1, LOCIB=11) ! ! Modify the SLAP matrix data structure to YSMP-Column. !***FIRST EXECUTABLE STATEMENT DSDCG IERR = 0 IF( N.LT.1 .OR. NELT.LT.1 ) THEN IERR = 3 RETURN end if CALL DS2Y( N, NELT, IA, JA, A, ISYM ) ! ! Set up the work arrays. ! Compute the inverse of the diagonal of the matrix. This ! will be used as the preconditioner. LOCIW = LOCIB ! LOCD = LOCRB LOCR = LOCD + N LOCZ = LOCR + N LOCP = LOCZ + N LOCDZ = LOCP + N LOCW = LOCDZ + N ! ! Check the workspace allocations. CALL DCHKW( 'DSDCG', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR ) IF( IERR.NE.0 ) RETURN ! IWORK(4) = LOCD IWORK(9) = LOCIW IWORK(10) = LOCW ! CALL DSDS(N, NELT, IA, JA, A, ISYM, RWORK(LOCD)) ! ! Do the Preconditioned Conjugate Gradient. CALL DCG(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSDI, & ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR), & RWORK(LOCZ), RWORK(LOCP), RWORK(LOCDZ), RWORK, IWORK) RETURN !------------- LAST LINE OF DSDCG FOLLOWS ----------------------------- END SUBROUTINE DSICCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) !***BEGIN PROLOGUE DSICCG !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2B4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(DSICCG-D), ! Symmetric Linear system, Sparse, ! Iterative Precondition, Incomplete Cholesky !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Incomplete Cholesky Conjugate Gradient Sparse Ax=b Solver. ! Routine to solve a symmetric positive definite linear ! system Ax = b using the incomplete Cholesky ! Preconditioned Conjugate Gradient method. ! !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, LENW, IWORK(NEL+2*n+1), LENIW ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(NEL+5*N) ! ! CALL DSICCG(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) ! ! *Arguments: ! N :IN Integer. ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :INOUT Integer IA(NELT). ! JA :INOUT Integer JA(NELT). ! A :INOUT Integer A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See ``Description'', ! below. If the SLAP Triad format is chosen it is changed ! internally to the SLAP Column format. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IERR = 7 => Incomplete factorization broke down ! and was fudged. Resulting preconditioning may ! be less than the best. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! RWORK :WORK Double Precision RWORK(LENW). ! Double Precision array used for workspace. NEL is the ! number of non- ! zeros in the lower triangle of the matrix (including the ! diagonal) ! LENW :IN Integer. ! Length of the double precision workspace, RWORK. ! LENW >= NEL+5*N. ! IWORK :WORK Integer IWORK(LENIW). ! Integer array used for workspace. NEL is the number of non- ! zeros in the lower triangle of the matrix (including the ! diagonal). ! Upon return the following locations of IWORK hold information ! which may be of use to the user: ! IWORK(9) Amount of Integer workspace actually used. ! IWORK(10) Amount of Double Precision workspace actually used. ! LENIW :IN Integer. ! Length of the integer workspace, IWORK. LENIW >= NEL+N+11. ! ! *Description: ! This routine performs preconditioned conjugate gradient ! method on the symmetric positive definite linear system ! Ax=b. The preconditioner is the incomplete Cholesky (IC) ! factorization of the matrix A. See DSICS for details about ! the incomplete factorization algorithm. One should note ! here however, that the IC factorization is a slow process ! and that one should save factorizations for reuse, if ! possible. The MSOLVE operation (handled in DSLLTI) does ! vectorize on machines with hardware gather/scatter and is ! quite fast. ! ! The Sparse Linear Algebra Package (SLAP) utilizes two matrix ! data structures: 1) the SLAP Triad format or 2) the SLAP ! Column format. The user can hand this routine either of the ! of these data structures and SLAP will figure out which on ! is being used and act accordingly. ! ! =================== S L A P Triad format =================== ! ! This routine requires that the matrix A be stored in the ! SLAP Triad format. In this format only the non-zeros are ! stored. They may appear in *ANY* order. The user supplies ! three arrays of length NELT, where NELT is the number of ! non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For ! each non-zero the user puts the row and column index of that ! matrix element in the IA and JA arrays. The value of the ! non-zero matrix element is placed in the corresponding ! location of the A array. This is an extremely easy data ! structure to generate. On the other hand it is not too ! efficient on vector computers for the iterative solution of ! linear systems. Hence, SLAP changes this input data ! structure to the SLAP Column format for the iteration (but ! does not change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *Side Effects: ! The SLAP Triad format (IA, JA, A) is modified internally to be ! the SLAP Column format. See above. ! ! *See Also: ! DCG, DSLLTI !***REFERENCES 1. Louis Hageman \& David Young, ``Applied Iterative ! Methods'', Academic Press, New York (1981) ISBN ! 0-12-313340-8. ! 2. Concus, Golub \& O'Leary, ``A Generalized Conjugate ! Gradient Method for the Numerical Solution of ! Elliptic Partial Differential Equations,'' in Sparse ! Matrix Computations (Bunch \& Rose, Eds.), Academic ! Press, New York (1979). !***ROUTINES CALLED DS2Y, DCHKW, DSICS, XERRWV, DCG !***END PROLOGUE DSICCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL INTEGER ITMAX, ITER, IUNIT, LENW, IWORK(LENIW), LENIW DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW) EXTERNAL DSMV, DSLLTI PARAMETER (LOCRB=1, LOCIB=11) ! ! Change the SLAP input matrix IA, JA, A to SLAP-Column format. !***FIRST EXECUTABLE STATEMENT DSICCG IERR = 0 IF( N.LT.1 .OR. NELT.LT.1 ) THEN IERR = 3 RETURN end if CALL DS2Y( N, NELT, IA, JA, A, ISYM ) ! ! Count number of elements in lower triangle of the matrix. ! Then set up the work arrays. IF( ISYM.EQ.0 ) THEN NEL = (NELT + N)/2 ELSE NEL = NELT end if ! LOCJEL = LOCIB LOCIEL = LOCJEL + NEL LOCIW = LOCIEL + N + 1 ! LOCEL = LOCRB LOCDIN = LOCEL + NEL LOCR = LOCDIN + N LOCZ = LOCR + N LOCP = LOCZ + N LOCDZ = LOCP + N LOCW = LOCDZ + N ! ! Check the workspace allocations. CALL DCHKW( 'DSICCG', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR ) IF( IERR.NE.0 ) RETURN ! IWORK(1) = NEL IWORK(2) = LOCJEL IWORK(3) = LOCIEL IWORK(4) = LOCEL IWORK(5) = LOCDIN IWORK(9) = LOCIW IWORK(10) = LOCW ! ! Compute the Incomplete Cholesky decomposition. ! CALL DSICS(N, NELT, IA, JA, A, ISYM, NEL, IWORK(LOCIEL), & IWORK(LOCJEL), RWORK(LOCEL), RWORK(LOCDIN), & RWORK(LOCR), IERR ) IF( IERR.NE.0 ) THEN CALL XERRWV('DSICCG: Warning...IC factorization broke down '// & 'on step i1. Diagonal was set to unity and '// & 'factorization proceeded.', 113, 1, 1, 1, IERR, 0, & 0, 0.0, 0.0 ) IERR = 7 end if ! ! Do the Preconditioned Conjugate Gradient. CALL DCG(N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSLLTI, & ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR), & RWORK(LOCZ), RWORK(LOCP), RWORK(LOCDZ), RWORK(1), & IWORK(1)) RETURN !------------- LAST LINE OF DSICCG FOLLOWS ---------------------------- END FUNCTION ISDCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, & TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, & RWORK, IWORK, AK, BK, BNRM, SOLNRM) !***BEGIN PROLOGUE ISDCG !***REFER TO DCG, DSDCG, DSICCG !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2B4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(ISDCG-D), ! Linear system, Sparse, Stop Test !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Preconditioned Conjugate Gradient Stop Test. ! This routine calculates the stop test for the Conjugate ! Gradient iteration scheme. It returns a nonzero if the ! error estimate (the type of which is determined by ITOL) ! is less than the user specified tolerance TOL. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER ! INTEGER IERR, IUNIT, IWORK(USER DEFINED) ! DOUBLE PRECISION B(N), X(N), A(N), TOL, ERR, R(N), Z(N) ! DOUBLE PRECISION P(N), DZ(N), RWORK(USER DEFINED), AK, BK ! DOUBLE PRECISION BNRM, SOLNRM ! EXTERNAL MSOLVE ! ! IF( ISDCG(N, B, X, NELT, IA, JA, A, ISYM, MSOLVE, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, DZ, RWORK, IWORK, ! $ AK, BK, BNRM, SOLNRM) .NE. 0 ) THEN ITERATION DONE ! ! *Arguments: ! N :IN Integer. ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :IN Double Precision X(N). ! The current approximate solution vector. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :IN Integer IA(NELT). ! JA :IN Integer JA(NELT). ! A :IN Double Precision A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See ``Description'' ! in the DCG, DSDCG or DSICCG routines. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! MSOLVE :EXT External. ! Name of a routine which solves a linear system MZ = R for ! Z given R with the preconditioning matrix M (M is supplied via ! RWORK and IWORK arrays). The name of the MSOLVE routine must ! be declared external in the calling program. The calling ! sequence to MSLOVE is: ! CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! Where N is the number of unknowns, R is the right-hand side ! vector, and Z is the solution upon return. RWORK is a double ! precision ! array that can be used to pass necessary preconditioning ! information and/or workspace to MSOLVE. IWORK is an integer ! work array for the same purpose as RWORK. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than tol, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the ``exact'' ! solution or a very accurate approximation (one with an error ! much less than tol) through a common block, ! COMMON /SOLBLK/ SOLN( ) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than tol. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :IN Integer. ! The iteration for which to check for convergence. ! ERR :OUT Double Precision. ! Error estimate of error in the X(N) approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Error flag. IERR is set to 3 if ITOL is not on of the ! acceptable values, see above. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! R :IN Double Precision R(N). ! The residual R = B-AX. ! Z :WORK Double Precision Z(N). ! Workspace used to hold the pseudo-residual M Z = R. ! P :IN Double Precision P(N). ! The conjugate direction vector. ! DZ :WORK Double Precision DZ(N). ! Workspace used to hold temporary vector(s). ! RWORK :WORK Double Precision RWORK(USER DEFINABLE). ! Double Precision array that can be used by MSOLVE. ! IWORK :WORK Integer IWORK(USER DEFINABLE). ! Integer array that can be used by MSOLVE. ! BNRM :INOUT Double Precision. ! Norm of the right hand side. Type of norm depends on ITOL. ! Calculated only on the first call. ! SOLNRM :INOUT Double Precision. ! 2-Norm of the true solution, SOLN. Only computed and used ! if ITOL = 11. ! ! *Function Return Values: ! 0 : Error estimate (determined by ITOL) is *NOT* less than the ! specified tolerance, TOL. The iteration must continue. ! 1 : Error estimate (determined by ITOL) is less than the ! specified tolerance, TOL. The iteration can be considered ! complete. ! ! *Precision: Double Precision ! *See Also: ! DCG, DSDCG, DSICCG ! ! *Cautions: ! This routine will attempt to write to the fortran logical output ! unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that ! this logical unit must be attached to a file or terminal ! before calling this routine with a non-zero value for IUNIT. ! This routine does not check for the validity of a non-zero IUNIT ! unit number. !***REFERENCES (NONE) !***ROUTINES CALLED MSOLVE, DNRM2 !***COMMON BLOCKS SOLBLK !***END PROLOGUE ISDCG IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX INTEGER ITER, IERR, IUNIT, IWORK(*) DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N) DOUBLE PRECISION Z(N), P(N), DZ(N), RWORK(*) EXTERNAL MSOLVE COMMON /SOLBLK/ SOLN(1) ! !***FIRST EXECUTABLE STATEMENT ISDCG ISDCG = 0 ! IF( ITOL.EQ.1 ) THEN ! err = ||Residual||/||RightHandSide|| (2-Norms). IF(ITER .EQ. 0) BNRM = DNRM2(N, B, 1) ERR = DNRM2(N, R, 1)/BNRM ELSE IF( ITOL.EQ.2 ) THEN ! -1 -1 ! err = ||M Residual||/||M RightHandSide|| (2-Norms). IF(ITER .EQ. 0) THEN CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) BNRM = DNRM2(N, DZ, 1) end if ERR = DNRM2(N, Z, 1)/BNRM ELSE IF( ITOL.EQ.11 ) THEN ! err = ||x-TrueSolution||/||TrueSolution|| (2-Norms). IF(ITER .EQ. 0) SOLNRM = DNRM2(N, SOLN, 1) DZ(1:n) = X(1:n) - SOLN(1:n) ERR = DNRM2(N, DZ, 1)/SOLNRM ELSE ! ! If we get here ITOL is not one of the acceptable values. ERR = 1.0E10 IERR = 3 end if ! IF(IUNIT .NE. 0) THEN IF( ITER.EQ.0 ) THEN WRITE(IUNIT,1000) N, ITOL end if WRITE(IUNIT,1010) ITER, ERR, AK, BK end if IF(ERR .LE. TOL) ISDCG = 1 RETURN 1000 FORMAT(' Preconditioned Conjugate Gradient for ', & 'N, ITOL = ',I5, I5, & /' ITER',' Error Estimate',' Alpha', & ' Beta') 1010 FORMAT(1X,I4,1X,E16.7,1X,E16.7,1X,E16.7) !------------- LAST LINE OF ISDCG FOLLOWS ------------------------------ END SUBROUTINE DCGN(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MTTVEC, & MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, & ATP, ATZ, DZ, ATDZ, RWORK, IWORK) !***BEGIN PROLOGUE DCGN !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2A4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(DCGN-D), ! Non-Symmetric Linear system solve, Sparse, ! Iterative Precondition, Normal Equations. !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Preconditioned CG Sparse Ax=b Solver for Normal Equations. ! Routine to solve a general linear system Ax = b using the ! Preconditioned Conjugate Gradient method applied to the ! normal equations AA'y = b, x=A'y. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, IWORK(USER DEFINABLE) ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), Z(N) ! DOUBLE PRECISION P(N), ATP(N), ATZ(N), DZ(N), ATDZ(N) ! DOUBLE PRECISION RWORK(USER DEFINABLE) ! EXTERNAL MATVEC, MTTVEC, MSOLVE ! ! CALL DCGN(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MTTVEC, ! $ MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, ! $ Z, P, ATP, ATZ, DZ, ATDZ, RWORK, IWORK) ! ! *Arguments: ! N :IN Integer ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :IN Integer IA(NELT). ! JA :IN Integer JA(NELT). ! A :IN Double Precision A(NELT). ! These arrays contain the matrix data structure for A. ! It could take any form. See "Description", below ! for more late breaking details... ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! MATVEC :EXT External. ! Name of a routine which performs the matrix vector multiply ! y = A*X given A and X. The name of the MATVEC routine must ! be declared external in the calling program. The calling ! sequence to MATVEC is: ! CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM ) ! Where N is the number of unknowns, Y is the product A*X ! upon return X is an input vector, NELT is the number of ! non-zeros in the SLAP-Column IA, JA, A storage for the matrix ! A. ISYM is a flag which, if non-zero, denotes that A is ! symmetric and only the lower or upper triangle is stored. ! MTTVEC :EXT External. ! Name of a routine which performs the matrix transpose vector ! multiply y = A'*X given A and X (where ' denotes transpose). ! The name of the MTTVEC routine must be declared external in ! the calling program. The calling sequence to MTTVEC is the ! same as that for MATVEC, viz.: ! CALL MTTVEC( N, X, Y, NELT, IA, JA, A, ISYM ) ! Where N is the number of unknowns, Y is the product A'*X ! upon return X is an input vector, NELT is the number of ! non-zeros in the SLAP-Column IA, JA, A storage for the matrix ! A. ISYM is a flag which, if non-zero, denotes that A is ! symmetric and only the lower or upper triangle is stored. ! MSOLVE :EXT External. ! Name of a routine which solves a linear system MZ = R for ! Z given R with the preconditioning matrix M (M is supplied via ! RWORK and IWORK arrays). The name of the MSOLVE routine must ! be declared external in the calling program. The calling ! sequence to MSOLVE is: ! CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) ! Where N is the number of unknowns, R is the right-hand side ! vector, and Z is the solution upon return. RWORK is a ! double precision ! array that can be used to pass necessary preconditioning ! information and/or workspace to MSOLVE. IWORK is an integer ! work array for the same purpose as RWORK. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! R :WORK Double Precision R(N). ! Z :WORK Double Precision Z(N). ! P :WORK Double Precision P(N). ! ATP :WORK Double Precision ATP(N). ! ATZ :WORK Double Precision ATZ(N). ! DZ :WORK Double Precision DZ(N). ! ATDZ :WORK Double Precision ATDZ(N). ! RWORK :WORK Double Precision RWORK(USER DEFINABLE). ! Double Precision array that can be used by MSOLVE. ! IWORK :WORK Integer IWORK(USER DEFINABLE). ! Integer array that can be used by MSOLVE. ! ! *Description: ! This routine applies the preconditioned conjugate gradient ! (PCG) method to a non-symmetric system of equations Ax=b. To ! do this the normal equations are solved: ! AA' y = b, where x = A'y. ! In PCG method the iteration count is determined by condition ! -1 ! number of the matrix (M A). In the situation where the ! normal equations are used to solve a non-symmetric system ! the condition number depends on AA' and should therefore be ! much worse than that of A. This is the conventional wisdom. ! When one has a good preconditioner for AA' this may not hold. ! The latter is the situation when DCGN should be tried. ! ! If one is trying to solve a symmetric system, SCG should be ! used instead. ! ! This routine does not care what matrix data structure is ! used for A and M. It simply calls the MATVEC and MSOLVE ! routines, with the arguments as described above. The user ! could write any type of structure and the appropriate MATVEC ! and MSOLVE routines. It is assumed that A is stored in the ! IA, JA, A arrays in some fashion and that M (or INV(M)) is ! stored in IWORK and RWORK) in some fashion. The SLAP ! routines SSDCGN and SSLUCN are examples of this procedure. ! ! Two examples of matrix data structures are the: 1) SLAP ! Triad format and 2) SLAP Column format. ! ! =================== S L A P Triad format =================== ! ! In this format only the non-zeros are stored. They may ! appear in *ANY* order. The user supplies three arrays of ! length NELT, where NELT is the number of non-zeros in the ! matrix: (IA(NELT), JA(NELT), A(NELT)). For each non-zero ! the user puts the row and column index of that matrix ! element in the IA and JA arrays. The value of the non-zero ! matrix element is placed in the corresponding location of ! the A array. This is an extremely easy data structure to ! generate. On the other hand it is not too efficient on ! vector computers for the iterative solution of linear ! systems. Hence, SLAP changes this input data structure to ! the SLAP Column format for the iteration (but does not ! change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *See Also: ! DSDCGN, DSLUCN, ISDCGN !***REFERENCES (NONE) !***ROUTINES CALLED MATVEC, MTTVEC, MSOLVE, ISDCGN, ! DCOPY, DDOT, DAXPY, D1MACH !***END PROLOGUE DCGN IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX, ITER INTEGER IUNIT, IWORK(*) DOUBLE PRECISION B(N), X(N), A(N), R(N), Z(N), P(N) DOUBLE PRECISION ATP(N), ATZ(N), DZ(N), ATDZ(N), RWORK(*) EXTERNAL MATVEC, MTTVEC, MSOLVE ! ! Check user input. !***FIRST EXECUTABLE STATEMENT DCGN ITER = 0 IERR = 0 IF( N.LT.1 ) THEN IERR = 3 RETURN end if TOLMIN = 500.0 * epsilon ( tolmin ) IF( TOL.LT.TOLMIN ) THEN TOL = TOLMIN IERR = 4 end if ! Calculate initial residual and pseudo-residual, and check ! stopping criterion. CALL MATVEC(N, X, R, NELT, IA, JA, A, ISYM) R(1:n) = B(1:n) - R(1:n) CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) CALL MTTVEC(N, Z, ATZ, NELT, IA, JA, A, ISYM) ! IF( ISDCGN(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MTTVEC, MSOLVE, & ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, Z, P, ATP, ATZ, & DZ, ATDZ, RWORK, IWORK, AK, BK, BNRM, SOLNRM) .NE. 0 ) & GO TO 200 IF( IERR.NE.0 ) RETURN ! ! ***** iteration loop ***** ! DO 100 K=1,ITMAX ITER = K ! ! Calculate coefficient BK and direction vector P. BKNUM = DDOT(N, Z, 1, R, 1) IF( BKNUM.LE.0.0D0 ) THEN IERR = 6 RETURN end if IF(ITER .EQ. 1) THEN CALL DCOPY(N, Z, 1, P, 1) ELSE BK = BKNUM/BKDEN DO 20 I = 1, N P(I) = Z(I) + BK*P(I) 20 CONTINUE end if BKDEN = BKNUM ! ! Calculate coefficient AK, new iterate X, new residual R, ! and new pseudo-residual ATZ. IF(ITER .NE. 1) CALL DAXPY(N, BK, ATP, 1, ATZ, 1) CALL DCOPY(N, ATZ, 1, ATP, 1) AKDEN = DDOT(N, ATP, 1, ATP, 1) IF( AKDEN.LE.0.0D0 ) THEN IERR = 6 RETURN end if AK = BKNUM/AKDEN CALL DAXPY(N, AK, ATP, 1, X, 1) CALL MATVEC(N, ATP, Z, NELT, IA, JA, A, ISYM) CALL DAXPY(N, -AK, Z, 1, R, 1) CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) CALL MTTVEC(N, Z, ATZ, NELT, IA, JA, A, ISYM) ! ! check stopping criterion. IF( ISDCGN(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, MTTVEC, & MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, R, & Z, P, ATP, ATZ, DZ, ATDZ, RWORK, IWORK, AK, BK, BNRM, & SOLNRM) .NE. 0) GOTO 200 ! 100 CONTINUE ! ! ***** end of loop ***** ! ! stopping criterion not satisfied. ITER = ITMAX + 1 ! 200 RETURN !------------- LAST LINE OF DCGN FOLLOWS ---------------------------- END SUBROUTINE DSDCGN(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) !***BEGIN PROLOGUE DSDCGN !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2A4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(SSDCGN-D), ! Non-Symmetric Linear system solve, Sparse, ! Iterative Precondition !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Diagonally Scaled CG Sparse Ax=b Solver for Normal Eqn's. ! Routine to solve a general linear system Ax = b using ! diagonal scaling with the Conjugate Gradient method ! applied to the the normal equations, viz., AA'y = b, ! where x = A'y. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, LENW, IWORK, LENIW ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(8*N) ! ! CALL DSDCGN(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW) ! ! *Arguments: ! N :IN Integer ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :INOUT Integer IA(NELT). ! JA :INOUT Integer JA(NELT). ! A :INOUT Double Precision A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See "Description", ! below. If the SLAP Triad format is chosen it is changed ! internally to the SLAP Column format. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! RWORK :WORK Double Precision RWORK(LENW). ! Double Precision array used for workspace. ! LENW :IN Integer. ! Length of the double precision workspace, RWORK. ! LENW >= 8*N. ! IWORK :WORK Integer IWORK(LENIW). ! Used to hold pointers into the RWORK array. ! Upon return the following locations of IWORK hold information ! which may be of use to the user: ! IWORK(9) Amount of Integer workspace actually used. ! IWORK(10) Amount of Double Precision workspace actually used. ! LENIW :IN Integer. ! Length of the integer workspace, IWORK. LENIW >= 10. ! ! *Description: ! This routine is simply a driver for the DCGN routine. It ! calls the DSD2S routine to set up the preconditioning and ! then calls DCGN with the appropriate MATVEC and MSOLVE ! routines. ! ! The Sparse Linear Algebra Package (SLAP) utilizes two matrix ! data structures: 1) the SLAP Triad format or 2) the SLAP ! Column format. The user can hand this routine either of the ! of these data structures and SLAP will figure out which on ! is being used and act accordingly. ! ! =================== S L A P Triad format =================== ! ! This routine requires that the matrix A be stored in the ! SLAP Triad format. In this format only the non-zeros are ! stored. They may appear in *ANY* order. The user supplies ! three arrays of length NELT, where NELT is the number of ! non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For ! each non-zero the user puts the row and column index of that ! matrix element in the IA and JA arrays. The value of the ! non-zero matrix element is placed in the corresponding ! location of the A array. This is an extremely easy data ! structure to generate. On the other hand it is not too ! efficient on vector computers for the iterative solution of ! linear systems. Hence, SLAP changes this input data ! structure to the SLAP Column format for the iteration (but ! does not change it back). ! ! Here is an example of the SLAP Triad storage format for a ! 5x5 Matrix. Recall that the entries may appear in any order. ! ! 5x5 Matrix SLAP Triad format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21 ! |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2 ! | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! =================== S L A P Column format ================== ! This routine requires that the matrix A be stored in the ! SLAP Column format. In this format the non-zeros are stored ! counting down columns (except for the diagonal entry, which ! must appear first in each "column") and are stored in the ! double precision array A. In other words, for each column ! in the matrix put the diagonal entry in A. Then put in the ! other non-zero elements going down the column (except the ! diagonal) in order. The IA array holds the row index for ! each non-zero. The JA array holds the offsets into the IA, ! A arrays for the beginning of each column. That is, ! IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the ! ICOL-th column in IA and A. IA(JA(ICOL+1)-1), ! A(JA(ICOL+1)-1) points to the end of the ICOL-th column. ! Note that we always have JA(N+1) = NELT+1, where N is the ! number of columns in the matrix and NELT is the number of ! non-zeros in the matrix. ! ! Here is an example of the SLAP Column storage format for a ! 5x5 Matrix (in the A and IA arrays '|' denotes the end of a ! column): ! ! 5x5 Matrix SLAP Column format for 5x5 matrix on left. ! 1 2 3 4 5 6 7 8 9 10 11 ! |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35 ! |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3 ! | 0 0 33 0 35| JA: 1 4 6 8 9 12 ! | 0 0 0 44 0| ! |51 0 53 0 55| ! ! *Precision: Double Precision ! *Side Effects: ! The SLAP Triad format (IA, JA, A) is modified internally to be ! the SLAP Column format. See above. ! ! *See Also: ! DCGN, DSD2S, DSMV, DSMTV, DSDI !***REFERENCES (NONE) !***ROUTINES CALLED DS2Y, DCHKW, DSD2S, DCGN, DSMV, DSMTV, DSDI !***END PROLOGUE DSDCGN IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL INTEGER ITMAX, ITER, IERR, IUNIT, LENW, IWORK(LENIW), LENIW DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(LENW) EXTERNAL DSMV, DSMTV, DSDI PARAMETER (LOCRB=1, LOCIB=11) ! ! Modify the SLAP matrix data structure to YSMP-Column. !***FIRST EXECUTABLE STATEMENT DSDCGN IERR = 0 IF( N.LT.1 .OR. NELT.LT.1 ) THEN IERR = 3 RETURN end if CALL DS2Y( N, NELT, IA, JA, A, ISYM ) ! ! Set up the work arrays. ! Compute the inverse of the diagonal of AA'. This will be ! used as the preconditioner. LOCIW = LOCIB ! LOCD = LOCRB LOCR = LOCD + N LOCZ = LOCR + N LOCP = LOCZ + N LOCATP = LOCP + N LOCATZ = LOCATP + N LOCDZ = LOCATZ + N LOCATD = LOCDZ + N LOCW = LOCATD + N ! ! Check the workspace allocations. CALL DCHKW( 'DSDCGN', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR ) IF( IERR.NE.0 ) RETURN ! IWORK(4) = LOCD IWORK(9) = LOCIW IWORK(10) = LOCW ! CALL DSD2S(N, NELT, IA, JA, A, ISYM, RWORK(1)) ! ! Perform Conjugate Gradient algorithm on the normal equations. CALL DCGN( N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSMTV, DSDI, & ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR), & RWORK(LOCZ), RWORK(LOCP), RWORK(LOCATP), RWORK(LOCATZ), & RWORK(LOCDZ), RWORK(LOCATD), RWORK, IWORK ) ! IF( ITER.GT.ITMAX ) IERR = 2 RETURN !------------- LAST LINE OF DSDCGN FOLLOWS ---------------------------- END SUBROUTINE DSLUCN(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, & ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) !***BEGIN PROLOGUE DSLUCN !***DATE WRITTEN 890404 (YYMMDD) !***REVISION DATE 890404 (YYMMDD) !***CATEGORY NO. D2B4 !***KEYWORDS LIBRARY=SLATEC(SLAP), ! TYPE=DOUBLE PRECISION(SSLUCN-D), ! Non-Symmetric Linear system, Sparse, ! Iterative Incomplete LU Precondition !***AUTHOR Greenbaum, Anne, Courant Institute ! Seager, Mark K., (LLNL) ! Lawrence Livermore National Laboratory ! PO BOX 808, L-300 ! Livermore, CA 94550 (415) 423-3141 ! seager@lll-crg.llnl.gov !***PURPOSE Incomplete LU CG Sparse Ax=b Solver for Normal Equations. ! Routine to solve a general linear system Ax = b using the ! incomplete LU decomposition with the Conjugate Gradient ! method applied to the normal equations, viz., AA'y = b, ! x=A'y. !***DESCRIPTION ! *Usage: ! INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX ! INTEGER ITER, IERR, IUNIT, LENW, IWORK(NEL+NU+4*N+2), LENIW ! DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(NEL+NU+8*N) ! ! CALL DSLUCN(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL, ! $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW ) ! ! *Arguments: ! N :IN Integer ! Order of the Matrix. ! B :IN Double Precision B(N). ! Right-hand side vector. ! X :INOUT Double Precision X(N). ! On input X is your initial guess for solution vector. ! On output X is the final approximate solution. ! NELT :IN Integer. ! Number of Non-Zeros stored in A. ! IA :INOUT Integer IA(NELT). ! JA :INOUT Integer JA(NELT). ! A :INOUT Double Precision A(NELT). ! These arrays should hold the matrix A in either the SLAP ! Triad format or the SLAP Column format. See "Description", ! below. If the SLAP Triad format is chosen it is changed ! internally to the SLAP Column format. ! ISYM :IN Integer. ! Flag to indicate symmetric storage format. ! If ISYM=0, all nonzero entries of the matrix are stored. ! If ISYM=1, the matrix is symmetric, and only the upper ! or lower triangle of the matrix is stored. ! ITOL :IN Integer. ! Flag to indicate type of convergence criterion. ! If ITOL=1, iteration stops when the 2-norm of the residual ! divided by the 2-norm of the right-hand side is less than TOL. ! If ITOL=2, iteration stops when the 2-norm of M-inv times the ! residual divided by the 2-norm of M-inv times the right hand ! side is less than TOL, where M-inv is the inverse of the ! diagonal of A. ! ITOL=11 is often useful for checking and comparing different ! routines. For this case, the user must supply the "exact" ! solution or a very accurate approximation (one with an error ! much less than TOL) through a common block, ! COMMON /SOLBLK/ SOLN(1) ! if ITOL=11, iteration stops when the 2-norm of the difference ! between the iterative approximation and the user-supplied ! solution divided by the 2-norm of the user-supplied solution ! is less than TOL. Note that this requires the user to set up ! the "COMMON /SOLBLK/ SOLN(LENGTH)" in the calling routine. ! The routine with this declaration should be loaded before the ! stop test so that the correct length is used by the loader. ! This procedure is not standard Fortran and may not work ! correctly on your system (although it has worked on every ! system the authors have tried). If ITOL is not 11 then this ! common block is indeed standard Fortran. ! TOL :IN Double Precision. ! Convergence criterion, as described above. ! ITMAX :IN Integer. ! Maximum number of iterations. ! ITER :OUT Integer. ! Number of iterations required to reach convergence, or ! ITMAX+1 if convergence criterion could not be achieved in ! ITMAX iterations. ! ERR :OUT Double Precision. ! Error estimate of error in final approximate solution, as ! defined by ITOL. ! IERR :OUT Integer. ! Return error flag. ! IERR = 0 => All went well. ! IERR = 1 => Insufficient storage allocated ! for WORK or IWORK. ! IERR = 2 => Method failed to converge in ! ITMAX steps. ! IERR = 3 => Error in user input. Check input ! value of N, ITOL. ! IERR = 4 => User error tolerance set too tight. ! Reset to 500.0*D1MACH(3). Iteration proceeded. ! IERR = 5 => Preconditioning matrix, M, is not ! Positive Definite. $(r,z) < 0.0$. ! IERR = 6 => Matrix A is not Positive Definite. ! $(p,Ap) < 0.0$. ! IERR = 7 => Incomplete factorization broke down ! and was fudged. Resulting preconditioning may ! be less than the best. ! IUNIT :IN Integer. ! Unit number on which to write the error at each iteration, ! if this is desired for monitoring convergence. If unit ! number is 0, no writing will occur. ! RWORK :WORK Double Precision RWORK(LENW). ! Double Precision array used for workspace. NEL is the number ! of non- ! zeros in the lower triangle of the matrix (including the ! diagonal). NU is the number of nonzeros in the upper ! triangle of the matrix (including the diagonal). ! LENW :IN Integer. ! Length of the double precision workspace, RWORK. ! LENW >= NEL+NU+8*N. ! IWORK :WORK Integer IWORK(LENIW). ! Integer array used for workspace. NEL is the number of non- ! zeros in the lower triangle of the matrix (including the ! diagonal). NU is the number of nonzeros i