00001 C Work performed under the auspices of the U.S. Department of Energy 00002 C by Lawrence Livermore National Laboratory under contract number 00003 C W-7405-Eng-48. 00004 C 00005 SUBROUTINE DSPIGM (NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL, 00006 * MAXLP1, KMP, EPLIN, CJ, RES, IRES, NRE, PSOL, NPSL, Z, V, 00007 * HES, Q, LGMR, WP, IWP, WK, DL, RHOK, IFLAG, IRST, NRSTS, 00008 * RPAR, IPAR) 00009 C 00010 C***BEGIN PROLOGUE DSPIGM 00011 C***DATE WRITTEN 890101 (YYMMDD) 00012 C***REVISION DATE 900926 (YYMMDD) 00013 C***REVISION DATE 940927 Removed MNEWT and added RHOK in call list. 00014 C 00015 C 00016 C----------------------------------------------------------------------- 00017 C***DESCRIPTION 00018 C 00019 C This routine solves the linear system A * Z = R using a scaled 00020 C preconditioned version of the generalized minimum residual method. 00021 C An initial guess of Z = 0 is assumed. 00022 C 00023 C On entry 00024 C 00025 C NEQ = Problem size, passed to PSOL. 00026 C 00027 C TN = Current Value of T. 00028 C 00029 C Y = Array Containing current dependent variable vector. 00030 C 00031 C YPRIME = Array Containing current first derivative of Y. 00032 C 00033 C SAVR = Array containing current value of G(T,Y,YPRIME). 00034 C 00035 C R = The right hand side of the system A*Z = R. 00036 C R is also used as work space when computing 00037 C the final approximation and will therefore be 00038 C destroyed. 00039 C (R is the same as V(*,MAXL+1) in the call to DSPIGM.) 00040 C 00041 C WGHT = The vector of length NEQ containing the nonzero 00042 C elements of the diagonal scaling matrix. 00043 C 00044 C MAXL = The maximum allowable order of the matrix H. 00045 C 00046 C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES. 00047 C 00048 C KMP = The number of previous vectors the new vector, VNEW, 00049 C must be made orthogonal to. (KMP .LE. MAXL.) 00050 C 00051 C EPLIN = Tolerance on residuals R-A*Z in weighted rms norm. 00052 C 00053 C CJ = Scalar proportional to current value of 00054 C 1/(step size H). 00055 C 00056 C WK = Real work array used by routine DATV and PSOL. 00057 C 00058 C DL = Real work array used for calculation of the residual 00059 C norm RHO when the method is incomplete (KMP.LT.MAXL) 00060 C and/or when using restarting. 00061 C 00062 C WP = Real work array used by preconditioner PSOL. 00063 C 00064 C IWP = Integer work array used by preconditioner PSOL. 00065 C 00066 C IRST = Method flag indicating if restarting is being 00067 C performed. IRST .GT. 0 means restarting is active, 00068 C while IRST = 0 means restarting is not being used. 00069 C 00070 C NRSTS = Counter for the number of restarts on the current 00071 C call to DSPIGM. If NRSTS .GT. 0, then the residual 00072 C R is already scaled, and so scaling of R is not 00073 C necessary. 00074 C 00075 C 00076 C On Return 00077 C 00078 C Z = The final computed approximation to the solution 00079 C of the system A*Z = R. 00080 C 00081 C LGMR = The number of iterations performed and 00082 C the current order of the upper Hessenberg 00083 C matrix HES. 00084 C 00085 C NRE = The number of calls to RES (i.e. DATV) 00086 C 00087 C NPSL = The number of calls to PSOL. 00088 C 00089 C V = The neq by (LGMR+1) array containing the LGMR 00090 C orthogonal vectors V(*,1) to V(*,LGMR). 00091 C 00092 C HES = The upper triangular factor of the QR decomposition 00093 C of the (LGMR+1) by LGMR upper Hessenberg matrix whose 00094 C entries are the scaled inner-products of A*V(*,I) 00095 C and V(*,K). 00096 C 00097 C Q = Real array of length 2*MAXL containing the components 00098 C of the givens rotations used in the QR decomposition 00099 C of HES. It is loaded in DHEQR and used in DHELS. 00100 C 00101 C IRES = Error flag from RES. 00102 C 00103 C DL = Scaled preconditioned residual, 00104 C (D-inverse)*(P-inverse)*(R-A*Z). Only loaded when 00105 C performing restarts of the Krylov iteration. 00106 C 00107 C RHOK = Weighted norm of final preconditioned residual. 00108 C 00109 C IFLAG = Integer error flag.. 00110 C 0 Means convergence in LGMR iterations, LGMR.LE.MAXL. 00111 C 1 Means the convergence test did not pass in MAXL 00112 C iterations, but the new residual norm (RHO) is 00113 C .LT. the old residual norm (RNRM), and so Z is 00114 C computed. 00115 C 2 Means the convergence test did not pass in MAXL 00116 C iterations, new residual norm (RHO) .GE. old residual 00117 C norm (RNRM), and the initial guess, Z = 0, is 00118 C returned. 00119 C 3 Means there was a recoverable error in PSOL 00120 C caused by the preconditioner being out of date. 00121 C -1 Means there was an unrecoverable error in PSOL. 00122 C 00123 C----------------------------------------------------------------------- 00124 C***ROUTINES CALLED 00125 C PSOL, DNRM2, DSCAL, DATV, DORTH, DHEQR, DCOPY, DHELS, DAXPY 00126 C 00127 C***END PROLOGUE DSPIGM 00128 C 00129 INTEGER NEQ,MAXL,MAXLP1,KMP,IRES,NRE,NPSL,LGMR,IWP, 00130 1 IFLAG,IRST,NRSTS,IPAR 00131 DOUBLE PRECISION TN,Y,YPRIME,SAVR,R,WGHT,EPLIN,CJ,Z,V,HES,Q,WP,WK, 00132 1 DL,RHOK,RPAR 00133 DIMENSION Y(*), YPRIME(*), SAVR(*), R(*), WGHT(*), Z(*), 00134 1 V(NEQ,*), HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*), 00135 2 RPAR(*), IPAR(*) 00136 INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1 00137 DOUBLE PRECISION RNRM,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM 00138 EXTERNAL RES, PSOL 00139 C 00140 IER = 0 00141 IFLAG = 0 00142 LGMR = 0 00143 NPSL = 0 00144 NRE = 0 00145 C----------------------------------------------------------------------- 00146 C The initial guess for Z is 0. The initial residual is therefore 00147 C the vector R. Initialize Z to 0. 00148 C----------------------------------------------------------------------- 00149 DO 10 I = 1,NEQ 00150 10 Z(I) = 0.0D0 00151 C----------------------------------------------------------------------- 00152 C Apply inverse of left preconditioner to vector R if NRSTS .EQ. 0. 00153 C Form V(*,1), the scaled preconditioned right hand side. 00154 C----------------------------------------------------------------------- 00155 IF (NRSTS .EQ. 0) THEN 00156 CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, WK, CJ, WGHT, WP, IWP, 00157 1 R, EPLIN, IER, RPAR, IPAR) 00158 NPSL = 1 00159 IF (IER .NE. 0) GO TO 300 00160 DO 30 I = 1,NEQ 00161 30 V(I,1) = R(I)*WGHT(I) 00162 ELSE 00163 DO 35 I = 1,NEQ 00164 35 V(I,1) = R(I) 00165 ENDIF 00166 C----------------------------------------------------------------------- 00167 C Calculate norm of scaled vector V(*,1) and normalize it 00168 C If, however, the norm of V(*,1) (i.e. the norm of the preconditioned 00169 C residual) is .le. EPLIN, then return with Z=0. 00170 C----------------------------------------------------------------------- 00171 RNRM = DNRM2 (NEQ, V, 1) 00172 IF (RNRM .LE. EPLIN) THEN 00173 RHOK = RNRM 00174 RETURN 00175 ENDIF 00176 TEM = 1.0D0/RNRM 00177 CALL DSCAL (NEQ, TEM, V(1,1), 1) 00178 C----------------------------------------------------------------------- 00179 C Zero out the HES array. 00180 C----------------------------------------------------------------------- 00181 DO 65 J = 1,MAXL 00182 DO 60 I = 1,MAXLP1 00183 60 HES(I,J) = 0.0D0 00184 65 CONTINUE 00185 C----------------------------------------------------------------------- 00186 C Main loop to compute the vectors V(*,2) to V(*,MAXL). 00187 C The running product PROD is needed for the convergence test. 00188 C----------------------------------------------------------------------- 00189 PROD = 1.0D0 00190 DO 90 LL = 1,MAXL 00191 LGMR = LL 00192 C----------------------------------------------------------------------- 00193 C Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is 00194 C the matrix A with scaling and inverse preconditioner factors applied. 00195 C Call routine DORTH to orthogonalize the new vector VNEW = V(*,LL+1). 00196 C call routine DHEQR to update the factors of HES. 00197 C----------------------------------------------------------------------- 00198 CALL DATV (NEQ, Y, TN, YPRIME, SAVR, V(1,LL), WGHT, Z, 00199 1 RES, IRES, PSOL, V(1,LL+1), WK, WP, IWP, CJ, EPLIN, 00200 1 IER, NRE, NPSL, RPAR, IPAR) 00201 IF (IRES .LT. 0) RETURN 00202 IF (IER .NE. 0) GO TO 300 00203 CALL DORTH (V(1,LL+1), V, HES, NEQ, LL, MAXLP1, KMP, SNORMW) 00204 HES(LL+1,LL) = SNORMW 00205 CALL DHEQR (HES, MAXLP1, LL, Q, INFO, LL) 00206 IF (INFO .EQ. LL) GO TO 120 00207 C----------------------------------------------------------------------- 00208 C Update RHO, the estimate of the norm of the residual R - A*ZL. 00209 C If KMP .LT. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not 00210 C necessarily orthogonal for LL .GT. KMP. The vector DL must then 00211 C be computed, and its norm used in the calculation of RHO. 00212 C----------------------------------------------------------------------- 00213 PROD = PROD*Q(2*LL) 00214 RHO = ABS(PROD*RNRM) 00215 IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN 00216 IF (LL .EQ. KMP+1) THEN 00217 CALL DCOPY (NEQ, V(1,1), 1, DL, 1) 00218 DO 75 I = 1,KMP 00219 IP1 = I + 1 00220 I2 = I*2 00221 S = Q(I2) 00222 C = Q(I2-1) 00223 DO 70 K = 1,NEQ 00224 70 DL(K) = S*DL(K) + C*V(K,IP1) 00225 75 CONTINUE 00226 ENDIF 00227 S = Q(2*LL) 00228 C = Q(2*LL-1)/SNORMW 00229 LLP1 = LL + 1 00230 DO 80 K = 1,NEQ 00231 80 DL(K) = S*DL(K) + C*V(K,LLP1) 00232 DLNRM = DNRM2 (NEQ, DL, 1) 00233 RHO = RHO*DLNRM 00234 ENDIF 00235 C----------------------------------------------------------------------- 00236 C Test for convergence. If passed, compute approximation ZL. 00237 C If failed and LL .LT. MAXL, then continue iterating. 00238 C----------------------------------------------------------------------- 00239 IF (RHO .LE. EPLIN) GO TO 200 00240 IF (LL .EQ. MAXL) GO TO 100 00241 C----------------------------------------------------------------------- 00242 C Rescale so that the norm of V(1,LL+1) is one. 00243 C----------------------------------------------------------------------- 00244 TEM = 1.0D0/SNORMW 00245 CALL DSCAL (NEQ, TEM, V(1,LL+1), 1) 00246 90 CONTINUE 00247 100 CONTINUE 00248 IF (RHO .LT. RNRM) GO TO 150 00249 120 CONTINUE 00250 IFLAG = 2 00251 DO 130 I = 1,NEQ 00252 130 Z(I) = 0.D0 00253 RETURN 00254 150 IFLAG = 1 00255 C----------------------------------------------------------------------- 00256 C The tolerance was not met, but the residual norm was reduced. 00257 C If performing restarting (IRST .gt. 0) calculate the residual vector 00258 C RL and store it in the DL array. If the incomplete version is 00259 C being used (KMP .lt. MAXL) then DL has already been calculated. 00260 C----------------------------------------------------------------------- 00261 IF (IRST .GT. 0) THEN 00262 IF (KMP .EQ. MAXL) THEN 00263 C 00264 C Calculate DL from the V(I)'s. 00265 C 00266 CALL DCOPY (NEQ, V(1,1), 1, DL, 1) 00267 MAXLM1 = MAXL - 1 00268 DO 175 I = 1,MAXLM1 00269 IP1 = I + 1 00270 I2 = I*2 00271 S = Q(I2) 00272 C = Q(I2-1) 00273 DO 170 K = 1,NEQ 00274 170 DL(K) = S*DL(K) + C*V(K,IP1) 00275 175 CONTINUE 00276 S = Q(2*MAXL) 00277 C = Q(2*MAXL-1)/SNORMW 00278 DO 180 K = 1,NEQ 00279 180 DL(K) = S*DL(K) + C*V(K,MAXLP1) 00280 ENDIF 00281 C 00282 C Scale DL by RNRM*PROD to obtain the residual RL. 00283 C 00284 TEM = RNRM*PROD 00285 CALL DSCAL(NEQ, TEM, DL, 1) 00286 ENDIF 00287 C----------------------------------------------------------------------- 00288 C Compute the approximation ZL to the solution. 00289 C Since the vector Z was used as work space, and the initial guess 00290 C of the Newton correction is zero, Z must be reset to zero. 00291 C----------------------------------------------------------------------- 00292 200 CONTINUE 00293 LL = LGMR 00294 LLP1 = LL + 1 00295 DO 210 K = 1,LLP1 00296 210 R(K) = 0.0D0 00297 R(1) = RNRM 00298 CALL DHELS (HES, MAXLP1, LL, Q, R) 00299 DO 220 K = 1,NEQ 00300 220 Z(K) = 0.0D0 00301 DO 230 I = 1,LL 00302 CALL DAXPY (NEQ, R(I), V(1,I), 1, Z, 1) 00303 230 CONTINUE 00304 DO 240 I = 1,NEQ 00305 240 Z(I) = Z(I)/WGHT(I) 00306 C Load RHO into RHOK. 00307 RHOK = RHO 00308 RETURN 00309 C----------------------------------------------------------------------- 00310 C This block handles error returns forced by routine PSOL. 00311 C----------------------------------------------------------------------- 00312 300 CONTINUE 00313 IF (IER .LT. 0) IFLAG = -1 00314 IF (IER .GT. 0) IFLAG = 3 00315 C 00316 RETURN 00317 C 00318 C------END OF SUBROUTINE DSPIGM----------------------------------------- 00319 END
1.7.1
