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 DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR, 00006 * SAVR,DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON, 00007 * S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW) 00008 C 00009 C***BEGIN PROLOGUE DNSK 00010 C***REFER TO DDASPK 00011 C***DATE WRITTEN 891219 (YYMMDD) 00012 C***REVISION DATE 900926 (YYMMDD) 00013 C***REVISION DATE 950126 (YYMMDD) 00014 C 00015 C 00016 C----------------------------------------------------------------------- 00017 C***DESCRIPTION 00018 C 00019 C DNSK solves a nonlinear system of 00020 C algebraic equations of the form 00021 C G(X,Y,YPRIME) = 0 for the unknown Y. 00022 C 00023 C The method used is a modified Newton scheme. 00024 C 00025 C The parameters represent 00026 C 00027 C X -- Independent variable. 00028 C Y -- Solution vector. 00029 C YPRIME -- Derivative of solution vector. 00030 C NEQ -- Number of unknowns. 00031 C RES -- External user-supplied subroutine 00032 C to evaluate the residual. See RES description 00033 C in DDASPK prologue. 00034 C PSOL -- External user-supplied routine to solve 00035 C a linear system using preconditioning. 00036 C See explanation inside DDASPK. 00037 C WT -- Vector of weights for error criterion. 00038 C RPAR,IPAR -- Real and integer arrays used for communication 00039 C between the calling program and external user 00040 C routines. They are not altered within DASPK. 00041 C SAVR -- Work vector for DNSK of length NEQ. 00042 C DELTA -- Work vector for DNSK of length NEQ. 00043 C E -- Error accumulation vector for DNSK of length NEQ. 00044 C WM,IWM -- Real and integer arrays storing 00045 C matrix information such as the matrix 00046 C of partial derivatives, permutation 00047 C vector, and various other information. 00048 C CJ -- Parameter always proportional to 1/H (step size). 00049 C SQRTN -- Square root of NEQ. 00050 C RSQRTN -- reciprical of square root of NEQ. 00051 C EPLIN -- Tolerance for linear system solver. 00052 C EPCON -- Tolerance to test for convergence of the Newton 00053 C iteration. 00054 C S -- Used for error convergence tests. 00055 C In the Newton iteration: S = RATE/(1.D0-RATE), 00056 C where RATE is the estimated rate of convergence 00057 C of the Newton iteration. 00058 C 00059 C The closer RATE is to 0., the faster the Newton 00060 C iteration is converging; the closer RATE is to 1., 00061 C the slower the Newton iteration is converging. 00062 C 00063 C The calling routine sends the initial value 00064 C of S to the Newton iteration. 00065 C CONFAC -- A residual scale factor to improve convergence. 00066 C TOLNEW -- Tolerance on the norm of Newton correction in 00067 C alternative Newton convergence test. 00068 C MULDEL -- A flag indicating whether or not to multiply 00069 C DELTA by CONFAC. 00070 C 0 ==> do not scale DELTA by CONFAC. 00071 C 1 ==> scale DELTA by CONFAC. 00072 C MAXIT -- Maximum allowed number of Newton iterations. 00073 C IRES -- Error flag returned from RES. See RES description 00074 C in DDASPK prologue. If IRES = -1, then IERNEW 00075 C will be set to 1. 00076 C If IRES < -1, then IERNEW will be set to -1. 00077 C IERSL -- Error flag for linear system solver. 00078 C See IERSL description in subroutine DSLVK. 00079 C If IERSL = 1, then IERNEW will be set to 1. 00080 C If IERSL < 0, then IERNEW will be set to -1. 00081 C IERNEW -- Error flag for Newton iteration. 00082 C 0 ==> Newton iteration converged. 00083 C 1 ==> recoverable error inside Newton iteration. 00084 C -1 ==> unrecoverable error inside Newton iteration. 00085 C----------------------------------------------------------------------- 00086 C 00087 C***ROUTINES CALLED 00088 C RES, DSLVK, DDWNRM 00089 C 00090 C***END PROLOGUE DNSK 00091 C 00092 C 00093 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00094 DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*),SAVR(*) 00095 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) 00096 EXTERNAL RES, PSOL 00097 C 00098 PARAMETER (LNNI=19, LNRE=12) 00099 C 00100 C Initialize Newton counter M and accumulation vector E. 00101 C 00102 M = 0 00103 DO 100 I=1,NEQ 00104 100 E(I) = 0.0D0 00105 C 00106 C Corrector loop. 00107 C 00108 300 CONTINUE 00109 IWM(LNNI) = IWM(LNNI) + 1 00110 C 00111 C If necessary, multiply residual by convergence factor. 00112 C 00113 IF (MULDEL .EQ. 1) THEN 00114 DO 320 I = 1,NEQ 00115 320 DELTA(I) = DELTA(I) * CONFAC 00116 ENDIF 00117 C 00118 C Save residual in SAVR. 00119 C 00120 DO 340 I = 1,NEQ 00121 340 SAVR(I) = DELTA(I) 00122 C 00123 C Compute a new iterate. Store the correction in DELTA. 00124 C 00125 CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM, 00126 * RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK, 00127 * RPAR, IPAR) 00128 IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 380 00129 C 00130 C Update Y, E, and YPRIME. 00131 C 00132 DO 360 I=1,NEQ 00133 Y(I) = Y(I) - DELTA(I) 00134 E(I) = E(I) - DELTA(I) 00135 360 YPRIME(I) = YPRIME(I) - CJ*DELTA(I) 00136 C 00137 C Test for convergence of the iteration. 00138 C 00139 DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR) 00140 IF (DELNRM .LE. TOLNEW) GO TO 370 00141 IF (M .EQ. 0) THEN 00142 OLDNRM = DELNRM 00143 ELSE 00144 RATE = (DELNRM/OLDNRM)**(1.0D0/M) 00145 IF (RATE .GT. 0.9D0) GO TO 380 00146 S = RATE/(1.0D0 - RATE) 00147 ENDIF 00148 IF (S*DELNRM .LE. EPCON) GO TO 370 00149 C 00150 C The corrector has not yet converged. Update M and test whether 00151 C the maximum number of iterations have been tried. 00152 C 00153 M = M + 1 00154 IF (M .GE. MAXIT) GO TO 380 00155 C 00156 C Evaluate the residual, and go back to do another iteration. 00157 C 00158 IWM(LNRE) = IWM(LNRE) + 1 00159 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 00160 IF (IRES .LT. 0) GO TO 380 00161 GO TO 300 00162 C 00163 C The iteration has converged. 00164 C 00165 370 RETURN 00166 C 00167 C The iteration has not converged. Set IERNEW appropriately. 00168 C 00169 380 CONTINUE 00170 IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN 00171 IERNEW = -1 00172 ELSE 00173 IERNEW = 1 00174 ENDIF 00175 RETURN 00176 C 00177 C 00178 C------END OF SUBROUTINE DNSK------------------------------------------- 00179 END
1.7.1
