*DECK DPIGMR SUBROUTINE DPIGMR (N, R0, SR, SZ, JSCAL, MAXL, MAXLP1, KMP, NRSTS, + JPRE, MATVEC, MSOLVE, NMSL, Z, V, HES, Q, LGMR, RPAR, IPAR, WK, + DL, RHOL, NRMAX, B, BNRM, X, XL, ITOL, TOL, NELT, IA, JA, A, + ISYM, IUNIT, IFLAG, ERR) C***BEGIN PROLOGUE DPIGMR C***SUBSIDIARY C***PURPOSE Internal routine for DGMRES. C***LIBRARY SLATEC (SLAP) C***CATEGORY D2A4, D2B4 C***TYPE DOUBLE PRECISION (SPIGMR-S, DPIGMR-D) C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION, C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov C Hindmarsh, Alan, (LLNL), alanh@llnl.gov C Seager, Mark K., (LLNL), seager@llnl.gov C Lawrence Livermore National Laboratory C PO Box 808, L-60 C Livermore, CA 94550 (510) 423-3141 C***DESCRIPTION C This routine solves the linear system A * Z = R0 using a C scaled preconditioned version of the generalized minimum C residual method. An initial guess of Z = 0 is assumed. C C *Usage: C INTEGER N, JSCAL, MAXL, MAXLP1, KMP, NRSTS, JPRE, NMSL, LGMR C INTEGER IPAR(USER DEFINED), NRMAX, ITOL, NELT, IA(NELT), JA(NELT) C INTEGER ISYM, IUNIT, IFLAG C DOUBLE PRECISION R0(N), SR(N), SZ(N), Z(N), V(N,MAXLP1), C $ HES(MAXLP1,MAXL), Q(2*MAXL), RPAR(USER DEFINED), C $ WK(N), DL(N), RHOL, B(N), BNRM, X(N), XL(N), C $ TOL, A(NELT), ERR C EXTERNAL MATVEC, MSOLVE C C CALL DPIGMR(N, R0, SR, SZ, JSCAL, MAXL, MAXLP1, KMP, C $ NRSTS, JPRE, MATVEC, MSOLVE, NMSL, Z, V, HES, Q, LGMR, C $ RPAR, IPAR, WK, DL, RHOL, NRMAX, B, BNRM, X, XL, C $ ITOL, TOL, NELT, IA, JA, A, ISYM, IUNIT, IFLAG, ERR) C C *Arguments: C N :IN Integer C The order of the matrix A, and the lengths C of the vectors SR, SZ, R0 and Z. C R0 :IN Double Precision R0(N) C R0 = the right hand side of the system A*Z = R0. C R0 is also used as workspace when computing C the final approximation. C (R0 is the same as V(*,MAXL+1) in the call to DPIGMR.) C SR :IN Double Precision SR(N) C SR is a vector of length N containing the non-zero C elements of the diagonal scaling matrix for R0. C SZ :IN Double Precision SZ(N) C SZ is a vector of length N containing the non-zero C elements of the diagonal scaling matrix for Z. C JSCAL :IN Integer C A flag indicating whether arrays SR and SZ are used. C JSCAL=0 means SR and SZ are not used and the C algorithm will perform as if all C SR(i) = 1 and SZ(i) = 1. C JSCAL=1 means only SZ is used, and the algorithm C performs as if all SR(i) = 1. C JSCAL=2 means only SR is used, and the algorithm C performs as if all SZ(i) = 1. C JSCAL=3 means both SR and SZ are used. C MAXL :IN Integer C The maximum allowable order of the matrix H. C MAXLP1 :IN Integer C MAXPL1 = MAXL + 1, used for dynamic dimensioning of HES. C KMP :IN Integer C The number of previous vectors the new vector VNEW C must be made orthogonal to. (KMP .le. MAXL) C NRSTS :IN Integer C Counter for the number of restarts on the current C call to DGMRES. If NRSTS .gt. 0, then the residual C R0 is already scaled, and so scaling of it is C not necessary. C JPRE :IN Integer C Preconditioner type flag. C MATVEC :EXT External. C Name of a routine which performs the matrix vector multiply C Y = A*X given A and X. The name of the MATVEC routine must C be declared external in the calling program. The calling C sequence to MATVEC is: C CALL MATVEC(N, X, Y, NELT, IA, JA, A, ISYM) C where N is the number of unknowns, Y is the product A*X C upon return, X is an input vector, and NELT is the number of C non-zeros in the SLAP IA, JA, A storage for the matrix A. C ISYM is a flag which, if non-zero, denotes that A is C symmetric and only the lower or upper triangle is stored. C MSOLVE :EXT External. C Name of the routine which solves a linear system Mz = r for C z given r with the preconditioning matrix M (M is supplied via C RPAR and IPAR arrays. The name of the MSOLVE routine must C be declared external in the calling program. The calling C sequence to MSOLVE is: C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RPAR, IPAR) C Where N is the number of unknowns, R is the right-hand side C vector and Z is the solution upon return. NELT, IA, JA, A and C ISYM are defined as below. RPAR is a double precision array C that can be used to pass necessary preconditioning information C and/or workspace to MSOLVE. IPAR is an integer work array C for the same purpose as RPAR. C NMSL :OUT Integer C The number of calls to MSOLVE. C Z :OUT Double Precision Z(N) C The final computed approximation to the solution C of the system A*Z = R0. C V :OUT Double Precision V(N,MAXLP1) C The N by (LGMR+1) array containing the LGMR C orthogonal vectors V(*,1) to V(*,LGMR). C HES :OUT Double Precision HES(MAXLP1,MAXL) C The upper triangular factor of the QR decomposition C of the (LGMR+1) by LGMR upper Hessenberg matrix whose C entries are the scaled inner-products of A*V(*,I) C and V(*,K). C Q :OUT Double Precision Q(2*MAXL) C A double precision array of length 2*MAXL containing the C components of the Givens rotations used in the QR C decomposition of HES. It is loaded in DHEQR and used in C DHELS. C LGMR :OUT Integer C The number of iterations performed and C the current order of the upper Hessenberg C matrix HES. C RPAR :IN Double Precision RPAR(USER DEFINED) C Double Precision workspace passed directly to the MSOLVE C routine. C IPAR :IN Integer IPAR(USER DEFINED) C Integer workspace passed directly to the MSOLVE routine. C WK :IN Double Precision WK(N) C A double precision work array of length N used by routines C MATVEC and MSOLVE. C DL :INOUT Double Precision DL(N) C On input, a double precision work array of length N used for C calculation of the residual norm RHO when the method is C incomplete (KMP.lt.MAXL), and/or when using restarting. C On output, the scaled residual vector RL. It is only loaded C when performing restarts of the Krylov iteration. C RHOL :OUT Double Precision C A double precision scalar containing the norm of the final C residual. C NRMAX :IN Integer C The maximum number of restarts of the Krylov iteration. C NRMAX .gt. 0 means restarting is active, while C NRMAX = 0 means restarting is not being used. C B :IN Double Precision B(N) C The right hand side of the linear system A*X = b. C BNRM :IN Double Precision C The scaled norm of b. C X :IN Double Precision X(N) C The current approximate solution as of the last C restart. C XL :IN Double Precision XL(N) C An array of length N used to hold the approximate C solution X(L) when ITOL=11. C ITOL :IN Integer C A flag to indicate the type of convergence criterion C used. See the driver for its description. C TOL :IN Double Precision C The tolerance on residuals R0-A*Z in scaled norm. C NELT :IN Integer C The length of arrays IA, JA and A. C IA :IN Integer IA(NELT) C An integer array of length NELT containing matrix data. C It is passed directly to the MATVEC and MSOLVE routines. C JA :IN Integer JA(NELT) C An integer array of length NELT containing matrix data. C It is passed directly to the MATVEC and MSOLVE routines. C A :IN Double Precision A(NELT) C A double precision array of length NELT containing matrix C data. It is passed directly to the MATVEC and MSOLVE routines. C ISYM :IN Integer C A flag to indicate symmetric matrix storage. C If ISYM=0, all non-zero entries of the matrix are C stored. If ISYM=1, the matrix is symmetric and C only the upper or lower triangular part is stored. C IUNIT :IN Integer C The i/o unit number for writing intermediate residual C norm values. C IFLAG :OUT Integer C An integer error flag.. C 0 means convergence in LGMR iterations, LGMR.le.MAXL. C 1 means the convergence test did not pass in MAXL C iterations, but the residual norm is .lt. norm(R0), C and so Z is computed. C 2 means the convergence test did not pass in MAXL C iterations, residual .ge. norm(R0), and Z = 0. C ERR :OUT Double Precision. C Error estimate of error in final approximate solution, as C defined by ITOL. C C *Cautions: C This routine will attempt to write to the Fortran logical output C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that C this logical unit is attached to a file or terminal before calling C this routine with a non-zero value for IUNIT. This routine does C not check for the validity of a non-zero IUNIT unit number. C C***SEE ALSO DGMRES C***ROUTINES CALLED DAXPY, DCOPY, DHELS, DHEQR, DNRM2, DORTH, DRLCAL, C DSCAL, ISDGMR C***REVISION HISTORY (YYMMDD) C 890404 DATE WRITTEN C 890404 Previous REVISION DATE C 890915 Made changes requested at July 1989 CML Meeting. (MKS) C 890922 Numerous changes to prologue to make closer to SLATEC C standard. (FNF) C 890929 Numerous changes to reduce SP/DP differences. (FNF) C 910411 Prologue converted to Version 4.0 format. (BAB) C 910502 Removed MATVEC and MSOLVE from ROUTINES CALLED list. (FNF) C 910506 Made subsidiary to DGMRES. (FNF) C 920511 Added complete declaration section. (WRB) C***END PROLOGUE DPIGMR