octave-jordi: liboctave/cruft/odepack/sstode.f annotate

annotate liboctave/cruft/odepack/sstode.f @ 21003:f7e416862e90 stable

doc: fix spelling of "occurred". * etc/OLD-ChangeLogs/scripts-ChangeLog, libinterp/corefcn/error.cc, libinterp/corefcn/urlwrite.cc, libinterp/octave-value/ov-oncleanup.cc, liboctave/cruft/dassl/ddastp.f, liboctave/cruft/odepack/sstode.f, liboctave/cruft/odepack/stode.f, liboctave/util/oct-inttypes.cc, scripts/gui/errordlg.m, scripts/gui/warndlg.m: fix spelling of "occurred"

author	Rafael Laboissiere <rafael@laboissiere.net>
date	Tue, 29 Dec 2015 11:51:16 +0100
parents	648dabbb4c6b
children

rev	line source
7793 96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	1 SUBROUTINE SSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	2 1 WM, IWM, F, JAC, PJAC, SLVS)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	3 C***BEGIN PROLOGUE SSTODE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	4 C***SUBSIDIARY
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	5 C***PURPOSE Performs one step of an ODEPACK integration.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	6 C***TYPE SINGLE PRECISION (SSTODE-S, DSTODE-D)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	7 C***AUTHOR Hindmarsh, Alan C., (LLNL)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	8 C***DESCRIPTION
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	9 C
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	10 C SSTODE performs one step of the integration of an initial value
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	11 C problem for a system of ordinary differential equations.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	12 C Note: SSTODE is independent of the value of the iteration method
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	13 C indicator MITER, when this is .ne. 0, and hence is independent
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	14 C of the type of chord method used, or the Jacobian structure.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	15 C Communication with SSTODE is done with the following variables:
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	16 C
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	17 C NEQ = integer array containing problem size in NEQ(1), and
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	18 C passed as the NEQ argument in all calls to F and JAC.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	19 C Y = an array of length .ge. N used as the Y argument in
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	20 C all calls to F and JAC.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	21 C YH = an NYH by LMAX array containing the dependent variables
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	22 C and their approximate scaled derivatives, where
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	23 C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	24 C j-th derivative of y(i), scaled by h**j/factorial(j)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	25 C (j = 0,1,...,NQ). on entry for the first step, the first
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	26 C two columns of YH must be set from the initial values.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	27 C NYH = a constant integer .ge. N, the first dimension of YH.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	28 C YH1 = a one-dimensional array occupying the same space as YH.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	29 C EWT = an array of length N containing multiplicative weights
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	30 C for local error measurements. Local errors in Y(i) are
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	31 C compared to 1.0/EWT(i) in various error tests.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	32 C SAVF = an array of working storage, of length N.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	33 C Also used for input of YH(*,MAXORD+2) when JSTART = -1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	34 C and MAXORD .lt. the current order NQ.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	35 C ACOR = a work array of length N, used for the accumulated
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	36 C corrections. On a successful return, ACOR(i) contains
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	37 C the estimated one-step local error in Y(i).
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	38 C WM,IWM = real and integer work arrays associated with matrix
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	39 C operations in chord iteration (MITER .ne. 0).
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	40 C PJAC = name of routine to evaluate and preprocess Jacobian matrix
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	41 C and P = I - hel0JAC, if a chord method is being used.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	42 C SLVS = name of routine to solve linear system in chord iteration.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	43 C CCMAX = maximum relative change in h*el0 before PJAC is called.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	44 C H = the step size to be attempted on the next step.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	45 C H is altered by the error control algorithm during the
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	46 C problem. H can be either positive or negative, but its
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	47 C sign must remain constant throughout the problem.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	48 C HMIN = the minimum absolute value of the step size h to be used.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	49 C HMXI = inverse of the maximum absolute value of h to be used.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	50 C HMXI = 0.0 is allowed and corresponds to an infinite hmax.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	51 C HMIN and HMXI may be changed at any time, but will not
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	52 C take effect until the next change of h is considered.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	53 C TN = the independent variable. TN is updated on each step taken.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	54 C JSTART = an integer used for input only, with the following
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	55 C values and meanings:
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	56 C 0 perform the first step.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	57 C .gt.0 take a new step continuing from the last.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	58 C -1 take the next step with a new value of H, MAXORD,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	59 C N, METH, MITER, and/or matrix parameters.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	60 C -2 take the next step with a new value of H,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	61 C but with other inputs unchanged.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	62 C On return, JSTART is set to 1 to facilitate continuation.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	63 C KFLAG = a completion code with the following meanings:
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	64 C 0 the step was succesful.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	65 C -1 the requested error could not be achieved.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	66 C -2 corrector convergence could not be achieved.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	67 C -3 fatal error in PJAC or SLVS.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	68 C A return with KFLAG = -1 or -2 means either
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	69 C abs(H) = HMIN or 10 consecutive failures occurred.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	70 C On a return with KFLAG negative, the values of TN and
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	71 C the YH array are as of the beginning of the last
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	72 C step, and H is the last step size attempted.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	73 C MAXORD = the maximum order of integration method to be allowed.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	74 C MAXCOR = the maximum number of corrector iterations allowed.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	75 C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	76 C MXNCF = maximum number of convergence failures allowed.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	77 C METH/MITER = the method flags. See description in driver.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	78 C N = the number of first-order differential equations.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	79 C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	80 C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	81 C
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	82 C***SEE ALSO SLSODE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	83 C***ROUTINES CALLED SCFODE, SVNORM
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	84 C***COMMON BLOCKS SLS001
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	85 C***REVISION HISTORY (YYMMDD)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	86 C 791129 DATE WRITTEN
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	87 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	88 C 890503 Minor cosmetic changes. (FNF)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	89 C 930809 Renamed to allow single/double precision versions. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	90 C 010413 Reduced size of Common block /SLS001/. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	91 C 031105 Restored 'own' variables to Common block /SLS001/, to
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	92 C enable interrupt/restart feature. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	93 C***END PROLOGUE SSTODE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	94 C**End
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	95 EXTERNAL F, JAC, PJAC, SLVS
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	96 INTEGER NEQ, NYH, IWM
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	97 REAL Y, YH, YH1, EWT, SAVF, ACOR, WM
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	98 DIMENSION NEQ(), Y(), YH(NYH,), YH1(), EWT(), SAVF(),
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	99 1 ACOR(), WM(), IWM(*)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	100 INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	101 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	102 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	103 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	104 INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	105 REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	106 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	107 REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	108 1 R, RH, RHDN, RHSM, RHUP, TOLD, SVNORM
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	109 COMMON /SLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	110 1 HOLD, RMAX, TESCO(3,12),
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	111 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	112 3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	113 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	114 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	115 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	116 C
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	117 C***FIRST EXECUTABLE STATEMENT SSTODE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	118 KFLAG = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	119 TOLD = TN
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	120 NCF = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	121 IERPJ = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	122 IERSL = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	123 JCUR = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	124 ICF = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	125 DELP = 0.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	126 IF (JSTART .GT. 0) GO TO 200
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	127 IF (JSTART .EQ. -1) GO TO 100
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	128 IF (JSTART .EQ. -2) GO TO 160
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	129 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	130 C On the first call, the order is set to 1, and other variables are
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	131 C initialized. RMAX is the maximum ratio by which H can be increased
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	132 C in a single step. It is initially 1.E4 to compensate for the small
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	133 C initial H, but then is normally equal to 10. If a failure
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	134 C occurs (in corrector convergence or error test), RMAX is set to 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	135 C for the next increase.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	136 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	137 LMAX = MAXORD + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	138 NQ = 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	139 L = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	140 IALTH = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	141 RMAX = 10000.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	142 RC = 0.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	143 EL0 = 1.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	144 CRATE = 0.7E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	145 HOLD = H
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	146 MEO = METH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	147 NSLP = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	148 IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	149 IRET = 3
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	150 GO TO 140
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	151 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	152 C The following block handles preliminaries needed when JSTART = -1.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	153 C IPUP is set to MITER to force a matrix update.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	154 C If an order increase is about to be considered (IALTH = 1),
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	155 C IALTH is reset to 2 to postpone consideration one more step.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	156 C If the caller has changed METH, SCFODE is called to reset
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	157 C the coefficients of the method.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	158 C If the caller has changed MAXORD to a value less than the current
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	159 C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	160 C If H is to be changed, YH must be rescaled.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	161 C If H or METH is being changed, IALTH is reset to L = NQ + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	162 C to prevent further changes in H for that many steps.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	163 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	164 100 IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	165 LMAX = MAXORD + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	166 IF (IALTH .EQ. 1) IALTH = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	167 IF (METH .EQ. MEO) GO TO 110
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	168 CALL SCFODE (METH, ELCO, TESCO)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	169 MEO = METH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	170 IF (NQ .GT. MAXORD) GO TO 120
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	171 IALTH = L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	172 IRET = 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	173 GO TO 150
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	174 110 IF (NQ .LE. MAXORD) GO TO 160
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	175 120 NQ = MAXORD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	176 L = LMAX
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	177 DO 125 I = 1,L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	178 125 EL(I) = ELCO(I,NQ)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	179 NQNYH = NQ*NYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	180 RC = RC*EL(1)/EL0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	181 EL0 = EL(1)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	182 CONIT = 0.5E0/(NQ+2)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	183 DDN = SVNORM (N, SAVF, EWT)/TESCO(1,L)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	184 EXDN = 1.0E0/L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	185 RHDN = 1.0E0/(1.3E0DDN*EXDN + 0.0000013E0)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	186 RH = MIN(RHDN,1.0E0)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	187 IREDO = 3
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	188 IF (H .EQ. HOLD) GO TO 170
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	189 RH = MIN(RH,ABS(H/HOLD))
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	190 H = HOLD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	191 GO TO 175
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	192 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	193 C SCFODE is called to get all the integration coefficients for the
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	194 C current METH. Then the EL vector and related constants are reset
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	195 C whenever the order NQ is changed, or at the start of the problem.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	196 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	197 140 CALL SCFODE (METH, ELCO, TESCO)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	198 150 DO 155 I = 1,L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	199 155 EL(I) = ELCO(I,NQ)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	200 NQNYH = NQ*NYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	201 RC = RC*EL(1)/EL0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	202 EL0 = EL(1)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	203 CONIT = 0.5E0/(NQ+2)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	204 GO TO (160, 170, 200), IRET
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	205 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	206 C If H is being changed, the H ratio RH is checked against
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	207 C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	208 C L = NQ + 1 to prevent a change of H for that many steps, unless
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	209 C forced by a convergence or error test failure.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	210 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	211 160 IF (H .EQ. HOLD) GO TO 200
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	212 RH = H/HOLD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	213 H = HOLD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	214 IREDO = 3
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	215 GO TO 175
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	216 170 RH = MAX(RH,HMIN/ABS(H))
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	217 175 RH = MIN(RH,RMAX)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	218 RH = RH/MAX(1.0E0,ABS(H)HMXIRH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	219 R = 1.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	220 DO 180 J = 2,L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	221 R = R*RH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	222 DO 180 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	223 180 YH(I,J) = YH(I,J)*R
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	224 H = H*RH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	225 RC = RC*RH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	226 IALTH = L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	227 IF (IREDO .EQ. 0) GO TO 690
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	228 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	229 C This section computes the predicted values by effectively
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	230 C multiplying the YH array by the Pascal Triangle matrix.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	231 C RC is the ratio of new to old values of the coefficient H*EL(1).
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	232 C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	233 C to force PJAC to be called, if a Jacobian is involved.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	234 C In any case, PJAC is called at least every MSBP steps.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	235 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	236 200 IF (ABS(RC-1.0E0) .GT. CCMAX) IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	237 IF (NST .GE. NSLP+MSBP) IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	238 TN = TN + H
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	239 I1 = NQNYH + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	240 DO 215 JB = 1,NQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	241 I1 = I1 - NYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	242 Cdir$ ivdep
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	243 DO 210 I = I1,NQNYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	244 210 YH1(I) = YH1(I) + YH1(I+NYH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	245 215 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	246 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	247 C Up to MAXCOR corrector iterations are taken. A convergence test is
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	248 C made on the R.M.S. norm of each correction, weighted by the error
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	249 C weight vector EWT. The sum of the corrections is accumulated in the
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	250 C vector ACOR(i). The YH array is not altered in the corrector loop.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	251 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	252 220 M = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	253 DO 230 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	254 230 Y(I) = YH(I,1)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	255 CALL F (NEQ, TN, Y, SAVF)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	256 NFE = NFE + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	257 IF (IPUP .LE. 0) GO TO 250
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	258 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	259 C If indicated, the matrix P = I - hel(1)J is reevaluated and
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	260 C preprocessed before starting the corrector iteration. IPUP is set
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	261 C to 0 as an indicator that this has been done.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	262 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	263 CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	264 IPUP = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	265 RC = 1.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	266 NSLP = NST
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	267 CRATE = 0.7E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	268 IF (IERPJ .NE. 0) GO TO 430
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	269 250 DO 260 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	270 260 ACOR(I) = 0.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	271 270 IF (MITER .NE. 0) GO TO 350
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	272 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	273 C In the case of functional iteration, update Y directly from
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	274 C the result of the last function evaluation.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	275 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	276 DO 290 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	277 SAVF(I) = H*SAVF(I) - YH(I,2)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	278 290 Y(I) = SAVF(I) - ACOR(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	279 DEL = SVNORM (N, Y, EWT)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	280 DO 300 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	281 Y(I) = YH(I,1) + EL(1)*SAVF(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	282 300 ACOR(I) = SAVF(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	283 GO TO 400
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	284 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	285 C In the case of the chord method, compute the corrector error,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	286 C and solve the linear system with that as right-hand side and
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	287 C P as coefficient matrix.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	288 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	289 350 DO 360 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	290 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	291 CALL SLVS (WM, IWM, Y, SAVF)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	292 IF (IERSL .LT. 0) GO TO 430
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	293 IF (IERSL .GT. 0) GO TO 410
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	294 DEL = SVNORM (N, Y, EWT)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	295 DO 380 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	296 ACOR(I) = ACOR(I) + Y(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	297 380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	298 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	299 C Test for convergence. If M.gt.0, an estimate of the convergence
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	300 C rate constant is stored in CRATE, and this is used in the test.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	301 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	302 400 IF (M .NE. 0) CRATE = MAX(0.2E0*CRATE,DEL/DELP)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	303 DCON = DELMIN(1.0E0,1.5E0CRATE)/(TESCO(2,NQ)*CONIT)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	304 IF (DCON .LE. 1.0E0) GO TO 450
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	305 M = M + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	306 IF (M .EQ. MAXCOR) GO TO 410
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	307 IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	308 DELP = DEL
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	309 CALL F (NEQ, TN, Y, SAVF)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	310 NFE = NFE + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	311 GO TO 270
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	312 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	313 C The corrector iteration failed to converge.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	314 C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	315 C the next try. Otherwise the YH array is retracted to its values
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	316 C before prediction, and H is reduced, if possible. If H cannot be
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	317 C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	318 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	319 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	320 ICF = 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	321 IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	322 GO TO 220
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	323 430 ICF = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	324 NCF = NCF + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	325 RMAX = 2.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	326 TN = TOLD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	327 I1 = NQNYH + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	328 DO 445 JB = 1,NQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	329 I1 = I1 - NYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	330 Cdir$ ivdep
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	331 DO 440 I = I1,NQNYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	332 440 YH1(I) = YH1(I) - YH1(I+NYH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	333 445 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	334 IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	335 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	336 IF (NCF .EQ. MXNCF) GO TO 670
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	337 RH = 0.25E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	338 IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	339 IREDO = 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	340 GO TO 170
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	341 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	342 C The corrector has converged. JCUR is set to 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	343 C to signal that the Jacobian involved may need updating later.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	344 C The local error test is made and control passes to statement 500
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	345 C if it fails.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	346 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	347 450 JCUR = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	348 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	349 IF (M .GT. 0) DSM = SVNORM (N, ACOR, EWT)/TESCO(2,NQ)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	350 IF (DSM .GT. 1.0E0) GO TO 500
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	351 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	352 C After a successful step, update the YH array.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	353 C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	354 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	355 C use in a possible order increase on the next step.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	356 C If a change in H is considered, an increase or decrease in order
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	357 C by one is considered also. A change in H is made only if it is by a
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	358 C factor of at least 1.1. If not, IALTH is set to 3 to prevent
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	359 C testing for that many steps.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	360 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	361 KFLAG = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	362 IREDO = 0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	363 NST = NST + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	364 HU = H
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	365 NQU = NQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	366 DO 470 J = 1,L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	367 DO 470 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	368 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	369 IALTH = IALTH - 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	370 IF (IALTH .EQ. 0) GO TO 520
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	371 IF (IALTH .GT. 1) GO TO 700
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	372 IF (L .EQ. LMAX) GO TO 700
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	373 DO 490 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	374 490 YH(I,LMAX) = ACOR(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	375 GO TO 700
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	376 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	377 C The error test failed. KFLAG keeps track of multiple failures.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	378 C Restore TN and the YH array to their previous values, and prepare
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	379 C to try the step again. Compute the optimum step size for this or
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	380 C one lower order. After 2 or more failures, H is forced to decrease
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	381 C by a factor of 0.2 or less.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	382 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	383 500 KFLAG = KFLAG - 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	384 TN = TOLD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	385 I1 = NQNYH + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	386 DO 515 JB = 1,NQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	387 I1 = I1 - NYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	388 Cdir$ ivdep
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	389 DO 510 I = I1,NQNYH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	390 510 YH1(I) = YH1(I) - YH1(I+NYH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	391 515 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	392 RMAX = 2.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	393 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	394 IF (KFLAG .LE. -3) GO TO 640
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	395 IREDO = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	396 RHUP = 0.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	397 GO TO 540
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	398 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	399 C Regardless of the success or failure of the step, factors
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	400 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	401 C at order NQ - 1, order NQ, or order NQ + 1, respectively.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	402 C In the case of failure, RHUP = 0.0 to avoid an order increase.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	403 C The largest of these is determined and the new order chosen
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	404 C accordingly. If the order is to be increased, we compute one
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	405 C additional scaled derivative.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	406 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	407 520 RHUP = 0.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	408 IF (L .EQ. LMAX) GO TO 540
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	409 DO 530 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	410 530 SAVF(I) = ACOR(I) - YH(I,LMAX)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	411 DUP = SVNORM (N, SAVF, EWT)/TESCO(3,NQ)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	412 EXUP = 1.0E0/(L+1)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	413 RHUP = 1.0E0/(1.4E0DUP*EXUP + 0.0000014E0)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	414 540 EXSM = 1.0E0/L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	415 RHSM = 1.0E0/(1.2E0DSM*EXSM + 0.0000012E0)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	416 RHDN = 0.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	417 IF (NQ .EQ. 1) GO TO 560
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	418 DDN = SVNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	419 EXDN = 1.0E0/NQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	420 RHDN = 1.0E0/(1.3E0DDN*EXDN + 0.0000013E0)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	421 560 IF (RHSM .GE. RHUP) GO TO 570
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	422 IF (RHUP .GT. RHDN) GO TO 590
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	423 GO TO 580
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	424 570 IF (RHSM .LT. RHDN) GO TO 580
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	425 NEWQ = NQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	426 RH = RHSM
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	427 GO TO 620
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	428 580 NEWQ = NQ - 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	429 RH = RHDN
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	430 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	431 GO TO 620
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	432 590 NEWQ = L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	433 RH = RHUP
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	434 IF (RH .LT. 1.1E0) GO TO 610
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	435 R = EL(L)/L
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	436 DO 600 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	437 600 YH(I,NEWQ+1) = ACOR(I)*R
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	438 GO TO 630
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	439 610 IALTH = 3
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	440 GO TO 700
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	441 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	442 IF (KFLAG .LE. -2) RH = MIN(RH,0.2E0)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	443 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	444 C If there is a change of order, reset NQ, l, and the coefficients.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	445 C In any case H is reset according to RH and the YH array is rescaled.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	446 C Then exit from 690 if the step was OK, or redo the step otherwise.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	447 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	448 IF (NEWQ .EQ. NQ) GO TO 170
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	449 630 NQ = NEWQ
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	450 L = NQ + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	451 IRET = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	452 GO TO 150
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	453 C-----------------------------------------------------------------------
21003 f7e416862e90 doc: fix spelling of "occurred". Rafael Laboissiere <rafael@laboissiere.net> parents: 15271 diff changeset	454 C Control reaches this section if 3 or more failures have occurred.
7793 96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	455 C If 10 failures have occurred, exit with KFLAG = -1.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	456 C It is assumed that the derivatives that have accumulated in the
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	457 C YH array have errors of the wrong order. Hence the first
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	458 C derivative is recomputed, and the order is set to 1. Then
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	459 C H is reduced by a factor of 10, and the step is retried,
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	460 C until it succeeds or H reaches HMIN.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	461 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	462 640 IF (KFLAG .EQ. -10) GO TO 660
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	463 RH = 0.1E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	464 RH = MAX(HMIN/ABS(H),RH)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	465 H = H*RH
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	466 DO 645 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	467 645 Y(I) = YH(I,1)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	468 CALL F (NEQ, TN, Y, SAVF)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	469 NFE = NFE + 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	470 DO 650 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	471 650 YH(I,2) = H*SAVF(I)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	472 IPUP = MITER
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	473 IALTH = 5
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	474 IF (NQ .EQ. 1) GO TO 200
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	475 NQ = 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	476 L = 2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	477 IRET = 3
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	478 GO TO 150
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	479 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	480 C All returns are made through this section. H is saved in HOLD
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	481 C to allow the caller to change H on the next step.
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	482 C-----------------------------------------------------------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	483 660 KFLAG = -1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	484 GO TO 720
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	485 670 KFLAG = -2
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	486 GO TO 720
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	487 680 KFLAG = -3
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	488 GO TO 720
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	489 690 RMAX = 10.0E0
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	490 700 R = 1.0E0/TESCO(2,NQU)
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	491 DO 710 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	492 710 ACOR(I) = ACOR(I)*R
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	493 720 HOLD = H
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	494 JSTART = 1
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	495 RETURN
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	496 C----------------------- END OF SUBROUTINE SSTODE ----------------------
96ba591be50f Add some more support for single precision to libcruft functions David Bateman <dbateman@free.fr> parents: diff changeset	497 END

Mercurial > hg > octave-jordi

annotate liboctave/cruft/odepack/sstode.f @ 21003:f7e416862e90 stable