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 DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR, 00006 * SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON, 00007 * RATEMX,MAXIT,STPTOL,ICNFLG,ICNSTR,IERNEW) 00008 C 00009 C***BEGIN PROLOGUE DNSIK 00010 C***REFER TO DDASPK 00011 C***DATE WRITTEN 940701 (YYMMDD) 00012 C***REVISION DATE 950714 (YYMMDD) 00013 C 00014 C 00015 C----------------------------------------------------------------------- 00016 C***DESCRIPTION 00017 C 00018 C DNSIK solves a nonlinear system of algebraic equations of the 00019 C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in 00020 C the initial conditions. 00021 C 00022 C The method used is a Newton scheme combined with a linesearch 00023 C algorithm, using Krylov iterative linear system methods. 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 ICOPT -- Initial condition option chosen (1 or 2). 00032 C ID -- Array of dimension NEQ, which must be initialized 00033 C if ICOPT = 1. See DDASIC. 00034 C RES -- External user-supplied subroutine 00035 C to evaluate the residual. See RES description 00036 C in DDASPK prologue. 00037 C PSOL -- External user-supplied routine to solve 00038 C a linear system using preconditioning. 00039 C See explanation inside DDASPK. 00040 C WT -- Vector of weights for error criterion. 00041 C RPAR,IPAR -- Real and integer arrays used for communication 00042 C between the calling program and external user 00043 C routines. They are not altered within DASPK. 00044 C SAVR -- Work vector for DNSIK of length NEQ. 00045 C DELTA -- Residual vector on entry, and work vector of 00046 C length NEQ for DNSIK. 00047 C R -- Work vector for DNSIK of length NEQ. 00048 C YIC,YPIC -- Work vectors for DNSIK, each of length NEQ. 00049 C PWK -- Work vector for DNSIK of length NEQ. 00050 C WM,IWM -- Real and integer arrays storing 00051 C matrix information such as the matrix 00052 C of partial derivatives, permutation 00053 C vector, and various other information. 00054 C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2). 00055 C SQRTN -- Square root of NEQ. 00056 C RSQRTN -- reciprical of square root of NEQ. 00057 C EPLIN -- Tolerance for linear system solver. 00058 C EPCON -- Tolerance to test for convergence of the Newton 00059 C iteration. 00060 C RATEMX -- Maximum convergence rate for which Newton iteration 00061 C is considered converging. 00062 C MAXIT -- Maximum allowed number of Newton iterations. 00063 C STPTOL -- Tolerance used in calculating the minimum lambda 00064 C value allowed. 00065 C ICNFLG -- Integer scalar. If nonzero, then constraint 00066 C violations in the proposed new approximate solution 00067 C will be checked for, and the maximum step length 00068 C will be adjusted accordingly. 00069 C ICNSTR -- Integer array of length NEQ containing flags for 00070 C checking constraints. 00071 C IERNEW -- Error flag for Newton iteration. 00072 C 0 ==> Newton iteration converged. 00073 C 1 ==> failed to converge, but RATE .lt. 1. 00074 C 2 ==> failed to converge, RATE .gt. RATEMX. 00075 C 3 ==> other recoverable error. 00076 C -1 ==> unrecoverable error inside Newton iteration. 00077 C----------------------------------------------------------------------- 00078 C 00079 C***ROUTINES CALLED 00080 C DFNRMK, DSLVK, DDWNRM, DLINSK, DCOPY 00081 C 00082 C***END PROLOGUE DNSIK 00083 C 00084 C 00085 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00086 DIMENSION Y(*),YPRIME(*),WT(*),ID(*),DELTA(*),R(*),SAVR(*) 00087 DIMENSION YIC(*),YPIC(*),PWK(*),WM(*),IWM(*), RPAR(*),IPAR(*) 00088 DIMENSION ICNSTR(*) 00089 EXTERNAL RES, PSOL 00090 C 00091 PARAMETER (LNNI=19, LNPS=21, LLOCWP=29, LLCIWP=30) 00092 PARAMETER (LLSOFF=35, LSTOL=14) 00093 C 00094 C 00095 C Initializations. M is the Newton iteration counter. 00096 C 00097 LSOFF = IWM(LLSOFF) 00098 M = 0 00099 RATE = 1.0D0 00100 LWP = IWM(LLOCWP) 00101 LIWP = IWM(LLCIWP) 00102 RLX = 0.4D0 00103 C 00104 C Save residual in SAVR. 00105 C 00106 CALL DCOPY (NEQ, DELTA, 1, SAVR, 1) 00107 C 00108 C Compute norm of (P-inverse)*(residual). 00109 C 00110 CALL DFNRMK (NEQ, Y, X, YPRIME, SAVR, R, CJ, WT, SQRTN, RSQRTN, 00111 * RES, IRES, PSOL, 1, IER, FNRM, EPLIN, WM(LWP), IWM(LIWP), 00112 * PWK, RPAR, IPAR) 00113 IWM(LNPS) = IWM(LNPS) + 1 00114 IF (IER .NE. 0) THEN 00115 IERNEW = 3 00116 RETURN 00117 ENDIF 00118 C 00119 C Return now if residual norm is .le. EPCON. 00120 C 00121 IF (FNRM .LE. EPCON) RETURN 00122 C 00123 C Newton iteration loop. 00124 C 00125 300 CONTINUE 00126 IWM(LNNI) = IWM(LNNI) + 1 00127 C 00128 C Compute a new step vector DELTA. 00129 C 00130 CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM, 00131 * RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK, 00132 * RPAR, IPAR) 00133 IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 390 00134 C 00135 C Get norm of DELTA. Return now if DELTA is zero. 00136 C 00137 DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR) 00138 IF (DELNRM .EQ. 0.0D0) RETURN 00139 C 00140 C Call linesearch routine for global strategy and set RATE. 00141 C 00142 OLDFNM = FNRM 00143 C 00144 CALL DLINSK (NEQ, Y, X, YPRIME, SAVR, CJ, DELTA, DELNRM, WT, 00145 * SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL, WM, IWM, 00146 * RHOK, FNRM, ICOPT, ID, WM(LWP), IWM(LIWP), R, EPLIN, YIC, YPIC, 00147 * PWK, ICNFLG, ICNSTR, RLX, RPAR, IPAR) 00148 C 00149 RATE = FNRM/OLDFNM 00150 C 00151 C Check for error condition from linesearch. 00152 IF (IRET .NE. 0) GO TO 390 00153 C 00154 C Test for convergence of the iteration, and return or loop. 00155 C 00156 IF (FNRM .LE. EPCON) RETURN 00157 C 00158 C The iteration has not yet converged. Update M. 00159 C Test whether the maximum number of iterations have been tried. 00160 C 00161 M=M+1 00162 IF(M .GE. MAXIT) GO TO 380 00163 C 00164 C Copy the residual SAVR to DELTA and loop for another iteration. 00165 C 00166 CALL DCOPY (NEQ, SAVR, 1, DELTA, 1) 00167 GO TO 300 00168 C 00169 C The maximum number of iterations was done. Set IERNEW and return. 00170 C 00171 380 IF (RATE .LE. RATEMX) THEN 00172 IERNEW = 1 00173 ELSE 00174 IERNEW = 2 00175 ENDIF 00176 RETURN 00177 C 00178 390 IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN 00179 IERNEW = -1 00180 ELSE 00181 IERNEW = 3 00182 IF (IRES .EQ. 0 .AND. IERSL .EQ. 1 .AND. M .GE. 2 00183 1 .AND. RATE .LT. 1.0D0) IERNEW = 1 00184 ENDIF 00185 RETURN 00186 C 00187 C 00188 C----------------------- END OF SUBROUTINE DNSIK------------------------ 00189 END
1.7.1
