Math-Lsoda
view release on metacpan or search on metacpan
C MITER = 3 means Newton iteration with Preconditioned
C Conjugate Gradient method (PCG)
C for the linear systems.
C MITER = 4 means Newton iteration with scaled Preconditioned
C Conjugate Gradient method (PCGS)
C for the linear systems.
C MITER = 9 means Newton iteration with only the
C user-supplied PSOL routine called (no Krylov
C iteration) for the linear systems.
C JPRE is ignored, and PSOL is called with LR = 0.
C See comments in the introduction about the choice of MITER.
C If MITER .ge. 1, the user must supply routines JAC and PSOL
C (the names are arbitrary) as described above.
C For MITER = 0, dummy arguments can be used.
C-----------------------------------------------------------------------
C Optional Inputs.
C
C The following is a list of the optional inputs provided for in the
C call sequence. (See also Part 2.) For each such input variable,
C this table lists its name as used in this documentation, its
C location in the call sequence, its meaning, and the default value.
C The use of any of these inputs requires IOPT = 1, and in that
C case all of these inputs are examined. A value of zero for any
C of these optional inputs will cause the default value to be used.
C Thus to use a subset of the optional inputs, simply preload
C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
C then set those of interest to nonzero values.
C
C Name Location Meaning and Default Value
C
C H0 RWORK(5) the step size to be attempted on the first step.
C The default value is determined by the solver.
C
C HMAX RWORK(6) the maximum absolute step size allowed.
C The default value is infinite.
C
C HMIN RWORK(7) the minimum absolute step size allowed.
C The default value is 0. (This lower bound is not
C enforced on the final step before reaching TCRIT
C when ITASK = 4 or 5.)
C
C DELT RWORK(8) convergence test constant in Krylov iteration
C algorithm. The default is .05.
C
C MAXORD IWORK(5) the maximum order to be allowed. The default
C value is 12 if METH = 1, and 5 if METH = 2.
C If MAXORD exceeds the default value, it will
C be reduced to the default value.
C If MAXORD is changed during the problem, it may
C cause the current order to be reduced.
C
C MXSTEP IWORK(6) maximum number of (internally defined) steps
C allowed during one call to the solver.
C The default value is 500.
C
C MXHNIL IWORK(7) maximum number of messages printed (per problem)
C warning that T + H = T on a step (H = step size).
C This must be positive to result in a non-default
C value. The default value is 10.
C
C MAXL IWORK(8) maximum number of iterations in the SPIOM, SPIGMR,
C PCG, or PCGS algorithm (.le. NEQ).
C The default is MAXL = MIN(5,NEQ).
C
C KMP IWORK(9) number of vectors on which orthogonalization
C is done in SPIOM or SPIGMR algorithm (.le. MAXL).
C The default is KMP = MAXL.
C Note: When KMP .lt. MAXL and MF = 22, the length
C of RWORK must be defined accordingly. See
C the definition of RWORK above.
C-----------------------------------------------------------------------
C Optional Outputs.
C
C As optional additional output from DLSODPK, the variables listed
C below are quantities related to the performance of DLSODPK
C which are available to the user. These are communicated by way of
C the work arrays, but also have internal mnemonic names as shown.
C Except where stated otherwise, all of these outputs are defined
C on any successful return from DLSODPK, and on any return with
C ISTATE = -1, -2, -4, -5, -6, or -7. On an illegal input return
C (ISTATE = -3), they will be unchanged from their existing values
C (if any), except possibly for TOLSF, LENRW, and LENIW.
C On any error return, outputs relevant to the error will be defined,
C as noted below.
C
C Name Location Meaning
C
C HU RWORK(11) the step size in t last used (successfully).
C
C HCUR RWORK(12) the step size to be attempted on the next step.
C
C TCUR RWORK(13) the current value of the independent variable
C which the solver has actually reached, i.e. the
C current internal mesh point in t. On output, TCUR
C will always be at least as far as the argument
C T, but may be farther (if interpolation was done).
C
C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
C computed when a request for too much accuracy was
C detected (ISTATE = -3 if detected at the start of
C the problem, ISTATE = -2 otherwise). If ITOL is
C left unaltered but RTOL and ATOL are uniformly
C scaled up by a factor of TOLSF for the next call,
C then the solver is deemed likely to succeed.
C (The user may also ignore TOLSF and alter the
C tolerance parameters in any other way appropriate.)
C
C NST IWORK(11) the number of steps taken for the problem so far.
C
C NFE IWORK(12) the number of f evaluations for the problem so far.
C
C NPE IWORK(13) the number of calls to JAC so far (for Jacobian
C evaluation associated with preconditioning).
C
C NQU IWORK(14) the method order last used (successfully).
C
C NQCUR IWORK(15) the order to be attempted on the next step.
C
C IMXER IWORK(16) the index of the component of largest magnitude in
C the weighted local error vector ( E(i)/EWT(i) ),
C on an error return with ISTATE = -4 or -5.
C
C LENRW IWORK(17) the length of RWORK actually required.
C This is defined on normal returns and on an illegal
C input return for insufficient storage.
C
C LENIW IWORK(18) the length of IWORK actually required.
C This is defined on normal returns and on an illegal
C input return for insufficient storage.
C
C NNI IWORK(19) number of nonlinear iterations so far (each of
C which calls an iterative linear solver).
C
C NLI IWORK(20) number of linear iterations so far.
C Note: A measure of the success of algorithm is
C the average number of linear iterations per
C nonlinear iteration, given by NLI/NNI.
C If this is close to MAXL, MAXL may be too small.
C
C NPS IWORK(21) number of preconditioning solve operations
C (PSOL calls) so far.
C
C NCFN IWORK(22) number of convergence failures of the nonlinear
C (Newton) iteration so far.
C Note: A measure of success is the overall
C rate of nonlinear convergence failures, NCFN/NST.
C
C NCFL IWORK(23) number of convergence failures of the linear
C iteration so far.
C Note: A measure of success is the overall
C rate of linear convergence failures, NCFL/NNI.
C
C The following two arrays are segments of the RWORK array which
C may also be of interest to the user as optional outputs.
C For each array, the table below gives its internal name,
C its base address in RWORK, and its description.
C
C Name Base Address Description
C
C YH 21 the Nordsieck history array, of size NYH by
C (NQCUR + 1), where NYH is the initial value
C of NEQ. For j = 0,1,...,NQCUR, column j+1
C of YH contains HCUR**j/factorial(j) times
C the j-th derivative of the interpolating
C polynomial currently representing the solution,
C evaluated at t = TCUR.
C
C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
C corrections on each step, scaled on output
C to represent the estimated local error in y
C on the last step. This is the vector E in
C the description of the error control. It is
C defined only on a successful return from
C DLSODPK.
C
C-----------------------------------------------------------------------
C Part 2. Other Routines Callable.
C
C The following are optional calls which the user may make to
C gain additional capabilities in conjunction with DLSODPK.
C (The routines XSETUN and XSETF are designed to conform to the
C SLATEC error handling package.)
C
C Form of Call Function
C CALL XSETUN(LUN) Set the logical unit number, LUN, for
C output of messages from DLSODPK, if
C the default is not desired.
C The default value of lun is 6.
C
C CALL XSETF(MFLAG) Set a flag to control the printing of
C messages by DLSODPK.
C MFLAG = 0 means do not print. (Danger:
C This risks losing valuable information.)
C MFLAG = 1 means print (the default).
C
C Either of the above calls may be made at
NWARN = 0
GO TO (210, 250, 220, 230, 240), ITASK
210 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
IF (IFLAG .NE. 0) GO TO 627
T = TOUT
GO TO 420
220 TP = TN - HU*(1.0D0 + 100.0D0*UROUND)
IF ((TP - TOUT)*H .GT. 0.0D0) GO TO 623
IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
GO TO 400
230 TCRIT = RWORK(1)
IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
IF ((TCRIT - TOUT)*H .LT. 0.0D0) GO TO 625
IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 245
CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
IF (IFLAG .NE. 0) GO TO 627
T = TOUT
GO TO 420
240 TCRIT = RWORK(1)
IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
245 HMX = ABS(TN) + ABS(H)
IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
IF (IHIT) GO TO 400
TNEXT = TN + H*(1.0D0 + 4.0D0*UROUND)
IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250
H = (TCRIT - TN)*(1.0D0 - 4.0D0*UROUND)
IF (ISTATE .EQ. 2) JSTART = -2
C-----------------------------------------------------------------------
C Block E.
C The next block is normally executed for all calls and contains
C the call to the one-step core integrator DSTODPK.
C
C This is a looping point for the integration steps.
C
C First check for too many steps being taken,
C Check for poor Newton/Krylov method performance, update EWT (if not
C at start of problem), check for too much accuracy being requested,
C and check for H below the roundoff level in T.
C-----------------------------------------------------------------------
250 CONTINUE
IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
NSTD = NST - NSLAST
NNID = NNI - NNI0
IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 255
AVDIM = REAL(NLI - NLI0)/REAL(NNID)
RCFN = REAL(NCFN - NCFN0)/REAL(NSTD)
RCFL = REAL(NCFL - NCFL0)/REAL(NNID)
LAVD = AVDIM .GT. (MAXL - 0.05D0)
LCFN = RCFN .GT. 0.9D0
LCFL = RCFL .GT. 0.9D0
LWARN = LAVD .OR. LCFN .OR. LCFL
IF (.NOT.LWARN) GO TO 255
NWARN = NWARN + 1
IF (NWARN .GT. 10) GO TO 255
IF (LAVD) THEN
MSG='DLSODPK- Warning. Poor iterative algorithm performance seen '
CALL XERRWD (MSG, 60, 111, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
ENDIF
IF (LAVD) THEN
MSG=' at T = R1 by average no. of linear iterations = R2 '
CALL XERRWD (MSG, 60, 111, 0, 0, 0, 0, 2, TN, AVDIM)
ENDIF
IF (LCFN) THEN
MSG='DLSODPK- Warning. Poor iterative algorithm performance seen '
CALL XERRWD (MSG, 60, 112, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
ENDIF
IF (LCFN) THEN
MSG=' at T = R1 by nonlinear convergence failure rate = R2 '
CALL XERRWD (MSG, 60, 112, 0, 0, 0, 0, 2, TN, RCFN)
ENDIF
IF (LCFL) THEN
MSG='DLSODPK- Warning. Poor iterative algorithm performance seen '
CALL XERRWD (MSG, 60, 113, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
ENDIF
IF (LCFL) THEN
MSG=' at T = R1 by linear convergence failure rate = R2 '
CALL XERRWD (MSG, 60, 113, 0, 0, 0, 0, 2, TN, RCFL)
ENDIF
255 CONTINUE
CALL DEWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
DO 260 I = 1,N
IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 510
260 RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1)
270 TOLSF = UROUND*DVNORM (N, RWORK(LYH), RWORK(LEWT))
IF (TOLSF .LE. 1.0D0) GO TO 280
TOLSF = TOLSF*2.0D0
IF (NST .EQ. 0) GO TO 626
GO TO 520
280 IF ((TN + H) .NE. TN) GO TO 290
NHNIL = NHNIL + 1
IF (NHNIL .GT. MXHNIL) GO TO 290
MSG = 'DLSODPK- Warning..Internal T(=R1) and H(=R2) are '
CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG=' such that in the machine, T + H = T on the next step '
CALL XERRWD (MSG, 60, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' (H = step size). Solver will continue anyway.'
CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 2, TN, H)
IF (NHNIL .LT. MXHNIL) GO TO 290
MSG = 'DLSODPK- Above warning has been issued I1 times. '
CALL XERRWD (MSG, 50, 102, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' It will not be issued again for this problem.'
CALL XERRWD (MSG, 50, 102, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0)
290 CONTINUE
C-----------------------------------------------------------------------
C CALL DSTODPK(NEQ,Y,YH,NYH,YH,EWT,SAVF,SAVX,ACOR,WM,IWM,F,JAC,PSOL)
C-----------------------------------------------------------------------
CALL DSTODPK (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT),
1 RWORK(LSAVF), RWORK(LSAVX), RWORK(LACOR), RWORK(LWM),
2 IWORK(LIWM), F, JAC, PSOL)
KGO = 1 - KFLAG
GO TO (300, 530, 540, 550), KGO
C-----------------------------------------------------------------------
C Block F.
C The following block handles the case of a successful return from the
C core integrator (KFLAG = 0). Test for stop conditions.
C-----------------------------------------------------------------------
300 INIT = 1
GO TO (310, 400, 330, 340, 350), ITASK
C ITASK = 1. If TOUT has been reached, interpolate. -------------------
310 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
C to a point slightly past that root, so that DLSODKR can
C locate the root by interpolation.
C
C Subroutine G may access user-defined quantities in
C NEQ(2),... and Y(NEQ(1)+1),... if NEQ is an array
C (dimensioned in G) and/or Y has length exceeding NEQ(1).
C See the descriptions of NEQ and Y above.
C
C NG = number of constraint functions g(i). If there are none,
C set NG = 0, and pass a dummy name for G.
C
C JROOT = integer array of length NG. Used only for output.
C On a return with ISTATE = 3 (one or more roots found),
C JROOT(i) = 1 if g(i) has a root at t, or JROOT(i) = 0 if not.
C-----------------------------------------------------------------------
C Optional Inputs.
C
C The following is a list of the optional inputs provided for in the
C call sequence. (See also Part 2.) For each such input variable,
C this table lists its name as used in this documentation, its
C location in the call sequence, its meaning, and the default value.
C The use of any of these inputs requires IOPT = 1, and in that
C case all of these inputs are examined. A value of zero for any
C of these optional inputs will cause the default value to be used.
C Thus to use a subset of the optional inputs, simply preload
C locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and
C then set those of interest to nonzero values.
C
C Name Location Meaning and Default Value
C
C H0 RWORK(5) the step size to be attempted on the first step.
C The default value is determined by the solver.
C
C HMAX RWORK(6) the maximum absolute step size allowed.
C The default value is infinite.
C
C HMIN RWORK(7) the minimum absolute step size allowed.
C The default value is 0. (This lower bound is not
C enforced on the final step before reaching TCRIT
C when ITASK = 4 or 5.)
C
C DELT RWORK(8) convergence test constant in Krylov iteration
C algorithm. The default is .05.
C
C MAXORD IWORK(5) the maximum order to be allowed. The default
C value is 12 if METH = 1, and 5 if METH = 2.
C If MAXORD exceeds the default value, it will
C be reduced to the default value.
C If MAXORD is changed during the problem, it may
C cause the current order to be reduced.
C
C MXSTEP IWORK(6) maximum number of (internally defined) steps
C allowed during one call to the solver.
C The default value is 500.
C
C MXHNIL IWORK(7) maximum number of messages printed (per problem)
C warning that T + H = T on a step (H = step size).
C This must be positive to result in a non-default
C value. The default value is 10.
C
C MAXL IWORK(8) maximum number of iterations in the SPIOM, SPIGMR,
C PCG, or PCGS algorithm (.le. NEQ).
C The default is MAXL = MIN(5,NEQ).
C
C KMP IWORK(9) number of vectors on which orthogonalization
C is done in SPIOM or SPIGMR algorithm (.le. MAXL).
C The default is KMP = MAXL.
C Note: When KMP .lt. MAXL and MF = 22, the length
C of RWORK must be defined accordingly. See
C the definition of RWORK above.
C-----------------------------------------------------------------------
C Optional Outputs.
C
C As optional additional output from DLSODKR, the variables listed
C below are quantities related to the performance of DLSODKR
C which are available to the user. These are communicated by way of
C the work arrays, but also have internal mnemonic names as shown.
C Except where stated otherwise, all of these outputs are defined
C on any successful return from DLSODKR, and on any return with
C ISTATE = -1, -2, -4, -5, -6, or -7. On an illegal input return
C (ISTATE = -3), they will be unchanged from their existing values
C (if any), except possibly for TOLSF, LENRW, and LENIW.
C On any error return, outputs relevant to the error will be defined,
C as noted below.
C
C Name Location Meaning
C
C HU RWORK(11) the step size in t last used (successfully).
C
C HCUR RWORK(12) the step size to be attempted on the next step.
C
C TCUR RWORK(13) the current value of the independent variable
C which the solver has actually reached, i.e. the
C current internal mesh point in t. On output, TCUR
C will always be at least as far as the argument
C T, but may be farther (if interpolation was done).
C
C TOLSF RWORK(14) a tolerance scale factor, greater than 1.0,
C computed when a request for too much accuracy was
C detected (ISTATE = -3 if detected at the start of
C the problem, ISTATE = -2 otherwise). If ITOL is
C left unaltered but RTOL and ATOL are uniformly
C scaled up by a factor of TOLSF for the next call,
C then the solver is deemed likely to succeed.
C (The user may also ignore TOLSF and alter the
C tolerance parameters in any other way appropriate.)
C
C NGE IWORK(10) the number of g evaluations for the problem so far.
C
C NST IWORK(11) the number of steps taken for the problem so far.
C
C NFE IWORK(12) the number of f evaluations for the problem so far.
C
C NPE IWORK(13) the number of calls to JAC so far (for evaluation
C of preconditioners).
C
C NQU IWORK(14) the method order last used (successfully).
C
C NQCUR IWORK(15) the order to be attempted on the next step.
C
C IMXER IWORK(16) the index of the component of largest magnitude in
C the weighted local error vector ( E(i)/EWT(i) ),
C on an error return with ISTATE = -4 or -5.
C
C LENRW IWORK(17) the length of RWORK actually required.
C This is defined on normal returns and on an illegal
C input return for insufficient storage.
C
C LENIW IWORK(18) the length of IWORK actually required.
C This is defined on normal returns and on an illegal
C input return for insufficient storage.
C
C NNI IWORK(19) number of nonlinear iterations so far (each of
C which calls an iterative linear solver).
C
C NLI IWORK(20) number of linear iterations so far.
C Note: A measure of the success of algorithm is
C the average number of linear iterations per
C nonlinear iteration, given by NLI/NNI.
C If this is close to MAXL, MAXL may be too small.
C
C NPS IWORK(21) number of preconditioning solve operations
C (PSOL calls) so far.
C
C NCFN IWORK(22) number of convergence failures of the nonlinear
C (Newton) iteration so far.
C Note: A measure of success is the overall
C rate of nonlinear convergence failures, NCFN/NST.
C
C NCFL IWORK(23) number of convergence failures of the linear
C iteration so far.
C Note: A measure of success is the overall
C rate of linear convergence failures, NCFL/NNI.
C
C NSFI IWORK(24) number of functional iteration steps so far.
C Note: A measure of the extent to which the
C problem is nonstiff is the ratio NSFI/NST.
C
C NJEV IWORK(25) number of JAC calls with JOK = -1 so far
C (number of evaluations of Jacobian data).
C
C The following two arrays are segments of the RWORK array which
C may also be of interest to the user as optional outputs.
C For each array, the table below gives its internal name,
C its base address in RWORK, and its description.
C
C Name Base Address Description
C
C YH 21 + 3*NG the Nordsieck history array, of size NYH by
C (NQCUR + 1), where NYH is the initial value
C of NEQ. For j = 0,1,...,NQCUR, column j+1
C of YH contains HCUR**j/factorial(j) times
C the j-th derivative of the interpolating
C polynomial currently representing the solution,
C evaluated at t = TCUR.
C
C ACOR LENRW-NEQ+1 array of size NEQ used for the accumulated
C corrections on each step, scaled on output
C to represent the estimated local error in y
C on the last step. This is the vector E in
C the description of the error control. It is
C defined only on a successful return from
C DLSODKR.
C
C-----------------------------------------------------------------------
C Part 2. Other Routines Callable.
C
C The following are optional calls which the user may make to
C gain additional capabilities in conjunction with DLSODKR.
C (The routines XSETUN and XSETF are designed to conform to the
C SLATEC error handling package.)
C
C Form of Call Function
C CALL XSETUN(LUN) Set the logical unit number, LUN, for
C output of messages from DLSODKR, if
C the default is not desired.
C The default value of LUN is 6.
C
210 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
IF (IFLAG .NE. 0) GO TO 627
T = TOUT
GO TO 420
220 TP = TN - HU*(1.0D0 + 100.0D0*UROUND)
IF ((TP - TOUT)*H .GT. 0.0D0) GO TO 623
IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
GO TO 400
230 TCRIT = RWORK(1)
IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
IF ((TCRIT - TOUT)*H .LT. 0.0D0) GO TO 625
IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 245
CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
IF (IFLAG .NE. 0) GO TO 627
T = TOUT
GO TO 420
240 TCRIT = RWORK(1)
IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
245 HMX = ABS(TN) + ABS(H)
IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
IF (IHIT) T = TCRIT
IF (IRFP .EQ. 1 .AND. TLAST .NE. TN .AND. ITASK .EQ. 5) GO TO 400
IF (IHIT) GO TO 400
TNEXT = TN + H*(1.0D0 + 4.0D0*UROUND)
IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250
H = (TCRIT - TN)*(1.0D0 - 4.0D0*UROUND)
IF (ISTATE .EQ. 2) JSTART = -2
C-----------------------------------------------------------------------
C Block E.
C The next block is normally executed for all calls and contains
C the call to the one-step core integrator DSTOKA.
C
C This is a looping point for the integration steps.
C
C First check for too many steps being taken,
C check for poor Newton/Krylov method performance, update EWT (if not
C at start of problem), check for too much accuracy being requested,
C and check for H below the roundoff level in T.
C-----------------------------------------------------------------------
250 CONTINUE
IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
NSTD = NST - NSLAST
NNID = NNI - NNI0
IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 255
AVDIM = REAL(NLI - NLI0)/REAL(NNID)
RCFN = REAL(NCFN - NCFN0)/REAL(NSTD)
RCFL = REAL(NCFL - NCFL0)/REAL(NNID)
LAVD = AVDIM .GT. (MAXL - 0.05D0)
LCFN = RCFN .GT. 0.9D0
LCFL = RCFL .GT. 0.9D0
LWARN = LAVD .OR. LCFN .OR. LCFL
IF (.NOT.LWARN) GO TO 255
NWARN = NWARN + 1
IF (NWARN .GT. 10) GO TO 255
IF (LAVD) THEN
MSG='DLSODKR- Warning. Poor iterative algorithm performance seen '
CALL XERRWD (MSG, 60, 111, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
ENDIF
IF (LAVD) THEN
MSG=' at T = R1 by average no. of linear iterations = R2 '
CALL XERRWD (MSG, 60, 111, 0, 0, 0, 0, 2, TN, AVDIM)
ENDIF
IF (LCFN) THEN
MSG='DLSODKR- Warning. Poor iterative algorithm performance seen '
CALL XERRWD (MSG, 60, 112, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
ENDIF
IF (LCFN) THEN
MSG=' at T = R1 by nonlinear convergence failure rate = R2 '
CALL XERRWD (MSG, 60, 112, 0, 0, 0, 0, 2, TN, RCFN)
ENDIF
IF (LCFL) THEN
MSG='DLSODKR- Warning. Poor iterative algorithm performance seen '
CALL XERRWD (MSG, 60, 113, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
ENDIF
IF (LCFL) THEN
MSG=' at T = R1 by linear convergence failure rate = R2 '
CALL XERRWD (MSG, 60, 113, 0, 0, 0, 0, 2, TN, RCFL)
ENDIF
255 CONTINUE
CALL DEWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
DO 260 I = 1,N
IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 510
260 RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1)
270 TOLSF = UROUND*DVNORM (N, RWORK(LYH), RWORK(LEWT))
IF (TOLSF .LE. 1.0D0) GO TO 280
TOLSF = TOLSF*2.0D0
IF (NST .EQ. 0) GO TO 626
GO TO 520
280 IF ((TN + H) .NE. TN) GO TO 290
NHNIL = NHNIL + 1
IF (NHNIL .GT. MXHNIL) GO TO 290
MSG = 'DLSODKR- Warning.. Internal T(=R1) and H(=R2) are'
CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG=' such that in the machine, T + H = T on the next step '
CALL XERRWD (MSG, 60, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' (H = step size). Solver will continue anyway.'
CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 2, TN, H)
IF (NHNIL .LT. MXHNIL) GO TO 290
MSG = 'DLSODKR- Above warning has been issued I1 times. '
CALL XERRWD (MSG, 50, 102, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' It will not be issued again for this problem.'
CALL XERRWD (MSG, 50, 102, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0)
290 CONTINUE
C-----------------------------------------------------------------------
C CALL DSTOKA(NEQ,Y,YH,NYH,YH,EWT,SAVF,SAVX,ACOR,WM,IWM,F,JAC,PSOL)
C-----------------------------------------------------------------------
CALL DSTOKA (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT),
1 RWORK(LSAVF), RWORK(LSAVX), RWORK(LACOR), RWORK(LWM),
2 IWORK(LIWM), F, JAC, PSOL)
KGO = 1 - KFLAG
GO TO (300, 530, 540, 550), KGO
C-----------------------------------------------------------------------
C Block F.
C The following block handles the case of a successful return from the
C core integrator (KFLAG = 0).
C Call DRCHEK to check for a root within the last step.
C Then, if no root was found, check for stop conditions.
C-----------------------------------------------------------------------
300 INIT = 1
C
( run in 0.839 second using v1.01-cache-2.11-cpan-71847e10f99 )