1*DECK DXLCAL 2 SUBROUTINE DXLCAL (N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM, 3 + WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR, NELT, IA, JA, A, 4 + ISYM) 5C***BEGIN PROLOGUE DXLCAL 6C***SUBSIDIARY 7C***PURPOSE Internal routine for DGMRES. 8C***LIBRARY SLATEC (SLAP) 9C***CATEGORY D2A4, D2B4 10C***TYPE DOUBLE PRECISION (SXLCAL-S, DXLCAL-D) 11C***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION, 12C NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE 13C***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov 14C Hindmarsh, Alan, (LLNL), alanh@llnl.gov 15C Seager, Mark K., (LLNL), seager@llnl.gov 16C Lawrence Livermore National Laboratory 17C PO Box 808, L-60 18C Livermore, CA 94550 (510) 423-3141 19C***DESCRIPTION 20C This routine computes the solution XL, the current DGMRES 21C iterate, given the V(I)'s and the QR factorization of the 22C Hessenberg matrix HES. This routine is only called when 23C ITOL=11. 24C 25C *Usage: 26C INTEGER N, LGMR, MAXLP1, JSCAL, JPRE, NMSL, IPAR(USER DEFINED) 27C INTEGER NELT, IA(NELT), JA(NELT), ISYM 28C DOUBLE PRECISION X(N), XL(N), ZL(N), HES(MAXLP1,MAXL), Q(2*MAXL), 29C $ V(N,MAXLP1), R0NRM, WK(N), SZ(N), 30C $ RPAR(USER DEFINED), A(NELT) 31C EXTERNAL MSOLVE 32C 33C CALL DXLCAL(N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM, 34C $ WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR, 35C $ NELT, IA, JA, A, ISYM) 36C 37C *Arguments: 38C N :IN Integer 39C The order of the matrix A, and the lengths 40C of the vectors SR, SZ, R0 and Z. 41C LGMR :IN Integer 42C The number of iterations performed and 43C the current order of the upper Hessenberg 44C matrix HES. 45C X :IN Double Precision X(N) 46C The current approximate solution as of the last restart. 47C XL :OUT Double Precision XL(N) 48C An array of length N used to hold the approximate 49C solution X(L). 50C Warning: XL and ZL are the same array in the calling routine. 51C ZL :IN Double Precision ZL(N) 52C An array of length N used to hold the approximate 53C solution Z(L). 54C HES :IN Double Precision HES(MAXLP1,MAXL) 55C The upper triangular factor of the QR decomposition 56C of the (LGMR+1) by LGMR upper Hessenberg matrix whose 57C entries are the scaled inner-products of A*V(*,i) and V(*,k). 58C MAXLP1 :IN Integer 59C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES. 60C MAXL is the maximum allowable order of the matrix HES. 61C Q :IN Double Precision Q(2*MAXL) 62C A double precision array of length 2*MAXL containing the 63C components of the Givens rotations used in the QR 64C decomposition of HES. It is loaded in DHEQR. 65C V :IN Double Precision V(N,MAXLP1) 66C The N by(LGMR+1) array containing the LGMR 67C orthogonal vectors V(*,1) to V(*,LGMR). 68C R0NRM :IN Double Precision 69C The scaled norm of the initial residual for the 70C current call to DPIGMR. 71C WK :IN Double Precision WK(N) 72C A double precision work array of length N. 73C SZ :IN Double Precision SZ(N) 74C A vector of length N containing the non-zero 75C elements of the diagonal scaling matrix for Z. 76C JSCAL :IN Integer 77C A flag indicating whether arrays SR and SZ are used. 78C JSCAL=0 means SR and SZ are not used and the 79C algorithm will perform as if all 80C SR(i) = 1 and SZ(i) = 1. 81C JSCAL=1 means only SZ is used, and the algorithm 82C performs as if all SR(i) = 1. 83C JSCAL=2 means only SR is used, and the algorithm 84C performs as if all SZ(i) = 1. 85C JSCAL=3 means both SR and SZ are used. 86C JPRE :IN Integer 87C The preconditioner type flag. 88C MSOLVE :EXT External. 89C Name of the routine which solves a linear system Mz = r for 90C z given r with the preconditioning matrix M (M is supplied via 91C RPAR and IPAR arrays. The name of the MSOLVE routine must 92C be declared external in the calling program. The calling 93C sequence to MSOLVE is: 94C CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RPAR, IPAR) 95C Where N is the number of unknowns, R is the right-hand side 96C vector and Z is the solution upon return. NELT, IA, JA, A and 97C ISYM are defined as below. RPAR is a double precision array 98C that can be used to pass necessary preconditioning information 99C and/or workspace to MSOLVE. IPAR is an integer work array 100C for the same purpose as RPAR. 101C NMSL :IN Integer 102C The number of calls to MSOLVE. 103C RPAR :IN Double Precision RPAR(USER DEFINED) 104C Double Precision workspace passed directly to the MSOLVE 105C routine. 106C IPAR :IN Integer IPAR(USER DEFINED) 107C Integer workspace passed directly to the MSOLVE routine. 108C NELT :IN Integer 109C The length of arrays IA, JA and A. 110C IA :IN Integer IA(NELT) 111C An integer array of length NELT containing matrix data. 112C It is passed directly to the MATVEC and MSOLVE routines. 113C JA :IN Integer JA(NELT) 114C An integer array of length NELT containing matrix data. 115C It is passed directly to the MATVEC and MSOLVE routines. 116C A :IN Double Precision A(NELT) 117C A double precision array of length NELT containing matrix 118C data. 119C It is passed directly to the MATVEC and MSOLVE routines. 120C ISYM :IN Integer 121C A flag to indicate symmetric matrix storage. 122C If ISYM=0, all non-zero entries of the matrix are 123C stored. If ISYM=1, the matrix is symmetric and 124C only the upper or lower triangular part is stored. 125C 126C***SEE ALSO DGMRES 127C***ROUTINES CALLED DAXPY, DCOPY, DHELS 128C***REVISION HISTORY (YYMMDD) 129C 890404 DATE WRITTEN 130C 890404 Previous REVISION DATE 131C 890915 Made changes requested at July 1989 CML Meeting. (MKS) 132C 890922 Numerous changes to prologue to make closer to SLATEC 133C standard. (FNF) 134C 890929 Numerous changes to reduce SP/DP differences. (FNF) 135C 910411 Prologue converted to Version 4.0 format. (BAB) 136C 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF) 137C 910506 Made subsidiary to DGMRES. (FNF) 138C 920511 Added complete declaration section. (WRB) 139C***END PROLOGUE DXLCAL 140C The following is for optimized compilation on LLNL/LTSS Crays. 141CLLL. OPTIMIZE 142C .. Scalar Arguments .. 143 DOUBLE PRECISION R0NRM 144 INTEGER ISYM, JPRE, JSCAL, LGMR, MAXLP1, N, NELT, NMSL 145C .. Array Arguments .. 146 DOUBLE PRECISION A(NELT), HES(MAXLP1,*), Q(*), RPAR(*), SZ(*), 147 + V(N,*), WK(N), X(N), XL(N), ZL(N) 148 INTEGER IA(NELT), IPAR(*), JA(NELT) 149C .. Subroutine Arguments .. 150 EXTERNAL MSOLVE 151C .. Local Scalars .. 152 INTEGER I, K, LL, LLP1 153C .. External Subroutines .. 154 EXTERNAL DAXPY, DCOPY, DHELS 155C***FIRST EXECUTABLE STATEMENT DXLCAL 156 LL = LGMR 157 LLP1 = LL + 1 158 DO 10 K = 1,LLP1 159 WK(K) = 0 160 10 CONTINUE 161 WK(1) = R0NRM 162 CALL DHELS(HES, MAXLP1, LL, Q, WK) 163 DO 20 K = 1,N 164 ZL(K) = 0 165 20 CONTINUE 166 DO 30 I = 1,LL 167 CALL DAXPY(N, WK(I), V(1,I), 1, ZL, 1) 168 30 CONTINUE 169 IF ((JSCAL .EQ. 1) .OR.(JSCAL .EQ. 3)) THEN 170 DO 40 K = 1,N 171 ZL(K) = ZL(K)/SZ(K) 172 40 CONTINUE 173 ENDIF 174 IF (JPRE .GT. 0) THEN 175 CALL DCOPY(N, ZL, 1, WK, 1) 176 CALL MSOLVE(N, WK, ZL, NELT, IA, JA, A, ISYM, RPAR, IPAR) 177 NMSL = NMSL + 1 178 ENDIF 179C calculate XL from X and ZL. 180 DO 50 K = 1,N 181 XL(K) = X(K) + ZL(K) 182 50 CONTINUE 183 RETURN 184C------------- LAST LINE OF DXLCAL FOLLOWS ---------------------------- 185 END 186