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 DLINSK (NEQ, Y, T, YPRIME, SAVR, CJ, P, PNRM, WT, 00006 * SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL, WM, IWM, 00007 * RHOK, FNRM, ICOPT, ID, WP, IWP, R, EPLIN, YNEW, YPNEW, PWK, 00008 * ICNFLG, ICNSTR, RLX, RPAR, IPAR) 00009 C 00010 C***BEGIN PROLOGUE DLINSK 00011 C***REFER TO DNSIK 00012 C***DATE WRITTEN 940830 (YYMMDD) 00013 C***REVISION DATE 951006 (Arguments SQRTN, RSQRTN added.) 00014 C***REVISION DATE 960129 Moved line RL = ONE to top block. 00015 C 00016 C 00017 C----------------------------------------------------------------------- 00018 C***DESCRIPTION 00019 C 00020 C DLINSK uses a linesearch algorithm to calculate a new (Y,YPRIME) 00021 C pair (YNEW,YPNEW) such that 00022 C 00023 C f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME) + 00024 C ALPHA*RL*RHOK*RHOK , 00025 C 00026 C where 0 < RL <= 1, and RHOK is the scaled preconditioned norm of 00027 C the final residual vector in the Krylov iteration. 00028 C Here, f(y,y') is defined as 00029 C 00030 C f(y,y') = (1/2)*norm( (P-inverse)*G(t,y,y') )**2 , 00031 C 00032 C where norm() is the weighted RMS vector norm, G is the DAE 00033 C system residual function, and P is the preconditioner used 00034 C in the Krylov iteration. 00035 C 00036 C In addition to the parameters defined elsewhere, we have 00037 C 00038 C SAVR -- Work array of length NEQ, containing the residual 00039 C vector G(t,y,y') on return. 00040 C P -- Approximate Newton step used in backtracking. 00041 C PNRM -- Weighted RMS norm of P. 00042 C LSOFF -- Flag showing whether the linesearch algorithm is 00043 C to be invoked. 0 means do the linesearch, 00044 C 1 means turn off linesearch. 00045 C STPTOL -- Tolerance used in calculating the minimum lambda 00046 C value allowed. 00047 C ICNFLG -- Integer scalar. If nonzero, then constraint violations 00048 C in the proposed new approximate solution will be 00049 C checked for, and the maximum step length will be 00050 C adjusted accordingly. 00051 C ICNSTR -- Integer array of length NEQ containing flags for 00052 C checking constraints. 00053 C RHOK -- Weighted norm of preconditioned Krylov residual. 00054 C RLX -- Real scalar restricting update size in DCNSTR. 00055 C YNEW -- Array of length NEQ used to hold the new Y in 00056 C performing the linesearch. 00057 C YPNEW -- Array of length NEQ used to hold the new YPRIME in 00058 C performing the linesearch. 00059 C PWK -- Work vector of length NEQ for use in PSOL. 00060 C Y -- Array of length NEQ containing the new Y (i.e.,=YNEW). 00061 C YPRIME -- Array of length NEQ containing the new YPRIME 00062 C (i.e.,=YPNEW). 00063 C FNRM -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the 00064 C current (Y,YPRIME) on input and output. 00065 C R -- Work space length NEQ for residual vector. 00066 C IRET -- Return flag. 00067 C IRET=0 means that a satisfactory (Y,YPRIME) was found. 00068 C IRET=1 means that the routine failed to find a new 00069 C (Y,YPRIME) that was sufficiently distinct from 00070 C the current (Y,YPRIME) pair. 00071 C IRET=2 means a failure in RES or PSOL. 00072 C----------------------------------------------------------------------- 00073 C 00074 C***ROUTINES CALLED 00075 C DFNRMK, DYYPNW, DCOPY 00076 C 00077 C***END PROLOGUE DLINSK 00078 C 00079 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00080 EXTERNAL RES, PSOL 00081 DIMENSION Y(*), YPRIME(*), P(*), WT(*), SAVR(*), R(*), ID(*) 00082 DIMENSION WM(*), IWM(*), YNEW(*), YPNEW(*), PWK(*), ICNSTR(*) 00083 DIMENSION WP(*), IWP(*), RPAR(*), IPAR(*) 00084 CHARACTER MSG*80 00085 C 00086 PARAMETER (LNRE=12, LNPS=21, LKPRIN=31) 00087 C 00088 SAVE ALPHA, ONE, TWO 00089 DATA ALPHA/1.0D-4/, ONE/1.0D0/, TWO/2.0D0/ 00090 C 00091 KPRIN=IWM(LKPRIN) 00092 F1NRM = (FNRM*FNRM)/TWO 00093 RATIO = ONE 00094 C 00095 IF (KPRIN .GE. 2) THEN 00096 MSG = '------ IN ROUTINE DLINSK-- PNRM = (R1) )' 00097 CALL XERRWD(MSG, 40, 921, 0, 0, 0, 0, 1, PNRM, 0.0D0) 00098 ENDIF 00099 TAU = PNRM 00100 IVIO = 0 00101 RL = ONE 00102 C----------------------------------------------------------------------- 00103 C Check for violations of the constraints, if any are imposed. 00104 C If any violations are found, the step vector P is rescaled, and the 00105 C constraint check is repeated, until no violations are found. 00106 C----------------------------------------------------------------------- 00107 IF (ICNFLG .NE. 0) THEN 00108 10 CONTINUE 00109 CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW) 00110 CALL DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR) 00111 IF (IRET .EQ. 1) THEN 00112 IVIO = 1 00113 RATIO1 = TAU/PNRM 00114 RATIO = RATIO*RATIO1 00115 DO 20 I = 1,NEQ 00116 20 P(I) = P(I)*RATIO1 00117 PNRM = TAU 00118 IF (KPRIN .GE. 2) THEN 00119 MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)' 00120 CALL XERRWD(MSG, 50, 922, 0, 1, IVAR, 0, 1, PNRM, 0.0D0) 00121 ENDIF 00122 IF (PNRM .LE. STPTOL) THEN 00123 IRET = 1 00124 RETURN 00125 ENDIF 00126 GO TO 10 00127 ENDIF 00128 ENDIF 00129 C 00130 SLPI = (-TWO*F1NRM + RHOK*RHOK)*RATIO 00131 RLMIN = STPTOL/PNRM 00132 IF (LSOFF .EQ. 0 .AND. KPRIN .GE. 2) THEN 00133 MSG = '------ MIN. LAMBDA = (R1)' 00134 CALL XERRWD(MSG, 25, 923, 0, 0, 0, 0, 1, RLMIN, 0.0D0) 00135 ENDIF 00136 C----------------------------------------------------------------------- 00137 C Begin iteration to find RL value satisfying alpha-condition. 00138 C Update YNEW and YPNEW, then compute norm of new scaled residual and 00139 C perform alpha condition test. 00140 C----------------------------------------------------------------------- 00141 100 CONTINUE 00142 CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW) 00143 CALL DFNRMK (NEQ, YNEW, T, YPNEW, SAVR, R, CJ, WT, SQRTN, RSQRTN, 00144 * RES, IRES, PSOL, 0, IER, FNRMP, EPLIN, WP, IWP, PWK, RPAR, IPAR) 00145 IWM(LNRE) = IWM(LNRE) + 1 00146 IF (IRES .GE. 0) IWM(LNPS) = IWM(LNPS) + 1 00147 IF (IRES .NE. 0 .OR. IER .NE. 0) THEN 00148 IRET = 2 00149 RETURN 00150 ENDIF 00151 IF (LSOFF .EQ. 1) GO TO 150 00152 C 00153 F1NRMP = FNRMP*FNRMP/TWO 00154 IF (KPRIN .GE. 2) THEN 00155 MSG = '------ LAMBDA = (R1)' 00156 CALL XERRWD(MSG, 20, 924, 0, 0, 0, 0, 1, RL, 0.0D0) 00157 MSG = '------ NORM(F1) = (R1), NORM(F1NEW) = (R2)' 00158 CALL XERRWD(MSG, 43, 925, 0, 0, 0, 0, 2, F1NRM, F1NRMP) 00159 ENDIF 00160 IF (F1NRMP .GT. F1NRM + ALPHA*SLPI*RL) GO TO 200 00161 C----------------------------------------------------------------------- 00162 C Alpha-condition is satisfied, or linesearch is turned off. 00163 C Copy YNEW,YPNEW to Y,YPRIME and return. 00164 C----------------------------------------------------------------------- 00165 150 IRET = 0 00166 CALL DCOPY(NEQ, YNEW, 1, Y, 1) 00167 CALL DCOPY(NEQ, YPNEW, 1, YPRIME, 1) 00168 FNRM = FNRMP 00169 IF (KPRIN .GE. 1) THEN 00170 MSG = '------ LEAVING ROUTINE DLINSK, FNRM = (R1)' 00171 CALL XERRWD(MSG, 42, 926, 0, 0, 0, 0, 1, FNRM, 0.0D0) 00172 ENDIF 00173 RETURN 00174 C----------------------------------------------------------------------- 00175 C Alpha-condition not satisfied. Perform backtrack to compute new RL 00176 C value. If RL is less than RLMIN, i.e. no satisfactory YNEW,YPNEW can 00177 C be found sufficiently distinct from Y,YPRIME, then return IRET = 1. 00178 C----------------------------------------------------------------------- 00179 200 CONTINUE 00180 IF (RL .LT. RLMIN) THEN 00181 IRET = 1 00182 RETURN 00183 ENDIF 00184 C 00185 RL = RL/TWO 00186 GO TO 100 00187 C 00188 C----------------------- END OF SUBROUTINE DLINSK ---------------------- 00189 END
1.7.1
