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