Mercurial > octave-nkf
changeset 7793:96ba591be50f
Add some more support for single precision to libcruft functions
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sun, 11 May 2008 22:51:50 +0200 |
| parents | 39c1026191e9 |
| children | 2b458dfe31ae |
| files | libcruft/ChangeLog libcruft/odepack/Makefile.in libcruft/odepack/scfode.f libcruft/odepack/sewset.f libcruft/odepack/sintdy.f libcruft/odepack/slsode.f libcruft/odepack/sprepj.f libcruft/odepack/ssolsy.f libcruft/odepack/sstode.f libcruft/odepack/svnorm.f libcruft/ordered-qz/Makefile.in libcruft/ordered-qz/sexchqz.f libcruft/ordered-qz/ssubsp.f libcruft/quadpack/Makefile.in libcruft/quadpack/qagi.f libcruft/quadpack/qagie.f libcruft/quadpack/qagp.f libcruft/quadpack/qagpe.f libcruft/quadpack/qelg.f libcruft/quadpack/qk15i.f libcruft/quadpack/qk21.f libcruft/quadpack/qpsrt.f liboctave/ChangeLog |
| diffstat | 23 files changed, 5410 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Wed May 21 09:36:46 2008 -0400 +++ b/libcruft/ChangeLog Sun May 11 22:51:50 2008 +0200 @@ -1,3 +1,17 @@ +2008-05-21 David Bateman <dbateman@free.fr> + + * odepack/scfode.f, odepack/sewset.f, odepack/sintdy.f, + odepack/slsode.f, odepack/sprepj.f, odepack/ssolsy.f, + odepack/sstode.f, odepack/svnorm.f: New files. + * odepack/Makefile.in (FSRC): Add them. + + * ordered-qz/sexchqz.f, ordered-qz/ssubsp.f: New files. + * ordered-qz/Makefile.in (FSRC): Add them. + * quadpack/qagi.f, quadpack/qagie.f, quadpack/qagp.f, + quadpack/qagpe.f, quadpack/qelg.f, quadpack/qk15i.f, + quadpack/qk21.f, quadpack/qpsrt.f: New files. + * quadpack/Makefile.in (FSRC): Add them. + 2008-05-20 Jaroslav Hajek <highegg@gmail.com> * qrupdate/cch1dn.f, qrupdate/cchinx.f, qrupdate/cqhqr.f,
--- a/libcruft/odepack/Makefile.in Wed May 21 09:36:46 2008 -0400 +++ b/libcruft/odepack/Makefile.in Sun May 11 22:51:50 2008 +0200 @@ -26,7 +26,8 @@ EXTERNAL_DISTFILES = $(DISTFILES) -FSRC = cfode.f dlsode.f ewset.f intdy.f prepj.f solsy.f stode.f vnorm.f +FSRC = cfode.f dlsode.f ewset.f intdy.f prepj.f solsy.f stode.f vnorm.f \ + scfode.f sewset.f sintdy.f slsode.f sprepj.f ssolsy.f sstode.f svnorm.f include $(TOPDIR)/Makeconf
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/scfode.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,127 @@ + SUBROUTINE SCFODE (METH, ELCO, TESCO) +C***BEGIN PROLOGUE SCFODE +C***SUBSIDIARY +C***PURPOSE Set ODE integrator coefficients. +C***TYPE SINGLE PRECISION (SCFODE-S, DCFODE-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C SCFODE is called by the integrator routine to set coefficients +C needed there. The coefficients for the current method, as +C given by the value of METH, are set for all orders and saved. +C The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2. +C (A smaller value of the maximum order is also allowed.) +C SCFODE is called once at the beginning of the problem, +C and is not called again unless and until METH is changed. +C +C The ELCO array contains the basic method coefficients. +C The coefficients el(i), 1 .le. i .le. nq+1, for the method of +C order nq are stored in ELCO(i,nq). They are given by a genetrating +C polynomial, i.e., +C l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq. +C For the implicit Adams methods, l(x) is given by +C dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0. +C For the BDF methods, l(x) is given by +C l(x) = (x+1)*(x+2)* ... *(x+nq)/K, +C where K = factorial(nq)*(1 + 1/2 + ... + 1/nq). +C +C The TESCO array contains test constants used for the +C local error test and the selection of step size and/or order. +C At order nq, TESCO(k,nq) is used for the selection of step +C size at order nq - 1 if k = 1, at order nq if k = 2, and at order +C nq + 1 if k = 3. +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED (NONE) +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890503 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C***END PROLOGUE SCFODE +C**End + INTEGER METH + INTEGER I, IB, NQ, NQM1, NQP1 + REAL ELCO, TESCO + REAL AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ, + 1 RQFAC, RQ1FAC, TSIGN, XPIN + DIMENSION ELCO(13,12), TESCO(3,12) + DIMENSION PC(12) +C +C***FIRST EXECUTABLE STATEMENT SCFODE + GO TO (100, 200), METH +C + 100 ELCO(1,1) = 1.0E0 + ELCO(2,1) = 1.0E0 + TESCO(1,1) = 0.0E0 + TESCO(2,1) = 2.0E0 + TESCO(1,2) = 1.0E0 + TESCO(3,12) = 0.0E0 + PC(1) = 1.0E0 + RQFAC = 1.0E0 + DO 140 NQ = 2,12 +C----------------------------------------------------------------------- +C The PC array will contain the coefficients of the polynomial +C p(x) = (x+1)*(x+2)*...*(x+nq-1). +C Initially, p(x) = 1. +C----------------------------------------------------------------------- + RQ1FAC = RQFAC + RQFAC = RQFAC/NQ + NQM1 = NQ - 1 + FNQM1 = NQM1 + NQP1 = NQ + 1 +C Form coefficients of p(x)*(x+nq-1). ---------------------------------- + PC(NQ) = 0.0E0 + DO 110 IB = 1,NQM1 + I = NQP1 - IB + 110 PC(I) = PC(I-1) + FNQM1*PC(I) + PC(1) = FNQM1*PC(1) +C Compute integral, -1 to 0, of p(x) and x*p(x). ----------------------- + PINT = PC(1) + XPIN = PC(1)/2.0E0 + TSIGN = 1.0E0 + DO 120 I = 2,NQ + TSIGN = -TSIGN + PINT = PINT + TSIGN*PC(I)/I + 120 XPIN = XPIN + TSIGN*PC(I)/(I+1) +C Store coefficients in ELCO and TESCO. -------------------------------- + ELCO(1,NQ) = PINT*RQ1FAC + ELCO(2,NQ) = 1.0E0 + DO 130 I = 2,NQ + 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I + AGAMQ = RQFAC*XPIN + RAGQ = 1.0E0/AGAMQ + TESCO(2,NQ) = RAGQ + IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1 + TESCO(3,NQM1) = RAGQ + 140 CONTINUE + RETURN +C + 200 PC(1) = 1.0E0 + RQ1FAC = 1.0E0 + DO 230 NQ = 1,5 +C----------------------------------------------------------------------- +C The PC array will contain the coefficients of the polynomial +C p(x) = (x+1)*(x+2)*...*(x+nq). +C Initially, p(x) = 1. +C----------------------------------------------------------------------- + FNQ = NQ + NQP1 = NQ + 1 +C Form coefficients of p(x)*(x+nq). ------------------------------------ + PC(NQP1) = 0.0E0 + DO 210 IB = 1,NQ + I = NQ + 2 - IB + 210 PC(I) = PC(I-1) + FNQ*PC(I) + PC(1) = FNQ*PC(1) +C Store coefficients in ELCO and TESCO. -------------------------------- + DO 220 I = 1,NQP1 + 220 ELCO(I,NQ) = PC(I)/PC(2) + ELCO(2,NQ) = 1.0E0 + TESCO(1,NQ) = RQ1FAC + TESCO(2,NQ) = NQP1/ELCO(1,NQ) + TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ) + RQ1FAC = RQ1FAC/FNQ + 230 CONTINUE + RETURN +C----------------------- END OF SUBROUTINE SCFODE ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/sewset.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,47 @@ + SUBROUTINE SEWSET (N, ITOL, RTOL, ATOL, YCUR, EWT) +C***BEGIN PROLOGUE SEWSET +C***SUBSIDIARY +C***PURPOSE Set error weight vector. +C***TYPE SINGLE PRECISION (SEWSET-S, DEWSET-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C This subroutine sets the error weight vector EWT according to +C EWT(i) = RTOL(i)*ABS(YCUR(i)) + ATOL(i), i = 1,...,N, +C with the subscript on RTOL and/or ATOL possibly replaced by 1 above, +C depending on the value of ITOL. +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED (NONE) +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890503 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C***END PROLOGUE SEWSET +C**End + INTEGER N, ITOL + INTEGER I + REAL RTOL, ATOL, YCUR, EWT + DIMENSION RTOL(*), ATOL(*), YCUR(N), EWT(N) +C +C***FIRST EXECUTABLE STATEMENT SEWSET + GO TO (10, 20, 30, 40), ITOL + 10 CONTINUE + DO 15 I = 1,N + 15 EWT(I) = RTOL(1)*ABS(YCUR(I)) + ATOL(1) + RETURN + 20 CONTINUE + DO 25 I = 1,N + 25 EWT(I) = RTOL(1)*ABS(YCUR(I)) + ATOL(I) + RETURN + 30 CONTINUE + DO 35 I = 1,N + 35 EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(1) + RETURN + 40 CONTINUE + DO 45 I = 1,N + 45 EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(I) + RETURN +C----------------------- END OF SUBROUTINE SEWSET ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/sintdy.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,107 @@ + SUBROUTINE SINTDY (T, K, YH, NYH, DKY, IFLAG) +C***BEGIN PROLOGUE SINTDY +C***SUBSIDIARY +C***PURPOSE Interpolate solution derivatives. +C***TYPE SINGLE PRECISION (SINTDY-S, DINTDY-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C SINTDY computes interpolated values of the K-th derivative of the +C dependent variable vector y, and stores it in DKY. This routine +C is called within the package with K = 0 and T = TOUT, but may +C also be called by the user for any K up to the current order. +C (See detailed instructions in the usage documentation.) +C +C The computed values in DKY are gotten by interpolation using the +C Nordsieck history array YH. This array corresponds uniquely to a +C vector-valued polynomial of degree NQCUR or less, and DKY is set +C to the K-th derivative of this polynomial at T. +C The formula for DKY is: +C q +C DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1) +C j=K +C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR. +C The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are +C communicated by COMMON. The above sum is done in reverse order. +C IFLAG is returned negative if either K or T is out of bounds. +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED XERRWV +C***COMMON BLOCKS SLS001 +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890503 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C 010412 Reduced size of Common block /SLS001/. (ACH) +C 031105 Restored 'own' variables to Common block /SLS001/, to +C enable interrupt/restart feature. (ACH) +C 050427 Corrected roundoff decrement in TP. (ACH) +C***END PROLOGUE SINTDY +C**End + INTEGER K, NYH, IFLAG + REAL T, YH, DKY + DIMENSION YH(NYH,*), DKY(*) + INTEGER IOWND, IOWNS, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + REAL ROWNS, + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + COMMON /SLS001/ ROWNS(209), + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 2 IOWND(6), IOWNS(6), + 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1 + REAL C, R, S, TP + CHARACTER*80 MSG +C +C***FIRST EXECUTABLE STATEMENT SINTDY + IFLAG = 0 + IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80 + TP = TN - HU - 100.0E0*UROUND*SIGN(ABS(TN) + ABS(HU), HU) + IF ((T-TP)*(T-TN) .GT. 0.0E0) GO TO 90 +C + S = (T - TN)/H + IC = 1 + IF (K .EQ. 0) GO TO 15 + JJ1 = L - K + DO 10 JJ = JJ1,NQ + 10 IC = IC*JJ + 15 C = IC + DO 20 I = 1,N + 20 DKY(I) = C*YH(I,L) + IF (K .EQ. NQ) GO TO 55 + JB2 = NQ - K + DO 50 JB = 1,JB2 + J = NQ - JB + JP1 = J + 1 + IC = 1 + IF (K .EQ. 0) GO TO 35 + JJ1 = JP1 - K + DO 30 JJ = JJ1,J + 30 IC = IC*JJ + 35 C = IC + DO 40 I = 1,N + 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) + 50 CONTINUE + IF (K .EQ. 0) RETURN + 55 R = H**(-K) + DO 60 I = 1,N + 60 DKY(I) = R*DKY(I) + RETURN +C + 80 MSG = 'SINTDY- K (=I1) illegal ' + CALL XERRWV (MSG, 30, 51, 0, 1, K, 0, 0, 0.0E0, 0.0E0) + IFLAG = -1 + RETURN + 90 MSG = 'SINTDY- T (=R1) illegal ' + CALL XERRWV (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0E0) + MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) ' + CALL XERRWV (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN) + IFLAG = -2 + RETURN +C----------------------- END OF SUBROUTINE SINTDY ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/slsode.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,1755 @@ +*DECK SLSODE + SUBROUTINE SLSODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, + 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) + EXTERNAL F, JAC + INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, MF + REAL Y, T, TOUT, RTOL, ATOL, RWORK + DIMENSION NEQ(*), Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW) +C***BEGIN PROLOGUE SLSODE +C***PURPOSE Livermore Solver for Ordinary Differential Equations. +C SLSODE solves the initial-value problem for stiff or +C nonstiff systems of first-order ODE's, +C dy/dt = f(t,y), or, in component form, +C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(N)), i=1,...,N. +C***CATEGORY I1A +C***TYPE SINGLE PRECISION (SLSODE-S, DLSODE-D) +C***KEYWORDS ORDINARY DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM, +C STIFF, NONSTIFF +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C Center for Applied Scientific Computing, L-561 +C Lawrence Livermore National Laboratory +C Livermore, CA 94551. +C***DESCRIPTION +C +C NOTE: The "Usage" and "Arguments" sections treat only a subset of +C available options, in condensed fashion. The options +C covered and the information supplied will support most +C standard uses of SLSODE. +C +C For more sophisticated uses, full details on all options are +C given in the concluding section, headed "Long Description." +C A synopsis of the SLSODE Long Description is provided at the +C beginning of that section; general topics covered are: +C - Elements of the call sequence; optional input and output +C - Optional supplemental routines in the SLSODE package +C - internal COMMON block +C +C *Usage: +C Communication between the user and the SLSODE package, for normal +C situations, is summarized here. This summary describes a subset +C of the available options. See "Long Description" for complete +C details, including optional communication, nonstandard options, +C and instructions for special situations. +C +C A sample program is given in the "Examples" section. +C +C Refer to the argument descriptions for the definitions of the +C quantities that appear in the following sample declarations. +C +C For MF = 10, +C PARAMETER (LRW = 20 + 16*NEQ, LIW = 20) +C For MF = 21 or 22, +C PARAMETER (LRW = 22 + 9*NEQ + NEQ**2, LIW = 20 + NEQ) +C For MF = 24 or 25, +C PARAMETER (LRW = 22 + 10*NEQ + (2*ML+MU)*NEQ, +C * LIW = 20 + NEQ) +C +C EXTERNAL F, JAC +C INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK(LIW), +C * LIW, MF +C REAL Y(NEQ), T, TOUT, RTOL, ATOL(ntol), RWORK(LRW) +C +C CALL SLSODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, +C * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) +C +C *Arguments: +C F :EXT Name of subroutine for right-hand-side vector f. +C This name must be declared EXTERNAL in calling +C program. The form of F must be: +C +C SUBROUTINE F (NEQ, T, Y, YDOT) +C INTEGER NEQ +C REAL T, Y(*), YDOT(*) +C +C The inputs are NEQ, T, Y. F is to set +C +C YDOT(i) = f(i,T,Y(1),Y(2),...,Y(NEQ)), +C i = 1, ..., NEQ . +C +C NEQ :IN Number of first-order ODE's. +C +C Y :INOUT Array of values of the y(t) vector, of length NEQ. +C Input: For the first call, Y should contain the +C values of y(t) at t = T. (Y is an input +C variable only if ISTATE = 1.) +C Output: On return, Y will contain the values at the +C new t-value. +C +C T :INOUT Value of the independent variable. On return it +C will be the current value of t (normally TOUT). +C +C TOUT :IN Next point where output is desired (.NE. T). +C +C ITOL :IN 1 or 2 according as ATOL (below) is a scalar or +C an array. +C +C RTOL :IN Relative tolerance parameter (scalar). +C +C ATOL :IN Absolute tolerance parameter (scalar or array). +C If ITOL = 1, ATOL need not be dimensioned. +C If ITOL = 2, ATOL must be dimensioned at least NEQ. +C +C The estimated local error in Y(i) will be controlled +C so as to be roughly less (in magnitude) than +C +C EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or +C EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2. +C +C Thus the local error test passes if, in each +C component, either the absolute error is less than +C ATOL (or ATOL(i)), or the relative error is less +C than RTOL. +C +C Use RTOL = 0.0 for pure absolute error control, and +C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative +C error control. Caution: Actual (global) errors may +C exceed these local tolerances, so choose them +C conservatively. +C +C ITASK :IN Flag indicating the task SLSODE is to perform. +C Use ITASK = 1 for normal computation of output +C values of y at t = TOUT. +C +C ISTATE:INOUT Index used for input and output to specify the state +C of the calculation. +C Input: +C 1 This is the first call for a problem. +C 2 This is a subsequent call. +C Output: +C 1 Nothing was done, as TOUT was equal to T. +C 2 SLSODE was successful (otherwise, negative). +C Note that ISTATE need not be modified after a +C successful return. +C -1 Excess work done on this call (perhaps wrong +C MF). +C -2 Excess accuracy requested (tolerances too +C small). +C -3 Illegal input detected (see printed message). +C -4 Repeated error test failures (check all +C inputs). +C -5 Repeated convergence failures (perhaps bad +C Jacobian supplied or wrong choice of MF or +C tolerances). +C -6 Error weight became zero during problem +C (solution component i vanished, and ATOL or +C ATOL(i) = 0.). +C +C IOPT :IN Flag indicating whether optional inputs are used: +C 0 No. +C 1 Yes. (See "Optional inputs" under "Long +C Description," Part 1.) +C +C RWORK :WORK Real work array of length at least: +C 20 + 16*NEQ for MF = 10, +C 22 + 9*NEQ + NEQ**2 for MF = 21 or 22, +C 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25. +C +C LRW :IN Declared length of RWORK (in user's DIMENSION +C statement). +C +C IWORK :WORK Integer work array of length at least: +C 20 for MF = 10, +C 20 + NEQ for MF = 21, 22, 24, or 25. +C +C If MF = 24 or 25, input in IWORK(1),IWORK(2) the +C lower and upper Jacobian half-bandwidths ML,MU. +C +C On return, IWORK contains information that may be +C of interest to the user: +C +C Name Location Meaning +C ----- --------- ----------------------------------------- +C NST IWORK(11) Number of steps taken for the problem so +C far. +C NFE IWORK(12) Number of f evaluations for the problem +C so far. +C NJE IWORK(13) Number of Jacobian evaluations (and of +C matrix LU decompositions) for the problem +C so far. +C NQU IWORK(14) Method order last used (successfully). +C LENRW IWORK(17) Length of RWORK actually required. This +C is defined on normal returns and on an +C illegal input return for insufficient +C storage. +C LENIW IWORK(18) Length of IWORK actually required. This +C is defined on normal returns and on an +C illegal input return for insufficient +C storage. +C +C LIW :IN Declared length of IWORK (in user's DIMENSION +C statement). +C +C JAC :EXT Name of subroutine for Jacobian matrix (MF = +C 21 or 24). If used, this name must be declared +C EXTERNAL in calling program. If not used, pass a +C dummy name. The form of JAC must be: +C +C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) +C INTEGER NEQ, ML, MU, NROWPD +C REAL T, Y(*), PD(NROWPD,*) +C +C See item c, under "Description" below for more +C information about JAC. +C +C MF :IN Method flag. Standard values are: +C 10 Nonstiff (Adams) method, no Jacobian used. +C 21 Stiff (BDF) method, user-supplied full Jacobian. +C 22 Stiff method, internally generated full +C Jacobian. +C 24 Stiff method, user-supplied banded Jacobian. +C 25 Stiff method, internally generated banded +C Jacobian. +C +C *Description: +C SLSODE solves the initial value problem for stiff or nonstiff +C systems of first-order ODE's, +C +C dy/dt = f(t,y) , +C +C or, in component form, +C +C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) +C (i = 1, ..., NEQ) . +C +C SLSODE is a package based on the GEAR and GEARB packages, and on +C the October 23, 1978, version of the tentative ODEPACK user +C interface standard, with minor modifications. +C +C The steps in solving such a problem are as follows. +C +C a. First write a subroutine of the form +C +C SUBROUTINE F (NEQ, T, Y, YDOT) +C INTEGER NEQ +C REAL T, Y(*), YDOT(*) +C +C which supplies the vector function f by loading YDOT(i) with +C f(i). +C +C b. Next determine (or guess) whether or not the problem is stiff. +C Stiffness occurs when the Jacobian matrix df/dy has an +C eigenvalue whose real part is negative and large in magnitude +C compared to the reciprocal of the t span of interest. If the +C problem is nonstiff, use method flag MF = 10. If it is stiff, +C there are four standard choices for MF, and SLSODE requires the +C Jacobian matrix in some form. This matrix is regarded either +C as full (MF = 21 or 22), or banded (MF = 24 or 25). In the +C banded case, SLSODE requires two half-bandwidth parameters ML +C and MU. These are, respectively, the widths of the lower and +C upper parts of the band, excluding the main diagonal. Thus the +C band consists of the locations (i,j) with +C +C i - ML <= j <= i + MU , +C +C and the full bandwidth is ML + MU + 1 . +C +C c. If the problem is stiff, you are encouraged to supply the +C Jacobian directly (MF = 21 or 24), but if this is not feasible, +C SLSODE will compute it internally by difference quotients (MF = +C 22 or 25). If you are supplying the Jacobian, write a +C subroutine of the form +C +C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) +C INTEGER NEQ, ML, MU, NRWOPD +C REAL T, Y(*), PD(NROWPD,*) +C +C which provides df/dy by loading PD as follows: +C - For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j), +C the partial derivative of f(i) with respect to y(j). (Ignore +C the ML and MU arguments in this case.) +C - For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with +C df(i)/dy(j); i.e., load the diagonal lines of df/dy into the +C rows of PD from the top down. +C - In either case, only nonzero elements need be loaded. +C +C d. Write a main program that calls subroutine SLSODE once for each +C point at which answers are desired. This should also provide +C for possible use of logical unit 6 for output of error messages +C by SLSODE. +C +C Before the first call to SLSODE, set ISTATE = 1, set Y and T to +C the initial values, and set TOUT to the first output point. To +C continue the integration after a successful return, simply +C reset TOUT and call SLSODE again. No other parameters need be +C reset. +C +C *Examples: +C The following is a simple example problem, with the coding needed +C for its solution by SLSODE. The problem is from chemical kinetics, +C and consists of the following three rate equations: +C +C dy1/dt = -.04*y1 + 1.E4*y2*y3 +C dy2/dt = .04*y1 - 1.E4*y2*y3 - 3.E7*y2**2 +C dy3/dt = 3.E7*y2**2 +C +C on the interval from t = 0.0 to t = 4.E10, with initial conditions +C y1 = 1.0, y2 = y3 = 0. The problem is stiff. +C +C The following coding solves this problem with SLSODE, using +C MF = 21 and printing results at t = .4, 4., ..., 4.E10. It uses +C ITOL = 2 and ATOL much smaller for y2 than for y1 or y3 because y2 +C has much smaller values. At the end of the run, statistical +C quantities of interest are printed. +C +C EXTERNAL FEX, JEX +C INTEGER IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(23), LIW, LRW, +C * MF, NEQ +C REAL ATOL(3), RTOL, RWORK(58), T, TOUT, Y(3) +C NEQ = 3 +C Y(1) = 1. +C Y(2) = 0. +C Y(3) = 0. +C T = 0. +C TOUT = .4 +C ITOL = 2 +C RTOL = 1.E-4 +C ATOL(1) = 1.E-6 +C ATOL(2) = 1.E-10 +C ATOL(3) = 1.E-6 +C ITASK = 1 +C ISTATE = 1 +C IOPT = 0 +C LRW = 58 +C LIW = 23 +C MF = 21 +C DO 40 IOUT = 1,12 +C CALL SLSODE (FEX, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, +C * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF) +C WRITE(6,20) T, Y(1), Y(2), Y(3) +C 20 FORMAT(' At t =',E12.4,' y =',3E14.6) +C IF (ISTATE .LT. 0) GO TO 80 +C 40 TOUT = TOUT*10. +C WRITE(6,60) IWORK(11), IWORK(12), IWORK(13) +C 60 FORMAT(/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4) +C STOP +C 80 WRITE(6,90) ISTATE +C 90 FORMAT(///' Error halt.. ISTATE =',I3) +C STOP +C END +C +C SUBROUTINE FEX (NEQ, T, Y, YDOT) +C INTEGER NEQ +C REAL T, Y(3), YDOT(3) +C YDOT(1) = -.04*Y(1) + 1.E4*Y(2)*Y(3) +C YDOT(3) = 3.E7*Y(2)*Y(2) +C YDOT(2) = -YDOT(1) - YDOT(3) +C RETURN +C END +C +C SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD) +C INTEGER NEQ, ML, MU, NRPD +C REAL T, Y(3), PD(NRPD,3) +C PD(1,1) = -.04 +C PD(1,2) = 1.E4*Y(3) +C PD(1,3) = 1.E4*Y(2) +C PD(2,1) = .04 +C PD(2,3) = -PD(1,3) +C PD(3,2) = 6.E7*Y(2) +C PD(2,2) = -PD(1,2) - PD(3,2) +C RETURN +C END +C +C The output from this program (on a Cray-1 in single precision) +C is as follows. +C +C At t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02 +C At t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02 +C At t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01 +C At t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01 +C At t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01 +C At t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01 +C At t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01 +C At t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01 +C At t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01 +C At t = 4.0000e+08 y = 5.494530e-06 2.197825e-11 9.999945e-01 +C At t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01 +C At t = 4.0000e+10 y = -7.170603e-08 -2.868241e-13 1.000000e+00 +C +C No. steps = 330, No. f-s = 405, No. J-s = 69 +C +C *Accuracy: +C The accuracy of the solution depends on the choice of tolerances +C RTOL and ATOL. Actual (global) errors may exceed these local +C tolerances, so choose them conservatively. +C +C *Cautions: +C The work arrays should not be altered between calls to SLSODE for +C the same problem, except possibly for the conditional and optional +C inputs. +C +C *Portability: +C Since NEQ is dimensioned inside SLSODE, some compilers may object +C to a call to SLSODE with NEQ a scalar variable. In this event, +C use DIMENSION NEQ(1). Similar remarks apply to RTOL and ATOL. +C +C Note to Cray users: +C For maximum efficiency, use the CFT77 compiler. Appropriate +C compiler optimization directives have been inserted for CFT77. +C +C *Reference: +C Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE +C Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds. +C (North-Holland, Amsterdam, 1983), pp. 55-64. +C +C *Long Description: +C The following complete description of the user interface to +C SLSODE consists of four parts: +C +C 1. The call sequence to subroutine SLSODE, which is a driver +C routine for the solver. This includes descriptions of both +C the call sequence arguments and user-supplied routines. +C Following these descriptions is a description of optional +C inputs available through the call sequence, and then a +C description of optional outputs in the work arrays. +C +C 2. Descriptions of other routines in the SLSODE package that may +C be (optionally) called by the user. These provide the ability +C to alter error message handling, save and restore the internal +C COMMON, and obtain specified derivatives of the solution y(t). +C +C 3. Descriptions of COMMON block to be declared in overlay or +C similar environments, or to be saved when doing an interrupt +C of the problem and continued solution later. +C +C 4. Description of two routines in the SLSODE package, either of +C which the user may replace with his own version, if desired. +C These relate to the measurement of errors. +C +C +C Part 1. Call Sequence +C ---------------------- +C +C Arguments +C --------- +C The call sequence parameters used for input only are +C +C F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF, +C +C and those used for both input and output are +C +C Y, T, ISTATE. +C +C The work arrays RWORK and IWORK are also used for conditional and +C optional inputs and optional outputs. (The term output here +C refers to the return from subroutine SLSODE to the user's calling +C program.) +C +C The legality of input parameters will be thoroughly checked on the +C initial call for the problem, but not checked thereafter unless a +C change in input parameters is flagged by ISTATE = 3 on input. +C +C The descriptions of the call arguments are as follows. +C +C F The name of the user-supplied subroutine defining the ODE +C system. The system must be put in the first-order form +C dy/dt = f(t,y), where f is a vector-valued function of +C the scalar t and the vector y. Subroutine F is to compute +C the function f. It is to have the form +C +C SUBROUTINE F (NEQ, T, Y, YDOT) +C REAL T, Y(*), YDOT(*) +C +C where NEQ, T, and Y are input, and the array YDOT = +C f(T,Y) is output. Y and YDOT are arrays of length NEQ. +C Subroutine F should not alter Y(1),...,Y(NEQ). F must be +C declared EXTERNAL in the calling program. +C +C Subroutine F may access user-defined quantities in +C NEQ(2),... and/or in Y(NEQ(1)+1),..., if NEQ is an array +C (dimensioned in F) and/or Y has length exceeding NEQ(1). +C See the descriptions of NEQ and Y below. +C +C If quantities computed in the F routine are needed +C externally to SLSODE, an extra call to F should be made +C for this purpose, for consistent and accurate results. +C If only the derivative dy/dt is needed, use SINTDY +C instead. +C +C NEQ The size of the ODE system (number of first-order +C ordinary differential equations). Used only for input. +C NEQ may be decreased, but not increased, during the +C problem. If NEQ is decreased (with ISTATE = 3 on input), +C the remaining components of Y should be left undisturbed, +C if these are to be accessed in F and/or JAC. +C +C Normally, NEQ is a scalar, and it is generally referred +C to as a scalar in this user interface description. +C However, NEQ may be an array, with NEQ(1) set to the +C system size. (The SLSODE package accesses only NEQ(1).) +C In either case, this parameter is passed as the NEQ +C argument in all calls to F and JAC. Hence, if it is an +C array, locations NEQ(2),... may be used to store other +C integer data and pass it to F and/or JAC. Subroutines +C F and/or JAC must include NEQ in a DIMENSION statement +C in that case. +C +C Y A real array for the vector of dependent variables, of +C length NEQ or more. Used for both input and output on +C the first call (ISTATE = 1), and only for output on +C other calls. On the first call, Y must contain the +C vector of initial values. On output, Y contains the +C computed solution vector, evaluated at T. If desired, +C the Y array may be used for other purposes between +C calls to the solver. +C +C This array is passed as the Y argument in all calls to F +C and JAC. Hence its length may exceed NEQ, and locations +C Y(NEQ+1),... may be used to store other real data and +C pass it to F and/or JAC. (The SLSODE package accesses +C only Y(1),...,Y(NEQ).) +C +C T The independent variable. On input, T is used only on +C the first call, as the initial point of the integration. +C On output, after each call, T is the value at which a +C computed solution Y is evaluated (usually the same as +C TOUT). On an error return, T is the farthest point +C reached. +C +C TOUT The next value of T at which a computed solution is +C desired. Used only for input. +C +C When starting the problem (ISTATE = 1), TOUT may be equal +C to T for one call, then should not equal T for the next +C call. For the initial T, an input value of TOUT .NE. T +C is used in order to determine the direction of the +C integration (i.e., the algebraic sign of the step sizes) +C and the rough scale of the problem. Integration in +C either direction (forward or backward in T) is permitted. +C +C If ITASK = 2 or 5 (one-step modes), TOUT is ignored +C after the first call (i.e., the first call with +C TOUT .NE. T). Otherwise, TOUT is required on every call. +C +C If ITASK = 1, 3, or 4, the values of TOUT need not be +C monotone, but a value of TOUT which backs up is limited +C to the current internal T interval, whose endpoints are +C TCUR - HU and TCUR. (See "Optional Outputs" below for +C TCUR and HU.) +C +C +C ITOL An indicator for the type of error control. See +C description below under ATOL. Used only for input. +C +C RTOL A relative error tolerance parameter, either a scalar or +C an array of length NEQ. See description below under +C ATOL. Input only. +C +C ATOL An absolute error tolerance parameter, either a scalar or +C an array of length NEQ. Input only. +C +C The input parameters ITOL, RTOL, and ATOL determine the +C error control performed by the solver. The solver will +C control the vector e = (e(i)) of estimated local errors +C in Y, according to an inequality of the form +C +C rms-norm of ( e(i)/EWT(i) ) <= 1, +C +C where +C +C EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i), +C +C and the rms-norm (root-mean-square norm) here is +C +C rms-norm(v) = SQRT(sum v(i)**2 / NEQ). +C +C Here EWT = (EWT(i)) is a vector of weights which must +C always be positive, and the values of RTOL and ATOL +C should all be nonnegative. The following table gives the +C types (scalar/array) of RTOL and ATOL, and the +C corresponding form of EWT(i). +C +C ITOL RTOL ATOL EWT(i) +C ---- ------ ------ ----------------------------- +C 1 scalar scalar RTOL*ABS(Y(i)) + ATOL +C 2 scalar array RTOL*ABS(Y(i)) + ATOL(i) +C 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL +C 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i) +C +C When either of these parameters is a scalar, it need not +C be dimensioned in the user's calling program. +C +C If none of the above choices (with ITOL, RTOL, and ATOL +C fixed throughout the problem) is suitable, more general +C error controls can be obtained by substituting +C user-supplied routines for the setting of EWT and/or for +C the norm calculation. See Part 4 below. +C +C If global errors are to be estimated by making a repeated +C run on the same problem with smaller tolerances, then all +C components of RTOL and ATOL (i.e., of EWT) should be +C scaled down uniformly. +C +C ITASK An index specifying the task to be performed. Input +C only. ITASK has the following values and meanings: +C 1 Normal computation of output values of y(t) at +C t = TOUT (by overshooting and interpolating). +C 2 Take one step only and return. +C 3 Stop at the first internal mesh point at or beyond +C t = TOUT and return. +C 4 Normal computation of output values of y(t) at +C t = TOUT but without overshooting t = TCRIT. TCRIT +C must be input as RWORK(1). TCRIT may be equal to or +C beyond TOUT, but not behind it in the direction of +C integration. This option is useful if the problem +C has a singularity at or beyond t = TCRIT. +C 5 Take one step, without passing TCRIT, and return. +C TCRIT must be input as RWORK(1). +C +C Note: If ITASK = 4 or 5 and the solver reaches TCRIT +C (within roundoff), it will return T = TCRIT (exactly) to +C indicate this (unless ITASK = 4 and TOUT comes before +C TCRIT, in which case answers at T = TOUT are returned +C first). +C +C ISTATE An index used for input and output to specify the state +C of the calculation. +C +C On input, the values of ISTATE are as follows: +C 1 This is the first call for the problem +C (initializations will be done). See "Note" below. +C 2 This is not the first call, and the calculation is to +C continue normally, with no change in any input +C parameters except possibly TOUT and ITASK. (If ITOL, +C RTOL, and/or ATOL are changed between calls with +C ISTATE = 2, the new values will be used but not +C tested for legality.) +C 3 This is not the first call, and the calculation is to +C continue normally, but with a change in input +C parameters other than TOUT and ITASK. Changes are +C allowed in NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, +C ML, MU, and any of the optional inputs except H0. +C (See IWORK description for ML and MU.) +C +C Note: A preliminary call with TOUT = T is not counted as +C a first call here, as no initialization or checking of +C input is done. (Such a call is sometimes useful for the +C purpose of outputting the initial conditions.) Thus the +C first call for which TOUT .NE. T requires ISTATE = 1 on +C input. +C +C On output, ISTATE has the following values and meanings: +C 1 Nothing was done, as TOUT was equal to T with +C ISTATE = 1 on input. +C 2 The integration was performed successfully. +C -1 An excessive amount of work (more than MXSTEP steps) +C was done on this call, before completing the +C requested task, but the integration was otherwise +C successful as far as T. (MXSTEP is an optional input +C and is normally 500.) To continue, the user may +C simply reset ISTATE to a value >1 and call again (the +C excess work step counter will be reset to 0). In +C addition, the user may increase MXSTEP to avoid this +C error return; see "Optional Inputs" below. +C -2 Too much accuracy was requested for the precision of +C the machine being used. This was detected before +C completing the requested task, but the integration +C was successful as far as T. To continue, the +C tolerance parameters must be reset, and ISTATE must +C be set to 3. The optional output TOLSF may be used +C for this purpose. (Note: If this condition is +C detected before taking any steps, then an illegal +C input return (ISTATE = -3) occurs instead.) +C -3 Illegal input was detected, before taking any +C integration steps. See written message for details. +C (Note: If the solver detects an infinite loop of +C calls to the solver with illegal input, it will cause +C the run to stop.) +C -4 There were repeated error-test failures on one +C attempted step, before completing the requested task, +C but the integration was successful as far as T. The +C problem may have a singularity, or the input may be +C inappropriate. +C -5 There were repeated convergence-test failures on one +C attempted step, before completing the requested task, +C but the integration was successful as far as T. This +C may be caused by an inaccurate Jacobian matrix, if +C one is being used. +C -6 EWT(i) became zero for some i during the integration. +C Pure relative error control (ATOL(i)=0.0) was +C requested on a variable which has now vanished. The +C integration was successful as far as T. +C +C Note: Since the normal output value of ISTATE is 2, it +C does not need to be reset for normal continuation. Also, +C since a negative input value of ISTATE will be regarded +C as illegal, a negative output value requires the user to +C change it, and possibly other inputs, before calling the +C solver again. +C +C IOPT An integer flag to specify whether any optional inputs +C are being used on this call. Input only. The optional +C inputs are listed under a separate heading below. +C 0 No optional inputs are being used. Default values +C will be used in all cases. +C 1 One or more optional inputs are being used. +C +C RWORK A real working array (single precision). The length of +C RWORK must be at least +C +C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM +C +C where +C NYH = the initial value of NEQ, +C MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a +C smaller value is given as an optional input), +C LWM = 0 if MITER = 0, +C LWM = NEQ**2 + 2 if MITER = 1 or 2, +C LWM = NEQ + 2 if MITER = 3, and +C LWM = (2*ML + MU + 1)*NEQ + 2 +C if MITER = 4 or 5. +C (See the MF description below for METH and MITER.) +C +C Thus if MAXORD has its default value and NEQ is constant, +C this length is: +C 20 + 16*NEQ for MF = 10, +C 22 + 16*NEQ + NEQ**2 for MF = 11 or 12, +C 22 + 17*NEQ for MF = 13, +C 22 + 17*NEQ + (2*ML + MU)*NEQ for MF = 14 or 15, +C 20 + 9*NEQ for MF = 20, +C 22 + 9*NEQ + NEQ**2 for MF = 21 or 22, +C 22 + 10*NEQ for MF = 23, +C 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25. +C +C The first 20 words of RWORK are reserved for conditional +C and optional inputs and optional outputs. +C +C The following word in RWORK is a conditional input: +C RWORK(1) = TCRIT, the critical value of t which the +C solver is not to overshoot. Required if ITASK +C is 4 or 5, and ignored otherwise. See ITASK. +C +C LRW The length of the array RWORK, as declared by the user. +C (This will be checked by the solver.) +C +C IWORK An integer work array. Its length must be at least +C 20 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or +C 20 + NEQ otherwise (MF = 11, 12, 14, 15, 21, 22, 24, 25). +C (See the MF description below for MITER.) The first few +C words of IWORK are used for conditional and optional +C inputs and optional outputs. +C +C The following two words in IWORK are conditional inputs: +C IWORK(1) = ML These are the lower and upper half- +C IWORK(2) = MU bandwidths, respectively, of the banded +C Jacobian, excluding the main diagonal. +C The band is defined by the matrix locations +C (i,j) with i - ML <= j <= i + MU. ML and MU +C must satisfy 0 <= ML,MU <= NEQ - 1. These are +C required if MITER is 4 or 5, and ignored +C otherwise. ML and MU may in fact be the band +C parameters for a matrix to which df/dy is only +C approximately equal. +C +C LIW The length of the array IWORK, as declared by the user. +C (This will be checked by the solver.) +C +C Note: The work arrays must not be altered between calls to SLSODE +C for the same problem, except possibly for the conditional and +C optional inputs, and except for the last 3*NEQ words of RWORK. +C The latter space is used for internal scratch space, and so is +C available for use by the user outside SLSODE between calls, if +C desired (but not for use by F or JAC). +C +C JAC The name of the user-supplied routine (MITER = 1 or 4) to +C compute the Jacobian matrix, df/dy, as a function of the +C scalar t and the vector y. (See the MF description below +C for MITER.) It is to have the form +C +C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) +C REAL T, Y(*), PD(NROWPD,*) +C +C where NEQ, T, Y, ML, MU, and NROWPD are input and the +C array PD is to be loaded with partial derivatives +C (elements of the Jacobian matrix) on output. PD must be +C given a first dimension of NROWPD. T and Y have the same +C meaning as in subroutine F. +C +C In the full matrix case (MITER = 1), ML and MU are +C ignored, and the Jacobian is to be loaded into PD in +C columnwise manner, with df(i)/dy(j) loaded into PD(i,j). +C +C In the band matrix case (MITER = 4), the elements within +C the band are to be loaded into PD in columnwise manner, +C with diagonal lines of df/dy loaded into the rows of PD. +C Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML +C and MU are the half-bandwidth parameters (see IWORK). +C The locations in PD in the two triangular areas which +C correspond to nonexistent matrix elements can be ignored +C or loaded arbitrarily, as they are overwritten by SLSODE. +C +C JAC need not provide df/dy exactly. A crude approximation +C (possibly with a smaller bandwidth) will do. +C +C In either case, PD is preset to zero by the solver, so +C that only the nonzero elements need be loaded by JAC. +C Each call to JAC is preceded by a call to F with the same +C arguments NEQ, T, and Y. Thus to gain some efficiency, +C intermediate quantities shared by both calculations may +C be saved in a user COMMON block by F and not recomputed +C by JAC, if desired. Also, JAC may alter the Y array, if +C desired. JAC must be declared EXTERNAL in the calling +C program. +C +C Subroutine JAC may access user-defined quantities in +C NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array +C (dimensioned in JAC) and/or Y has length exceeding +C NEQ(1). See the descriptions of NEQ and Y above. +C +C MF The method flag. Used only for input. The legal values +C of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, +C and 25. MF has decimal digits METH and MITER: +C MF = 10*METH + MITER . +C +C METH indicates the basic linear multistep method: +C 1 Implicit Adams method. +C 2 Method based on backward differentiation formulas +C (BDF's). +C +C MITER indicates the corrector iteration method: +C 0 Functional iteration (no Jacobian matrix is +C involved). +C 1 Chord iteration with a user-supplied full (NEQ by +C NEQ) Jacobian. +C 2 Chord iteration with an internally generated +C (difference quotient) full Jacobian (using NEQ +C extra calls to F per df/dy value). +C 3 Chord iteration with an internally generated +C diagonal Jacobian approximation (using one extra call +C to F per df/dy evaluation). +C 4 Chord iteration with a user-supplied banded Jacobian. +C 5 Chord iteration with an internally generated banded +C Jacobian (using ML + MU + 1 extra calls to F per +C df/dy evaluation). +C +C If MITER = 1 or 4, the user must supply a subroutine JAC +C (the name is arbitrary) as described above under JAC. +C For other values of MITER, a dummy argument 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 case +C all of these inputs are examined. A value of zero for any of +C 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, +C and then set those of interest to nonzero values. +C +C Name Location Meaning and default value +C ------ --------- ----------------------------------------------- +C H0 RWORK(5) Step size to be attempted on the first step. +C The default value is determined by the solver. +C HMAX RWORK(6) Maximum absolute step size allowed. The +C default value is infinite. +C HMIN RWORK(7) Minimum absolute step size allowed. The +C default value is 0. (This lower bound is not +C enforced on the final step before reaching +C TCRIT when ITASK = 4 or 5.) +C MAXORD IWORK(5) Maximum order to be allowed. The default value +C is 12 if METH = 1, and 5 if METH = 2. (See the +C MF description above for METH.) If MAXORD +C exceeds the default value, it will be reduced +C to the default value. If MAXORD is changed +C during the problem, it may cause the current +C order to be reduced. +C MXSTEP IWORK(6) Maximum number of (internally defined) steps +C allowed during one call to the solver. The +C default value is 500. +C MXHNIL IWORK(7) Maximum number of messages printed (per +C problem) warning that T + H = T on a step +C (H = step size). This must be positive to +C result in a nondefault value. The default +C value is 10. +C +C Optional Outputs +C ---------------- +C As optional additional output from SLSODE, the variables listed +C below are quantities related to the performance of SLSODE which +C are available to the user. These are communicated by way of the +C work arrays, but also have internal mnemonic names as shown. +C Except where stated otherwise, all of these outputs are defined on +C any successful return from SLSODE, and on any return with ISTATE = +C -1, -2, -4, -5, or -6. On an illegal input return (ISTATE = -3), +C they will be unchanged from their existing values (if any), except +C possibly for TOLSF, LENRW, and LENIW. On any error return, +C outputs relevant to the error will be defined, as noted below. +C +C Name Location Meaning +C ----- --------- ------------------------------------------------ +C HU RWORK(11) Step size in t last used (successfully). +C HCUR RWORK(12) Step size to be attempted on the next step. +C TCUR RWORK(13) Current value of the independent variable which +C the solver has actually reached, i.e., the +C current internal mesh point in t. On output, +C TCUR will always be at least as far as the +C argument T, but may be farther (if interpolation +C was done). +C TOLSF RWORK(14) Tolerance scale factor, greater than 1.0, +C computed when a request for too much accuracy +C was detected (ISTATE = -3 if detected at the +C start of the problem, ISTATE = -2 otherwise). +C If ITOL is left unaltered but RTOL and ATOL are +C uniformly scaled up by a factor of TOLSF for the +C next call, then the solver is deemed likely to +C succeed. (The user may also ignore TOLSF and +C alter the tolerance parameters in any other way +C appropriate.) +C NST IWORK(11) Number of steps taken for the problem so far. +C NFE IWORK(12) Number of F evaluations for the problem so far. +C NJE IWORK(13) Number of Jacobian evaluations (and of matrix LU +C decompositions) for the problem so far. +C NQU IWORK(14) Method order last used (successfully). +C NQCUR IWORK(15) Order to be attempted on the next step. +C IMXER IWORK(16) 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 LENRW IWORK(17) Length of RWORK actually required. This is +C defined on normal returns and on an illegal +C input return for insufficient storage. +C LENIW IWORK(18) Length of IWORK actually required. This is +C defined on normal returns and on an illegal +C input return for insufficient storage. +C +C The following two arrays are segments of the RWORK array which may +C also be of interest to the user as optional outputs. For each +C array, the table below gives its internal name, its base address +C 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 of +C NEQ. For j = 0,1,...,NQCUR, column j + 1 of +C YH contains HCUR**j/factorial(j) times the jth +C derivative of the interpolating polynomial +C currently representing the solution, evaluated +C at t = TCUR. +C ACOR LENRW-NEQ+1 Array of size NEQ used for the accumulated +C corrections on each step, scaled on output to +C represent the estimated local error in Y on +C the last step. This is the vector e in the +C description of the error control. It is +C defined only on successful return from SLSODE. +C +C +C Part 2. Other Callable Routines +C -------------------------------- +C +C The following are optional calls which the user may make to gain +C additional capabilities in conjunction with SLSODE. +C +C Form of call Function +C ------------------------ ---------------------------------------- +C CALL XSETUN(LUN) Set the logical unit number, LUN, for +C output of messages from SLSODE, if the +C default is not desired. The default +C value of LUN is 6. This call may be made +C at any time and will take effect +C immediately. +C CALL XSETF(MFLAG) Set a flag to control the printing of +C messages by SLSODE. MFLAG = 0 means do +C not print. (Danger: this risks losing +C valuable information.) MFLAG = 1 means +C print (the default). This call may be +C made at any time and will take effect +C immediately. +C CALL SSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the +C internal COMMON blocks used by SLSODE +C (see Part 3 below). RSAV must be a +C real array of length 218 or more, and +C ISAV must be an integer array of length +C 37 or more. JOB = 1 means save COMMON +C into RSAV/ISAV. JOB = 2 means restore +C COMMON from same. SSRCOM is useful if +C one is interrupting a run and restarting +C later, or alternating between two or +C more problems solved with SLSODE. +C CALL SINTDY(,,,,,) Provide derivatives of y, of various +C (see below) orders, at a specified point t, if +C desired. It may be called only after a +C successful return from SLSODE. Detailed +C instructions follow. +C +C Detailed instructions for using SINTDY +C -------------------------------------- +C The form of the CALL is: +C +C CALL SINTDY (T, K, RWORK(21), NYH, DKY, IFLAG) +C +C The input parameters are: +C +C T Value of independent variable where answers are +C desired (normally the same as the T last returned by +C SLSODE). For valid results, T must lie between +C TCUR - HU and TCUR. (See "Optional Outputs" above +C for TCUR and HU.) +C K Integer order of the derivative desired. K must +C satisfy 0 <= K <= NQCUR, where NQCUR is the current +C order (see "Optional Outputs"). The capability +C corresponding to K = 0, i.e., computing y(t), is +C already provided by SLSODE directly. Since +C NQCUR >= 1, the first derivative dy/dt is always +C available with SINTDY. +C RWORK(21) The base address of the history array YH. +C NYH Column length of YH, equal to the initial value of NEQ. +C +C The output parameters are: +C +C DKY Real array of length NEQ containing the computed value +C of the Kth derivative of y(t). +C IFLAG Integer flag, returned as 0 if K and T were legal, +C -1 if K was illegal, and -2 if T was illegal. +C On an error return, a message is also written. +C +C +C Part 3. Common Blocks +C ---------------------- +C +C If SLSODE is to be used in an overlay situation, the user must +C declare, in the primary overlay, the variables in: +C (1) the call sequence to SLSODE, +C (2) the internal COMMON block /SLS001/, of length 255 +C (218 single precision words followed by 37 integer words). +C +C If SLSODE is used on a system in which the contents of internal +C COMMON blocks are not preserved between calls, the user should +C declare the above COMMON block in his main program to insure that +C its contents are preserved. +C +C If the solution of a given problem by SLSODE is to be interrupted +C and then later continued, as when restarting an interrupted run or +C alternating between two or more problems, the user should save, +C following the return from the last SLSODE call prior to the +C interruption, the contents of the call sequence variables and the +C internal COMMON block, and later restore these values before the +C next SLSODE call for that problem. In addition, if XSETUN and/or +C XSETF was called for non-default handling of error messages, then +C these calls must be repeated. To save and restore the COMMON +C block, use subroutine SSRCOM (see Part 2 above). +C +C +C Part 4. Optionally Replaceable Solver Routines +C ----------------------------------------------- +C +C Below are descriptions of two routines in the SLSODE package which +C relate to the measurement of errors. Either routine can be +C replaced by a user-supplied version, if desired. However, since +C such a replacement may have a major impact on performance, it +C should be done only when absolutely necessary, and only with great +C caution. (Note: The means by which the package version of a +C routine is superseded by the user's version may be system- +C dependent.) +C +C SEWSET +C ------ +C The following subroutine is called just before each internal +C integration step, and sets the array of error weights, EWT, as +C described under ITOL/RTOL/ATOL above: +C +C SUBROUTINE SEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT) +C +C where NEQ, ITOL, RTOL, and ATOL are as in the SLSODE call +C sequence, YCUR contains the current dependent variable vector, +C and EWT is the array of weights set by SEWSET. +C +C If the user supplies this subroutine, it must return in EWT(i) +C (i = 1,...,NEQ) a positive quantity suitable for comparing errors +C in Y(i) to. The EWT array returned by SEWSET is passed to the +C SVNORM routine (see below), and also used by SLSODE in the +C computation of the optional output IMXER, the diagonal Jacobian +C approximation, and the increments for difference quotient +C Jacobians. +C +C In the user-supplied version of SEWSET, it may be desirable to use +C the current values of derivatives of y. Derivatives up to order NQ +C are available from the history array YH, described above under +C optional outputs. In SEWSET, YH is identical to the YCUR array, +C extended to NQ + 1 columns with a column length of NYH and scale +C factors of H**j/factorial(j). On the first call for the problem, +C given by NST = 0, NQ is 1 and H is temporarily set to 1.0. +C NYH is the initial value of NEQ. The quantities NQ, H, and NST +C can be obtained by including in SEWSET the statements: +C REAL RLS +C COMMON /SLS001/ RLS(218),ILS(37) +C NQ = ILS(33) +C NST = ILS(34) +C H = RLS(212) +C Thus, for example, the current value of dy/dt can be obtained as +C YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is unnecessary +C when NST = 0). +C +C SVNORM +C ------ +C SVNORM is a real function routine which computes the weighted +C root-mean-square norm of a vector v: +C +C d = SVNORM (n, v, w) +C +C where: +C n = the length of the vector, +C v = real array of length n containing the vector, +C w = real array of length n containing weights, +C d = SQRT( (1/n) * sum(v(i)*w(i))**2 ). +C +C SVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where +C EWT is as set by subroutine SEWSET. +C +C If the user supplies this function, it should return a nonnegative +C value of SVNORM suitable for use in the error control in SLSODE. +C None of the arguments should be altered by SVNORM. For example, a +C user-supplied SVNORM routine might: +C - Substitute a max-norm of (v(i)*w(i)) for the rms-norm, or +C - Ignore some components of v in the norm, with the effect of +C suppressing the error control on those components of Y. +C --------------------------------------------------------------------- +C***ROUTINES CALLED SEWSET, SINTDY, RUMACH, SSTODE, SVNORM, XERRWV +C***COMMON BLOCKS SLS001 +C***REVISION HISTORY (YYYYMMDD) +C 19791129 DATE WRITTEN +C 19791213 Minor changes to declarations; DELP init. in STODE. +C 19800118 Treat NEQ as array; integer declarations added throughout; +C minor changes to prologue. +C 19800306 Corrected TESCO(1,NQP1) setting in CFODE. +C 19800519 Corrected access of YH on forced order reduction; +C numerous corrections to prologues and other comments. +C 19800617 In main driver, added loading of SQRT(UROUND) in RWORK; +C minor corrections to main prologue. +C 19800923 Added zero initialization of HU and NQU. +C 19801218 Revised XERRWV routine; minor corrections to main prologue. +C 19810401 Minor changes to comments and an error message. +C 19810814 Numerous revisions: replaced EWT by 1/EWT; used flags +C JCUR, ICF, IERPJ, IERSL between STODE and subordinates; +C added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF; +C reorganized returns from STODE; reorganized type decls.; +C fixed message length in XERRWV; changed default LUNIT to 6; +C changed Common lengths; changed comments throughout. +C 19870330 Major update by ACH: corrected comments throughout; +C removed TRET from Common; rewrote EWSET with 4 loops; +C fixed t test in INTDY; added Cray directives in STODE; +C in STODE, fixed DELP init. and logic around PJAC call; +C combined routines to save/restore Common; +C passed LEVEL = 0 in error message calls (except run abort). +C 19890426 Modified prologue to SLATEC/LDOC format. (FNF) +C 19890501 Many improvements to prologue. (FNF) +C 19890503 A few final corrections to prologue. (FNF) +C 19890504 Minor cosmetic changes. (FNF) +C 19890510 Corrected description of Y in Arguments section. (FNF) +C 19890517 Minor corrections to prologue. (FNF) +C 19920514 Updated with prologue edited 891025 by G. Shaw for manual. +C 19920515 Converted source lines to upper case. (FNF) +C 19920603 Revised XERRWV calls using mixed upper-lower case. (ACH) +C 19920616 Revised prologue comment regarding CFT. (ACH) +C 19921116 Revised prologue comments regarding Common. (ACH). +C 19930326 Added comment about non-reentrancy. (FNF) +C 19930723 Changed R1MACH to RUMACH. (FNF) +C 19930801 Removed ILLIN and NTREP from Common (affects driver logic); +C minor changes to prologue and internal comments; +C changed Hollerith strings to quoted strings; +C changed internal comments to mixed case; +C replaced XERRWV with new version using character type; +C changed dummy dimensions from 1 to *. (ACH) +C 19930809 Changed to generic intrinsic names; changed names of +C subprograms and Common blocks to SLSODE etc. (ACH) +C 19930929 Eliminated use of REAL intrinsic; other minor changes. (ACH) +C 20010412 Removed all 'own' variables from Common block /SLS001/ +C (affects declarations in 6 routines). (ACH) +C 20010509 Minor corrections to prologue. (ACH) +C 20031105 Restored 'own' variables to Common block /SLS001/, to +C enable interrupt/restart feature. (ACH) +C 20031112 Added SAVE statements for data-loaded constants. +C +C*** END PROLOGUE SLSODE +C +C*Internal Notes: +C +C Other Routines in the SLSODE Package. +C +C In addition to Subroutine SLSODE, the SLSODE package includes the +C following subroutines and function routines: +C SINTDY computes an interpolated value of the y vector at t = TOUT. +C SSTODE is the core integrator, which does one step of the +C integration and the associated error control. +C SCFODE sets all method coefficients and test constants. +C SPREPJ computes and preprocesses the Jacobian matrix J = df/dy +C and the Newton iteration matrix P = I - h*l0*J. +C SSOLSY manages solution of linear system in chord iteration. +C SEWSET sets the error weight vector EWT before each step. +C SVNORM computes the weighted R.M.S. norm of a vector. +C SSRCOM is a user-callable routine to save and restore +C the contents of the internal Common block. +C DGETRF AND DGETRS ARE ROUTINES FROM LAPACK FOR SOLVING FULL +C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS. +C DGBTRF AND DGBTRS ARE ROUTINES FROM LAPACK FOR SOLVING BANDED +C LINEAR SYSTEMS. +C RUMACH computes the unit roundoff in a machine-independent manner. +C XERRWV, XSETUN, XSETF, IXSAV, IUMACH handle the printing of all +C error messages and warnings. XERRWV is machine-dependent. +C Note: SVNORM, RUMACH, IXSAV, and IUMACH are function routines. +C All the others are subroutines. +C +C**End +C +C Declare externals. + EXTERNAL SPREPJ, SSOLSY + REAL RUMACH, SVNORM +C +C Declare all other variables. + INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0, + 1 LENIW, LENRW, LENWM, ML, MORD, MU, MXHNL0, MXSTP0 + REAL ROWNS, + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + REAL ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI, + 1 TCRIT, TDIST, TNEXT, TOL, TOLSF, TP, SIZE, SUM, W0 + DIMENSION MORD(2) + LOGICAL IHIT + CHARACTER*80 MSG + SAVE MORD, MXSTP0, MXHNL0 +C----------------------------------------------------------------------- +C The following internal Common block contains +C (a) variables which are local to any subroutine but whose values must +C be preserved between calls to the routine ("own" variables), and +C (b) variables which are communicated between subroutines. +C The block SLS001 is declared in subroutines SLSODE, SINTDY, SSTODE, +C SPREPJ, and SSOLSY. +C Groups of variables are replaced by dummy arrays in the Common +C declarations in routines where those variables are not used. +C----------------------------------------------------------------------- + COMMON /SLS001/ ROWNS(209), + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 2 INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS(6), + 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU +C + DATA MORD(1),MORD(2)/12,5/, MXSTP0/500/, MXHNL0/10/ +C----------------------------------------------------------------------- +C Block A. +C This code block is executed on every call. +C It tests ISTATE and ITASK for legality and branches appropriately. +C If ISTATE .GT. 1 but the flag INIT shows that initialization has +C not yet been done, an error return occurs. +C If ISTATE = 1 and TOUT = T, return immediately. +C----------------------------------------------------------------------- +C +C***FIRST EXECUTABLE STATEMENT SLSODE + IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601 + IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602 + IF (ISTATE .EQ. 1) GO TO 10 + IF (INIT .EQ. 0) GO TO 603 + IF (ISTATE .EQ. 2) GO TO 200 + GO TO 20 + 10 INIT = 0 + IF (TOUT .EQ. T) RETURN +C----------------------------------------------------------------------- +C Block B. +C The next code block is executed for the initial call (ISTATE = 1), +C or for a continuation call with parameter changes (ISTATE = 3). +C It contains checking of all inputs and various initializations. +C +C First check legality of the non-optional inputs NEQ, ITOL, IOPT, +C MF, ML, and MU. +C----------------------------------------------------------------------- + 20 IF (NEQ(1) .LE. 0) GO TO 604 + IF (ISTATE .EQ. 1) GO TO 25 + IF (NEQ(1) .GT. N) GO TO 605 + 25 N = NEQ(1) + IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606 + IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 607 + METH = MF/10 + MITER = MF - 10*METH + IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608 + IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608 + IF (MITER .LE. 3) GO TO 30 + ML = IWORK(1) + MU = IWORK(2) + IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609 + IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610 + 30 CONTINUE +C Next process and check the optional inputs. -------------------------- + IF (IOPT .EQ. 1) GO TO 40 + MAXORD = MORD(METH) + MXSTEP = MXSTP0 + MXHNIL = MXHNL0 + IF (ISTATE .EQ. 1) H0 = 0.0E0 + HMXI = 0.0E0 + HMIN = 0.0E0 + GO TO 60 + 40 MAXORD = IWORK(5) + IF (MAXORD .LT. 0) GO TO 611 + IF (MAXORD .EQ. 0) MAXORD = 100 + MAXORD = MIN(MAXORD,MORD(METH)) + MXSTEP = IWORK(6) + IF (MXSTEP .LT. 0) GO TO 612 + IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0 + MXHNIL = IWORK(7) + IF (MXHNIL .LT. 0) GO TO 613 + IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0 + IF (ISTATE .NE. 1) GO TO 50 + H0 = RWORK(5) + IF ((TOUT - T)*H0 .LT. 0.0E0) GO TO 614 + 50 HMAX = RWORK(6) + IF (HMAX .LT. 0.0E0) GO TO 615 + HMXI = 0.0E0 + IF (HMAX .GT. 0.0E0) HMXI = 1.0E0/HMAX + HMIN = RWORK(7) + IF (HMIN .LT. 0.0E0) GO TO 616 +C----------------------------------------------------------------------- +C Set work array pointers and check lengths LRW and LIW. +C Pointers to segments of RWORK and IWORK are named by prefixing L to +C the name of the segment. E.g., the segment YH starts at RWORK(LYH). +C Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR. +C----------------------------------------------------------------------- + 60 LYH = 21 + IF (ISTATE .EQ. 1) NYH = N + LWM = LYH + (MAXORD + 1)*NYH + IF (MITER .EQ. 0) LENWM = 0 + IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2 + IF (MITER .EQ. 3) LENWM = N + 2 + IF (MITER .GE. 4) LENWM = (2*ML + MU + 1)*N + 2 + LEWT = LWM + LENWM + LSAVF = LEWT + N + LACOR = LSAVF + N + LENRW = LACOR + N - 1 + IWORK(17) = LENRW + LIWM = 1 + LENIW = 20 + N + IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20 + IWORK(18) = LENIW + IF (LENRW .GT. LRW) GO TO 617 + IF (LENIW .GT. LIW) GO TO 618 +C Check RTOL and ATOL for legality. ------------------------------------ + RTOLI = RTOL(1) + ATOLI = ATOL(1) + DO 70 I = 1,N + IF (ITOL .GE. 3) RTOLI = RTOL(I) + IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) + IF (RTOLI .LT. 0.0E0) GO TO 619 + IF (ATOLI .LT. 0.0E0) GO TO 620 + 70 CONTINUE + IF (ISTATE .EQ. 1) GO TO 100 +C If ISTATE = 3, set flag to signal parameter changes to SSTODE. ------- + JSTART = -1 + IF (NQ .LE. MAXORD) GO TO 90 +C MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. --------- + DO 80 I = 1,N + 80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1) +C Reload WM(1) = RWORK(LWM), since LWM may have changed. --------------- + 90 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) + IF (N .EQ. NYH) GO TO 200 +C NEQ was reduced. Zero part of YH to avoid undefined references. ----- + I1 = LYH + L*NYH + I2 = LYH + (MAXORD + 1)*NYH - 1 + IF (I1 .GT. I2) GO TO 200 + DO 95 I = I1,I2 + 95 RWORK(I) = 0.0E0 + GO TO 200 +C----------------------------------------------------------------------- +C Block C. +C The next block is for the initial call only (ISTATE = 1). +C It contains all remaining initializations, the initial call to F, +C and the calculation of the initial step size. +C The error weights in EWT are inverted after being loaded. +C----------------------------------------------------------------------- + 100 UROUND = RUMACH() + TN = T + IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 110 + TCRIT = RWORK(1) + IF ((TCRIT - TOUT)*(TOUT - T) .LT. 0.0E0) GO TO 625 + IF (H0 .NE. 0.0E0 .AND. (T + H0 - TCRIT)*H0 .GT. 0.0E0) + 1 H0 = TCRIT - T + 110 JSTART = 0 + IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) + NHNIL = 0 + NST = 0 + NJE = 0 + NSLAST = 0 + HU = 0.0E0 + NQU = 0 + CCMAX = 0.3E0 + MAXCOR = 3 + MSBP = 20 + MXNCF = 10 +C Initial call to F. (LF0 points to YH(*,2).) ------------------------- + LF0 = LYH + NYH + CALL F (NEQ, T, Y, RWORK(LF0)) + NFE = 1 +C Load the initial value vector in YH. --------------------------------- + DO 115 I = 1,N + 115 RWORK(I+LYH-1) = Y(I) +C Load and invert the EWT array. (H is temporarily set to 1.0.) ------- + NQ = 1 + H = 1.0E0 + CALL SEWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) + DO 120 I = 1,N + IF (RWORK(I+LEWT-1) .LE. 0.0E0) GO TO 621 + 120 RWORK(I+LEWT-1) = 1.0E0/RWORK(I+LEWT-1) +C----------------------------------------------------------------------- +C The coding below computes the step size, H0, to be attempted on the +C first step, unless the user has supplied a value for this. +C First check that TOUT - T differs significantly from zero. +C A scalar tolerance quantity TOL is computed, as MAX(RTOL(I)) +C if this is positive, or MAX(ATOL(I)/ABS(Y(I))) otherwise, adjusted +C so as to be between 100*UROUND and 1.0E-3. +C Then the computed value H0 is given by.. +C NEQ +C H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 ) +C 1 +C where w0 = MAX ( ABS(T), ABS(TOUT) ), +C f(i) = i-th component of initial value of f, +C ywt(i) = EWT(i)/TOL (a weight for y(i)). +C The sign of H0 is inferred from the initial values of TOUT and T. +C----------------------------------------------------------------------- + IF (H0 .NE. 0.0E0) GO TO 180 + TDIST = ABS(TOUT - T) + W0 = MAX(ABS(T),ABS(TOUT)) + IF (TDIST .LT. 2.0E0*UROUND*W0) GO TO 622 + TOL = RTOL(1) + IF (ITOL .LE. 2) GO TO 140 + DO 130 I = 1,N + 130 TOL = MAX(TOL,RTOL(I)) + 140 IF (TOL .GT. 0.0E0) GO TO 160 + ATOLI = ATOL(1) + DO 150 I = 1,N + IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) + AYI = ABS(Y(I)) + IF (AYI .NE. 0.0E0) TOL = MAX(TOL,ATOLI/AYI) + 150 CONTINUE + 160 TOL = MAX(TOL,100.0E0*UROUND) + TOL = MIN(TOL,0.001E0) + SUM = SVNORM (N, RWORK(LF0), RWORK(LEWT)) + SUM = 1.0E0/(TOL*W0*W0) + TOL*SUM**2 + H0 = 1.0E0/SQRT(SUM) + H0 = MIN(H0,TDIST) + H0 = SIGN(H0,TOUT-T) +C Adjust H0 if necessary to meet HMAX bound. --------------------------- + 180 RH = ABS(H0)*HMXI + IF (RH .GT. 1.0E0) H0 = H0/RH +C Load H with H0 and scale YH(*,2) by H0. ------------------------------ + H = H0 + DO 190 I = 1,N + 190 RWORK(I+LF0-1) = H0*RWORK(I+LF0-1) + GO TO 270 +C----------------------------------------------------------------------- +C Block D. +C The next code block is for continuation calls only (ISTATE = 2 or 3) +C and is to check stop conditions before taking a step. +C----------------------------------------------------------------------- + 200 NSLAST = NST + GO TO (210, 250, 220, 230, 240), ITASK + 210 IF ((TN - TOUT)*H .LT. 0.0E0) GO TO 250 + CALL SINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) + IF (IFLAG .NE. 0) GO TO 627 + T = TOUT + GO TO 420 + 220 TP = TN - HU*(1.0E0 + 100.0E0*UROUND) + IF ((TP - TOUT)*H .GT. 0.0E0) GO TO 623 + IF ((TN - TOUT)*H .LT. 0.0E0) GO TO 250 + GO TO 400 + 230 TCRIT = RWORK(1) + IF ((TN - TCRIT)*H .GT. 0.0E0) GO TO 624 + IF ((TCRIT - TOUT)*H .LT. 0.0E0) GO TO 625 + IF ((TN - TOUT)*H .LT. 0.0E0) GO TO 245 + CALL SINTDY (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.0E0) GO TO 624 + 245 HMX = ABS(TN) + ABS(H) + IHIT = ABS(TN - TCRIT) .LE. 100.0E0*UROUND*HMX + IF (IHIT) GO TO 400 + TNEXT = TN + H*(1.0E0 + 4.0E0*UROUND) + IF ((TNEXT - TCRIT)*H .LE. 0.0E0) GO TO 250 + H = (TCRIT - TN)*(1.0E0 - 4.0E0*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 SSTODE. +C +C This is a looping point for the integration steps. +C +C First check for too many steps being taken, update EWT (if not at +C start of problem), check for too much accuracy being requested, and +C check for H below the roundoff level in T. +C----------------------------------------------------------------------- + 250 CONTINUE + IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500 + CALL SEWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) + DO 260 I = 1,N + IF (RWORK(I+LEWT-1) .LE. 0.0E0) GO TO 510 + 260 RWORK(I+LEWT-1) = 1.0E0/RWORK(I+LEWT-1) + 270 TOLSF = UROUND*SVNORM (N, RWORK(LYH), RWORK(LEWT)) + IF (TOLSF .LE. 1.0E0) GO TO 280 + TOLSF = TOLSF*2.0E0 + 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 = 'SLSODE- Warning..internal T (=R1) and H (=R2) are' + CALL XERRWV (MSG, 50, 101, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG=' such that in the machine, T + H = T on the next step ' + CALL XERRWV (MSG, 60, 101, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' (H = step size). Solver will continue anyway' + CALL XERRWV (MSG, 50, 101, 0, 0, 0, 0, 2, TN, H) + IF (NHNIL .LT. MXHNIL) GO TO 290 + MSG = 'SLSODE- Above warning has been issued I1 times. ' + CALL XERRWV (MSG, 50, 102, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' It will not be issued again for this problem' + CALL XERRWV (MSG, 50, 102, 0, 1, MXHNIL, 0, 0, 0.0E0, 0.0E0) + 290 CONTINUE +C----------------------------------------------------------------------- +C CALL SSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,SPREPJ,SSOLSY) +C----------------------------------------------------------------------- + CALL SSTODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT), + 1 RWORK(LSAVF), RWORK(LACOR), RWORK(LWM), IWORK(LIWM), + 2 F, JAC, SPREPJ, SSOLSY) + KGO = 1 - KFLAG + GO TO (300, 530, 540), 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.0E0) GO TO 250 + CALL SINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) + T = TOUT + GO TO 420 +C ITASK = 3. Jump to exit if TOUT was reached. ------------------------ + 330 IF ((TN - TOUT)*H .GE. 0.0E0) GO TO 400 + GO TO 250 +C ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary. + 340 IF ((TN - TOUT)*H .LT. 0.0E0) GO TO 345 + CALL SINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) + T = TOUT + GO TO 420 + 345 HMX = ABS(TN) + ABS(H) + IHIT = ABS(TN - TCRIT) .LE. 100.0E0*UROUND*HMX + IF (IHIT) GO TO 400 + TNEXT = TN + H*(1.0E0 + 4.0E0*UROUND) + IF ((TNEXT - TCRIT)*H .LE. 0.0E0) GO TO 250 + H = (TCRIT - TN)*(1.0E0 - 4.0E0*UROUND) + JSTART = -2 + GO TO 250 +C ITASK = 5. See if TCRIT was reached and jump to exit. --------------- + 350 HMX = ABS(TN) + ABS(H) + IHIT = ABS(TN - TCRIT) .LE. 100.0E0*UROUND*HMX +C----------------------------------------------------------------------- +C Block G. +C The following block handles all successful returns from SLSODE. +C If ITASK .NE. 1, Y is loaded from YH and T is set accordingly. +C ISTATE is set to 2, and the optional outputs are loaded into the +C work arrays before returning. +C----------------------------------------------------------------------- + 400 DO 410 I = 1,N + 410 Y(I) = RWORK(I+LYH-1) + T = TN + IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420 + IF (IHIT) T = TCRIT + 420 ISTATE = 2 + RWORK(11) = HU + RWORK(12) = H + RWORK(13) = TN + IWORK(11) = NST + IWORK(12) = NFE + IWORK(13) = NJE + IWORK(14) = NQU + IWORK(15) = NQ + RETURN +C----------------------------------------------------------------------- +C Block H. +C The following block handles all unsuccessful returns other than +C those for illegal input. First the error message routine is called. +C If there was an error test or convergence test failure, IMXER is set. +C Then Y is loaded from YH and T is set to TN. The optional outputs +C are loaded into the work arrays before returning. +C----------------------------------------------------------------------- +C The maximum number of steps was taken before reaching TOUT. ---------- + 500 MSG = 'SLSODE- At current T (=R1), MXSTEP (=I1) steps ' + CALL XERRWV (MSG, 50, 201, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' taken on this call before reaching TOUT ' + CALL XERRWV (MSG, 50, 201, 0, 1, MXSTEP, 0, 1, TN, 0.0E0) + ISTATE = -1 + GO TO 580 +C EWT(I) .LE. 0.0 for some I (not at start of problem). ---------------- + 510 EWTI = RWORK(LEWT+I-1) + MSG = 'SLSODE- At T (=R1), EWT(I1) has become R2 .LE. 0.' + CALL XERRWV (MSG, 50, 202, 0, 1, I, 0, 2, TN, EWTI) + ISTATE = -6 + GO TO 580 +C Too much accuracy requested for machine precision. ------------------- + 520 MSG = 'SLSODE- At T (=R1), too much accuracy requested ' + CALL XERRWV (MSG, 50, 203, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' for precision of machine.. see TOLSF (=R2) ' + CALL XERRWV (MSG, 50, 203, 0, 0, 0, 0, 2, TN, TOLSF) + RWORK(14) = TOLSF + ISTATE = -2 + GO TO 580 +C KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. ----- + 530 MSG = 'SLSODE- At T(=R1) and step size H(=R2), the error' + CALL XERRWV (MSG, 50, 204, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' test failed repeatedly or with ABS(H) = HMIN' + CALL XERRWV (MSG, 50, 204, 0, 0, 0, 0, 2, TN, H) + ISTATE = -4 + GO TO 560 +C KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ---- + 540 MSG = 'SLSODE- At T (=R1) and step size H (=R2), the ' + CALL XERRWV (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' corrector convergence failed repeatedly ' + CALL XERRWV (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG = ' or with ABS(H) = HMIN ' + CALL XERRWV (MSG, 30, 205, 0, 0, 0, 0, 2, TN, H) + ISTATE = -5 +C Compute IMXER if relevant. ------------------------------------------- + 560 BIG = 0.0E0 + IMXER = 1 + DO 570 I = 1,N + SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1)) + IF (BIG .GE. SIZE) GO TO 570 + BIG = SIZE + IMXER = I + 570 CONTINUE + IWORK(16) = IMXER +C Set Y vector, T, and optional outputs. ------------------------------- + 580 DO 590 I = 1,N + 590 Y(I) = RWORK(I+LYH-1) + T = TN + RWORK(11) = HU + RWORK(12) = H + RWORK(13) = TN + IWORK(11) = NST + IWORK(12) = NFE + IWORK(13) = NJE + IWORK(14) = NQU + IWORK(15) = NQ + RETURN +C----------------------------------------------------------------------- +C Block I. +C The following block handles all error returns due to illegal input +C (ISTATE = -3), as detected before calling the core integrator. +C First the error message routine is called. If the illegal input +C is a negative ISTATE, the run is aborted (apparent infinite loop). +C----------------------------------------------------------------------- + 601 MSG = 'SLSODE- ISTATE (=I1) illegal ' + CALL XERRWV (MSG, 30, 1, 0, 1, ISTATE, 0, 0, 0.0E0, 0.0E0) + IF (ISTATE .LT. 0) GO TO 800 + GO TO 700 + 602 MSG = 'SLSODE- ITASK (=I1) illegal ' + CALL XERRWV (MSG, 30, 2, 0, 1, ITASK, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 603 MSG = 'SLSODE- ISTATE .GT. 1 but SLSODE not initialized ' + CALL XERRWV (MSG, 50, 3, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 604 MSG = 'SLSODE- NEQ (=I1) .LT. 1 ' + CALL XERRWV (MSG, 30, 4, 0, 1, NEQ(1), 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 605 MSG = 'SLSODE- ISTATE = 3 and NEQ increased (I1 to I2) ' + CALL XERRWV (MSG, 50, 5, 0, 2, N, NEQ(1), 0, 0.0E0, 0.0E0) + GO TO 700 + 606 MSG = 'SLSODE- ITOL (=I1) illegal ' + CALL XERRWV (MSG, 30, 6, 0, 1, ITOL, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 607 MSG = 'SLSODE- IOPT (=I1) illegal ' + CALL XERRWV (MSG, 30, 7, 0, 1, IOPT, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 608 MSG = 'SLSODE- MF (=I1) illegal ' + CALL XERRWV (MSG, 30, 8, 0, 1, MF, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 609 MSG = 'SLSODE- ML (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)' + CALL XERRWV (MSG, 50, 9, 0, 2, ML, NEQ(1), 0, 0.0E0, 0.0E0) + GO TO 700 + 610 MSG = 'SLSODE- MU (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)' + CALL XERRWV (MSG, 50, 10, 0, 2, MU, NEQ(1), 0, 0.0E0, 0.0E0) + GO TO 700 + 611 MSG = 'SLSODE- MAXORD (=I1) .LT. 0 ' + CALL XERRWV (MSG, 30, 11, 0, 1, MAXORD, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 612 MSG = 'SLSODE- MXSTEP (=I1) .LT. 0 ' + CALL XERRWV (MSG, 30, 12, 0, 1, MXSTEP, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 613 MSG = 'SLSODE- MXHNIL (=I1) .LT. 0 ' + CALL XERRWV (MSG, 30, 13, 0, 1, MXHNIL, 0, 0, 0.0E0, 0.0E0) + GO TO 700 + 614 MSG = 'SLSODE- TOUT (=R1) behind T (=R2) ' + CALL XERRWV (MSG, 40, 14, 0, 0, 0, 0, 2, TOUT, T) + MSG = ' Integration direction is given by H0 (=R1) ' + CALL XERRWV (MSG, 50, 14, 0, 0, 0, 0, 1, H0, 0.0E0) + GO TO 700 + 615 MSG = 'SLSODE- HMAX (=R1) .LT. 0.0 ' + CALL XERRWV (MSG, 30, 15, 0, 0, 0, 0, 1, HMAX, 0.0E0) + GO TO 700 + 616 MSG = 'SLSODE- HMIN (=R1) .LT. 0.0 ' + CALL XERRWV (MSG, 30, 16, 0, 0, 0, 0, 1, HMIN, 0.0E0) + GO TO 700 + 617 CONTINUE + MSG='SLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)' + CALL XERRWV (MSG, 60, 17, 0, 2, LENRW, LRW, 0, 0.0E0, 0.0E0) + GO TO 700 + 618 CONTINUE + MSG='SLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)' + CALL XERRWV (MSG, 60, 18, 0, 2, LENIW, LIW, 0, 0.0E0, 0.0E0) + GO TO 700 + 619 MSG = 'SLSODE- RTOL(I1) is R1 .LT. 0.0 ' + CALL XERRWV (MSG, 40, 19, 0, 1, I, 0, 1, RTOLI, 0.0E0) + GO TO 700 + 620 MSG = 'SLSODE- ATOL(I1) is R1 .LT. 0.0 ' + CALL XERRWV (MSG, 40, 20, 0, 1, I, 0, 1, ATOLI, 0.0E0) + GO TO 700 + 621 EWTI = RWORK(LEWT+I-1) + MSG = 'SLSODE- EWT(I1) is R1 .LE. 0.0 ' + CALL XERRWV (MSG, 40, 21, 0, 1, I, 0, 1, EWTI, 0.0E0) + GO TO 700 + 622 CONTINUE + MSG='SLSODE- TOUT (=R1) too close to T(=R2) to start integration' + CALL XERRWV (MSG, 60, 22, 0, 0, 0, 0, 2, TOUT, T) + GO TO 700 + 623 CONTINUE + MSG='SLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' + CALL XERRWV (MSG, 60, 23, 0, 1, ITASK, 0, 2, TOUT, TP) + GO TO 700 + 624 CONTINUE + MSG='SLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) ' + CALL XERRWV (MSG, 60, 24, 0, 0, 0, 0, 2, TCRIT, TN) + GO TO 700 + 625 CONTINUE + MSG='SLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' + CALL XERRWV (MSG, 60, 25, 0, 0, 0, 0, 2, TCRIT, TOUT) + GO TO 700 + 626 MSG = 'SLSODE- At start of problem, too much accuracy ' + CALL XERRWV (MSG, 50, 26, 0, 0, 0, 0, 0, 0.0E0, 0.0E0) + MSG=' requested for precision of machine.. See TOLSF (=R1) ' + CALL XERRWV (MSG, 60, 26, 0, 0, 0, 0, 1, TOLSF, 0.0E0) + RWORK(14) = TOLSF + GO TO 700 + 627 MSG = 'SLSODE- Trouble in SINTDY. ITASK = I1, TOUT = R1' + CALL XERRWV (MSG, 50, 27, 0, 1, ITASK, 0, 1, TOUT, 0.0E0) +C + 700 ISTATE = -3 + RETURN +C + 800 MSG = 'SLSODE- Run aborted.. apparent infinite loop ' + CALL XERRWV (MSG, 50, 303, 2, 0, 0, 0, 0, 0.0E0, 0.0E0) + RETURN +C----------------------- END OF SUBROUTINE SLSODE ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/sprepj.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,193 @@ + SUBROUTINE SPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, + 1 F, JAC) +C***BEGIN PROLOGUE SPREPJ +C***SUBSIDIARY +C***PURPOSE Compute and process Newton iteration matrix. +C***TYPE SINGLE PRECISION (SPREPJ-S, DPREPJ-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C SPREPJ is called by SSTODE to compute and process the matrix +C P = I - h*el(1)*J , where J is an approximation to the Jacobian. +C Here J is computed by the user-supplied routine JAC if +C MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5. +C If MITER = 3, a diagonal approximation to J is used. +C J is stored in WM and replaced by P. If MITER .ne. 3, P is then +C subjected to LU decomposition in preparation for later solution +C of linear systems with P as coefficient matrix. This is done +C by SGETRF if MITER = 1 or 2, and by SGBTRF if MITER = 4 or 5. +C +C In addition to variables described in SSTODE and SLSODE prologues, +C communication with SPREPJ uses the following: +C Y = array containing predicted values on entry. +C FTEM = work array of length N (ACOR in SSTODE). +C SAVF = array containing f evaluated at predicted y. +C WM = real work space for matrices. On output it contains the +C inverse diagonal matrix if MITER = 3 and the LU decomposition +C of P if MITER is 1, 2 , 4, or 5. +C Storage of matrix elements starts at WM(3). +C WM also contains the following matrix-related data: +C WM(1) = SQRT(UROUND), used in numerical Jacobian increments. +C WM(2) = H*EL0, saved for later use if MITER = 3. +C IWM = integer work space containing pivot information, starting at +C IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band +C parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. +C EL0 = EL(1) (input). +C IERPJ = output error flag, = 0 if no trouble, .gt. 0 if +C P matrix found to be singular. +C JCUR = output flag = 1 to indicate that the Jacobian matrix +C (or approximation) is now current. +C This routine also uses the COMMON variables EL0, H, TN, UROUND, +C MITER, N, NFE, and NJE. +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED SGBTRF, SGETRF, SVNORM +C***COMMON BLOCKS SLS001 +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890504 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C 010412 Reduced size of Common block /SLS001/. (ACH) +C 031105 Restored 'own' variables to Common block /SLS001/, to +C enable interrupt/restart feature. (ACH) +C***END PROLOGUE SPREPJ +C**End + EXTERNAL F, JAC + INTEGER NEQ, NYH, IWM + REAL Y, YH, EWT, FTEM, SAVF, WM + DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), + 1 WM(*), IWM(*) + INTEGER IOWND, IOWNS, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + REAL ROWNS, + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + COMMON /SLS001/ ROWNS(209), + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 2 IOWND(6), IOWNS(6), + 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, + 1 MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1 + REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, + 1 SVNORM +C +C***FIRST EXECUTABLE STATEMENT SPREPJ + NJE = NJE + 1 + IERPJ = 0 + JCUR = 1 + HL0 = H*EL0 + GO TO (100, 200, 300, 400, 500), MITER +C If MITER = 1, call JAC and multiply by scalar. ----------------------- + 100 LENP = N*N + DO 110 I = 1,LENP + 110 WM(I+2) = 0.0E0 + CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N) + CON = -HL0 + DO 120 I = 1,LENP + 120 WM(I+2) = WM(I+2)*CON + GO TO 240 +C If MITER = 2, make N calls to F to approximate J. -------------------- + 200 FAC = SVNORM (N, SAVF, EWT) + R0 = 1000.0E0*ABS(H)*UROUND*N*FAC + IF (R0 .EQ. 0.0E0) R0 = 1.0E0 + SRUR = WM(1) + J1 = 2 + DO 230 J = 1,N + YJ = Y(J) + R = MAX(SRUR*ABS(YJ),R0/EWT(J)) + Y(J) = Y(J) + R + FAC = -HL0/R + CALL F (NEQ, TN, Y, FTEM) + DO 220 I = 1,N + 220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC + Y(J) = YJ + J1 = J1 + N + 230 CONTINUE + NFE = NFE + N +C Add identity matrix. ------------------------------------------------- + 240 J = 3 + NP1 = N + 1 + DO 250 I = 1,N + WM(J) = WM(J) + 1.0E0 + 250 J = J + NP1 +C Do LU decomposition on P. -------------------------------------------- + CALL SGETRF (N, N, WM(3), N, IWM(21), IER) + IF (IER .NE. 0) IERPJ = 1 + RETURN +C If MITER = 3, construct a diagonal approximation to J and P. --------- + 300 WM(2) = HL0 + R = EL0*0.1E0 + DO 310 I = 1,N + 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) + CALL F (NEQ, TN, Y, WM(3)) + NFE = NFE + 1 + DO 320 I = 1,N + R0 = H*SAVF(I) - YH(I,2) + DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I)) + WM(I+2) = 1.0E0 + IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320 + IF (ABS(DI) .EQ. 0.0E0) GO TO 330 + WM(I+2) = 0.1E0*R0/DI + 320 CONTINUE + RETURN + 330 IERPJ = 1 + RETURN +C If MITER = 4, call JAC and multiply by scalar. ----------------------- + 400 ML = IWM(1) + MU = IWM(2) + ML3 = ML + 3 + MBAND = ML + MU + 1 + MEBAND = MBAND + ML + LENP = MEBAND*N + DO 410 I = 1,LENP + 410 WM(I+2) = 0.0E0 + CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND) + CON = -HL0 + DO 420 I = 1,LENP + 420 WM(I+2) = WM(I+2)*CON + GO TO 570 +C If MITER = 5, make MBAND calls to F to approximate J. ---------------- + 500 ML = IWM(1) + MU = IWM(2) + MBAND = ML + MU + 1 + MBA = MIN(MBAND,N) + MEBAND = MBAND + ML + MEB1 = MEBAND - 1 + SRUR = WM(1) + FAC = SVNORM (N, SAVF, EWT) + R0 = 1000.0E0*ABS(H)*UROUND*N*FAC + IF (R0 .EQ. 0.0E0) R0 = 1.0E0 + DO 560 J = 1,MBA + DO 530 I = J,N,MBAND + YI = Y(I) + R = MAX(SRUR*ABS(YI),R0/EWT(I)) + 530 Y(I) = Y(I) + R + CALL F (NEQ, TN, Y, FTEM) + DO 550 JJ = J,N,MBAND + Y(JJ) = YH(JJ,1) + YJJ = Y(JJ) + R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ)) + FAC = -HL0/R + I1 = MAX(JJ-MU,1) + I2 = MIN(JJ+ML,N) + II = JJ*MEB1 - ML + 2 + DO 540 I = I1,I2 + 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC + 550 CONTINUE + 560 CONTINUE + NFE = NFE + MBA +C Add identity matrix. ------------------------------------------------- + 570 II = MBAND + 2 + DO 580 I = 1,N + WM(II) = WM(II) + 1.0E0 + 580 II = II + MEBAND +C Do LU decomposition of P. -------------------------------------------- + CALL SGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER) + IF (IER .NE. 0) IERPJ = 1 + RETURN +C----------------------- END OF SUBROUTINE SPREPJ ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/ssolsy.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,91 @@ + SUBROUTINE SSOLSY (WM, IWM, X, TEM) +C***BEGIN PROLOGUE SSOLSY +C***SUBSIDIARY +C***PURPOSE ODEPACK linear system solver. +C***TYPE SINGLE PRECISION (SSOLSY-S, DSOLSY-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C This routine manages the solution of the linear system arising from +C a chord iteration. It is called if MITER .ne. 0. +C If MITER is 1 or 2, it calls SGETRF to accomplish this. +C If MITER = 3 it updates the coefficient h*EL0 in the diagonal +C matrix, and then computes the solution. +C If MITER is 4 or 5, it calls SGBTRS. +C Communication with SSOLSY uses the following variables: +C WM = real work space containing the inverse diagonal matrix if +C MITER = 3 and the LU decomposition of the matrix otherwise. +C Storage of matrix elements starts at WM(3). +C WM also contains the following matrix-related data: +C WM(1) = SQRT(UROUND) (not used here), +C WM(2) = HL0, the previous value of h*EL0, used if MITER = 3. +C IWM = integer work space containing pivot information, starting at +C IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band +C parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. +C X = the right-hand side vector on input, and the solution vector +C on output, of length N. +C TEM = vector of work space of length N, not used in this version. +C IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred. +C IERSL = 1 if a singular matrix arose with MITER = 3. +C This routine also uses the COMMON variables EL0, H, MITER, and N. +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED SGBTRS, SGETRS +C***COMMON BLOCKS SLS001 +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890503 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C 010412 Reduced size of Common block /SLS001/. (ACH) +C 031105 Restored 'own' variables to Common block /SLS001/, to +C enable interrupt/restart feature. (ACH) +C***END PROLOGUE SSOLSY +C**End + INTEGER IWM + REAL WM, X, TEM + DIMENSION WM(*), IWM(*), X(*), TEM(*) + INTEGER IOWND, IOWNS, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + REAL ROWNS, + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + COMMON /SLS001/ ROWNS(209), + 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 2 IOWND(6), IOWNS(6), + 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, MEBAND, ML, MU + REAL DI, HL0, PHL0, R +C +C***FIRST EXECUTABLE STATEMENT SSOLSY + IERSL = 0 + GO TO (100, 100, 300, 400, 400), MITER + 100 CALL SGETRS ( 'N', N, 1, WM(3), N, IWM(21), X, N, INLPCK) + RETURN +C + 300 PHL0 = WM(2) + HL0 = H*EL0 + WM(2) = HL0 + IF (HL0 .EQ. PHL0) GO TO 330 + R = HL0/PHL0 + DO 320 I = 1,N + DI = 1.0E0 - R*(1.0E0 - 1.0E0/WM(I+2)) + IF (ABS(DI) .EQ. 0.0E0) GO TO 390 + 320 WM(I+2) = 1.0E0/DI + 330 DO 340 I = 1,N + 340 X(I) = WM(I+2)*X(I) + RETURN + 390 IERSL = 1 + RETURN +C + 400 ML = IWM(1) + MU = IWM(2) + MEBAND = 2*ML + MU + 1 + CALL SGBTRS ( 'N', N, ML, MU, 1, WM(3), MEBAND, IWM(21), X, N, + * INLPCK) + RETURN +C----------------------- END OF SUBROUTINE SSOLSY ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/sstode.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,497 @@ + SUBROUTINE SSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, + 1 WM, IWM, F, JAC, PJAC, SLVS) +C***BEGIN PROLOGUE SSTODE +C***SUBSIDIARY +C***PURPOSE Performs one step of an ODEPACK integration. +C***TYPE SINGLE PRECISION (SSTODE-S, DSTODE-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C SSTODE performs one step of the integration of an initial value +C problem for a system of ordinary differential equations. +C Note: SSTODE is independent of the value of the iteration method +C indicator MITER, when this is .ne. 0, and hence is independent +C of the type of chord method used, or the Jacobian structure. +C Communication with SSTODE is done with the following variables: +C +C NEQ = integer array containing problem size in NEQ(1), and +C passed as the NEQ argument in all calls to F and JAC. +C Y = an array of length .ge. N used as the Y argument in +C all calls to F and JAC. +C YH = an NYH by LMAX array containing the dependent variables +C and their approximate scaled derivatives, where +C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate +C j-th derivative of y(i), scaled by h**j/factorial(j) +C (j = 0,1,...,NQ). on entry for the first step, the first +C two columns of YH must be set from the initial values. +C NYH = a constant integer .ge. N, the first dimension of YH. +C YH1 = a one-dimensional array occupying the same space as YH. +C EWT = an array of length N containing multiplicative weights +C for local error measurements. Local errors in Y(i) are +C compared to 1.0/EWT(i) in various error tests. +C SAVF = an array of working storage, of length N. +C Also used for input of YH(*,MAXORD+2) when JSTART = -1 +C and MAXORD .lt. the current order NQ. +C ACOR = a work array of length N, used for the accumulated +C corrections. On a successful return, ACOR(i) contains +C the estimated one-step local error in Y(i). +C WM,IWM = real and integer work arrays associated with matrix +C operations in chord iteration (MITER .ne. 0). +C PJAC = name of routine to evaluate and preprocess Jacobian matrix +C and P = I - h*el0*JAC, if a chord method is being used. +C SLVS = name of routine to solve linear system in chord iteration. +C CCMAX = maximum relative change in h*el0 before PJAC is called. +C H = the step size to be attempted on the next step. +C H is altered by the error control algorithm during the +C problem. H can be either positive or negative, but its +C sign must remain constant throughout the problem. +C HMIN = the minimum absolute value of the step size h to be used. +C HMXI = inverse of the maximum absolute value of h to be used. +C HMXI = 0.0 is allowed and corresponds to an infinite hmax. +C HMIN and HMXI may be changed at any time, but will not +C take effect until the next change of h is considered. +C TN = the independent variable. TN is updated on each step taken. +C JSTART = an integer used for input only, with the following +C values and meanings: +C 0 perform the first step. +C .gt.0 take a new step continuing from the last. +C -1 take the next step with a new value of H, MAXORD, +C N, METH, MITER, and/or matrix parameters. +C -2 take the next step with a new value of H, +C but with other inputs unchanged. +C On return, JSTART is set to 1 to facilitate continuation. +C KFLAG = a completion code with the following meanings: +C 0 the step was succesful. +C -1 the requested error could not be achieved. +C -2 corrector convergence could not be achieved. +C -3 fatal error in PJAC or SLVS. +C A return with KFLAG = -1 or -2 means either +C abs(H) = HMIN or 10 consecutive failures occurred. +C On a return with KFLAG negative, the values of TN and +C the YH array are as of the beginning of the last +C step, and H is the last step size attempted. +C MAXORD = the maximum order of integration method to be allowed. +C MAXCOR = the maximum number of corrector iterations allowed. +C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0). +C MXNCF = maximum number of convergence failures allowed. +C METH/MITER = the method flags. See description in driver. +C N = the number of first-order differential equations. +C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD, +C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON. +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED SCFODE, SVNORM +C***COMMON BLOCKS SLS001 +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890503 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C 010413 Reduced size of Common block /SLS001/. (ACH) +C 031105 Restored 'own' variables to Common block /SLS001/, to +C enable interrupt/restart feature. (ACH) +C***END PROLOGUE SSTODE +C**End + EXTERNAL F, JAC, PJAC, SLVS + INTEGER NEQ, NYH, IWM + REAL Y, YH, YH1, EWT, SAVF, ACOR, WM + DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), + 1 ACOR(*), WM(*), IWM(*) + INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, + 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ + REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, + 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND + REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, + 1 R, RH, RHDN, RHSM, RHUP, TOLD, SVNORM + COMMON /SLS001/ CONIT, CRATE, EL(13), ELCO(13,12), + 1 HOLD, RMAX, TESCO(3,12), + 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, + 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, + 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, + 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU +C +C***FIRST EXECUTABLE STATEMENT SSTODE + KFLAG = 0 + TOLD = TN + NCF = 0 + IERPJ = 0 + IERSL = 0 + JCUR = 0 + ICF = 0 + DELP = 0.0E0 + IF (JSTART .GT. 0) GO TO 200 + IF (JSTART .EQ. -1) GO TO 100 + IF (JSTART .EQ. -2) GO TO 160 +C----------------------------------------------------------------------- +C On the first call, the order is set to 1, and other variables are +C initialized. RMAX is the maximum ratio by which H can be increased +C in a single step. It is initially 1.E4 to compensate for the small +C initial H, but then is normally equal to 10. If a failure +C occurs (in corrector convergence or error test), RMAX is set to 2 +C for the next increase. +C----------------------------------------------------------------------- + LMAX = MAXORD + 1 + NQ = 1 + L = 2 + IALTH = 2 + RMAX = 10000.0E0 + RC = 0.0E0 + EL0 = 1.0E0 + CRATE = 0.7E0 + HOLD = H + MEO = METH + NSLP = 0 + IPUP = MITER + IRET = 3 + GO TO 140 +C----------------------------------------------------------------------- +C The following block handles preliminaries needed when JSTART = -1. +C IPUP is set to MITER to force a matrix update. +C If an order increase is about to be considered (IALTH = 1), +C IALTH is reset to 2 to postpone consideration one more step. +C If the caller has changed METH, SCFODE is called to reset +C the coefficients of the method. +C If the caller has changed MAXORD to a value less than the current +C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly. +C If H is to be changed, YH must be rescaled. +C If H or METH is being changed, IALTH is reset to L = NQ + 1 +C to prevent further changes in H for that many steps. +C----------------------------------------------------------------------- + 100 IPUP = MITER + LMAX = MAXORD + 1 + IF (IALTH .EQ. 1) IALTH = 2 + IF (METH .EQ. MEO) GO TO 110 + CALL SCFODE (METH, ELCO, TESCO) + MEO = METH + IF (NQ .GT. MAXORD) GO TO 120 + IALTH = L + IRET = 1 + GO TO 150 + 110 IF (NQ .LE. MAXORD) GO TO 160 + 120 NQ = MAXORD + L = LMAX + DO 125 I = 1,L + 125 EL(I) = ELCO(I,NQ) + NQNYH = NQ*NYH + RC = RC*EL(1)/EL0 + EL0 = EL(1) + CONIT = 0.5E0/(NQ+2) + DDN = SVNORM (N, SAVF, EWT)/TESCO(1,L) + EXDN = 1.0E0/L + RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) + RH = MIN(RHDN,1.0E0) + IREDO = 3 + IF (H .EQ. HOLD) GO TO 170 + RH = MIN(RH,ABS(H/HOLD)) + H = HOLD + GO TO 175 +C----------------------------------------------------------------------- +C SCFODE is called to get all the integration coefficients for the +C current METH. Then the EL vector and related constants are reset +C whenever the order NQ is changed, or at the start of the problem. +C----------------------------------------------------------------------- + 140 CALL SCFODE (METH, ELCO, TESCO) + 150 DO 155 I = 1,L + 155 EL(I) = ELCO(I,NQ) + NQNYH = NQ*NYH + RC = RC*EL(1)/EL0 + EL0 = EL(1) + CONIT = 0.5E0/(NQ+2) + GO TO (160, 170, 200), IRET +C----------------------------------------------------------------------- +C If H is being changed, the H ratio RH is checked against +C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to +C L = NQ + 1 to prevent a change of H for that many steps, unless +C forced by a convergence or error test failure. +C----------------------------------------------------------------------- + 160 IF (H .EQ. HOLD) GO TO 200 + RH = H/HOLD + H = HOLD + IREDO = 3 + GO TO 175 + 170 RH = MAX(RH,HMIN/ABS(H)) + 175 RH = MIN(RH,RMAX) + RH = RH/MAX(1.0E0,ABS(H)*HMXI*RH) + R = 1.0E0 + DO 180 J = 2,L + R = R*RH + DO 180 I = 1,N + 180 YH(I,J) = YH(I,J)*R + H = H*RH + RC = RC*RH + IALTH = L + IF (IREDO .EQ. 0) GO TO 690 +C----------------------------------------------------------------------- +C This section computes the predicted values by effectively +C multiplying the YH array by the Pascal Triangle matrix. +C RC is the ratio of new to old values of the coefficient H*EL(1). +C When RC differs from 1 by more than CCMAX, IPUP is set to MITER +C to force PJAC to be called, if a Jacobian is involved. +C In any case, PJAC is called at least every MSBP steps. +C----------------------------------------------------------------------- + 200 IF (ABS(RC-1.0E0) .GT. CCMAX) IPUP = MITER + IF (NST .GE. NSLP+MSBP) IPUP = MITER + TN = TN + H + I1 = NQNYH + 1 + DO 215 JB = 1,NQ + I1 = I1 - NYH +Cdir$ ivdep + DO 210 I = I1,NQNYH + 210 YH1(I) = YH1(I) + YH1(I+NYH) + 215 CONTINUE +C----------------------------------------------------------------------- +C Up to MAXCOR corrector iterations are taken. A convergence test is +C made on the R.M.S. norm of each correction, weighted by the error +C weight vector EWT. The sum of the corrections is accumulated in the +C vector ACOR(i). The YH array is not altered in the corrector loop. +C----------------------------------------------------------------------- + 220 M = 0 + DO 230 I = 1,N + 230 Y(I) = YH(I,1) + CALL F (NEQ, TN, Y, SAVF) + NFE = NFE + 1 + IF (IPUP .LE. 0) GO TO 250 +C----------------------------------------------------------------------- +C If indicated, the matrix P = I - h*el(1)*J is reevaluated and +C preprocessed before starting the corrector iteration. IPUP is set +C to 0 as an indicator that this has been done. +C----------------------------------------------------------------------- + CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC) + IPUP = 0 + RC = 1.0E0 + NSLP = NST + CRATE = 0.7E0 + IF (IERPJ .NE. 0) GO TO 430 + 250 DO 260 I = 1,N + 260 ACOR(I) = 0.0E0 + 270 IF (MITER .NE. 0) GO TO 350 +C----------------------------------------------------------------------- +C In the case of functional iteration, update Y directly from +C the result of the last function evaluation. +C----------------------------------------------------------------------- + DO 290 I = 1,N + SAVF(I) = H*SAVF(I) - YH(I,2) + 290 Y(I) = SAVF(I) - ACOR(I) + DEL = SVNORM (N, Y, EWT) + DO 300 I = 1,N + Y(I) = YH(I,1) + EL(1)*SAVF(I) + 300 ACOR(I) = SAVF(I) + GO TO 400 +C----------------------------------------------------------------------- +C In the case of the chord method, compute the corrector error, +C and solve the linear system with that as right-hand side and +C P as coefficient matrix. +C----------------------------------------------------------------------- + 350 DO 360 I = 1,N + 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) + CALL SLVS (WM, IWM, Y, SAVF) + IF (IERSL .LT. 0) GO TO 430 + IF (IERSL .GT. 0) GO TO 410 + DEL = SVNORM (N, Y, EWT) + DO 380 I = 1,N + ACOR(I) = ACOR(I) + Y(I) + 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) +C----------------------------------------------------------------------- +C Test for convergence. If M.gt.0, an estimate of the convergence +C rate constant is stored in CRATE, and this is used in the test. +C----------------------------------------------------------------------- + 400 IF (M .NE. 0) CRATE = MAX(0.2E0*CRATE,DEL/DELP) + DCON = DEL*MIN(1.0E0,1.5E0*CRATE)/(TESCO(2,NQ)*CONIT) + IF (DCON .LE. 1.0E0) GO TO 450 + M = M + 1 + IF (M .EQ. MAXCOR) GO TO 410 + IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410 + DELP = DEL + CALL F (NEQ, TN, Y, SAVF) + NFE = NFE + 1 + GO TO 270 +C----------------------------------------------------------------------- +C The corrector iteration failed to converge. +C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for +C the next try. Otherwise the YH array is retracted to its values +C before prediction, and H is reduced, if possible. If H cannot be +C reduced or MXNCF failures have occurred, exit with KFLAG = -2. +C----------------------------------------------------------------------- + 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 + ICF = 1 + IPUP = MITER + GO TO 220 + 430 ICF = 2 + NCF = NCF + 1 + RMAX = 2.0E0 + TN = TOLD + I1 = NQNYH + 1 + DO 445 JB = 1,NQ + I1 = I1 - NYH +Cdir$ ivdep + DO 440 I = I1,NQNYH + 440 YH1(I) = YH1(I) - YH1(I+NYH) + 445 CONTINUE + IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 + IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670 + IF (NCF .EQ. MXNCF) GO TO 670 + RH = 0.25E0 + IPUP = MITER + IREDO = 1 + GO TO 170 +C----------------------------------------------------------------------- +C The corrector has converged. JCUR is set to 0 +C to signal that the Jacobian involved may need updating later. +C The local error test is made and control passes to statement 500 +C if it fails. +C----------------------------------------------------------------------- + 450 JCUR = 0 + IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) + IF (M .GT. 0) DSM = SVNORM (N, ACOR, EWT)/TESCO(2,NQ) + IF (DSM .GT. 1.0E0) GO TO 500 +C----------------------------------------------------------------------- +C After a successful step, update the YH array. +C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1. +C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for +C use in a possible order increase on the next step. +C If a change in H is considered, an increase or decrease in order +C by one is considered also. A change in H is made only if it is by a +C factor of at least 1.1. If not, IALTH is set to 3 to prevent +C testing for that many steps. +C----------------------------------------------------------------------- + KFLAG = 0 + IREDO = 0 + NST = NST + 1 + HU = H + NQU = NQ + DO 470 J = 1,L + DO 470 I = 1,N + 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) + IALTH = IALTH - 1 + IF (IALTH .EQ. 0) GO TO 520 + IF (IALTH .GT. 1) GO TO 700 + IF (L .EQ. LMAX) GO TO 700 + DO 490 I = 1,N + 490 YH(I,LMAX) = ACOR(I) + GO TO 700 +C----------------------------------------------------------------------- +C The error test failed. KFLAG keeps track of multiple failures. +C Restore TN and the YH array to their previous values, and prepare +C to try the step again. Compute the optimum step size for this or +C one lower order. After 2 or more failures, H is forced to decrease +C by a factor of 0.2 or less. +C----------------------------------------------------------------------- + 500 KFLAG = KFLAG - 1 + TN = TOLD + I1 = NQNYH + 1 + DO 515 JB = 1,NQ + I1 = I1 - NYH +Cdir$ ivdep + DO 510 I = I1,NQNYH + 510 YH1(I) = YH1(I) - YH1(I+NYH) + 515 CONTINUE + RMAX = 2.0E0 + IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660 + IF (KFLAG .LE. -3) GO TO 640 + IREDO = 2 + RHUP = 0.0E0 + GO TO 540 +C----------------------------------------------------------------------- +C Regardless of the success or failure of the step, factors +C RHDN, RHSM, and RHUP are computed, by which H could be multiplied +C at order NQ - 1, order NQ, or order NQ + 1, respectively. +C In the case of failure, RHUP = 0.0 to avoid an order increase. +C The largest of these is determined and the new order chosen +C accordingly. If the order is to be increased, we compute one +C additional scaled derivative. +C----------------------------------------------------------------------- + 520 RHUP = 0.0E0 + IF (L .EQ. LMAX) GO TO 540 + DO 530 I = 1,N + 530 SAVF(I) = ACOR(I) - YH(I,LMAX) + DUP = SVNORM (N, SAVF, EWT)/TESCO(3,NQ) + EXUP = 1.0E0/(L+1) + RHUP = 1.0E0/(1.4E0*DUP**EXUP + 0.0000014E0) + 540 EXSM = 1.0E0/L + RHSM = 1.0E0/(1.2E0*DSM**EXSM + 0.0000012E0) + RHDN = 0.0E0 + IF (NQ .EQ. 1) GO TO 560 + DDN = SVNORM (N, YH(1,L), EWT)/TESCO(1,NQ) + EXDN = 1.0E0/NQ + RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) + 560 IF (RHSM .GE. RHUP) GO TO 570 + IF (RHUP .GT. RHDN) GO TO 590 + GO TO 580 + 570 IF (RHSM .LT. RHDN) GO TO 580 + NEWQ = NQ + RH = RHSM + GO TO 620 + 580 NEWQ = NQ - 1 + RH = RHDN + IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0 + GO TO 620 + 590 NEWQ = L + RH = RHUP + IF (RH .LT. 1.1E0) GO TO 610 + R = EL(L)/L + DO 600 I = 1,N + 600 YH(I,NEWQ+1) = ACOR(I)*R + GO TO 630 + 610 IALTH = 3 + GO TO 700 + 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610 + IF (KFLAG .LE. -2) RH = MIN(RH,0.2E0) +C----------------------------------------------------------------------- +C If there is a change of order, reset NQ, l, and the coefficients. +C In any case H is reset according to RH and the YH array is rescaled. +C Then exit from 690 if the step was OK, or redo the step otherwise. +C----------------------------------------------------------------------- + IF (NEWQ .EQ. NQ) GO TO 170 + 630 NQ = NEWQ + L = NQ + 1 + IRET = 2 + GO TO 150 +C----------------------------------------------------------------------- +C Control reaches this section if 3 or more failures have occured. +C If 10 failures have occurred, exit with KFLAG = -1. +C It is assumed that the derivatives that have accumulated in the +C YH array have errors of the wrong order. Hence the first +C derivative is recomputed, and the order is set to 1. Then +C H is reduced by a factor of 10, and the step is retried, +C until it succeeds or H reaches HMIN. +C----------------------------------------------------------------------- + 640 IF (KFLAG .EQ. -10) GO TO 660 + RH = 0.1E0 + RH = MAX(HMIN/ABS(H),RH) + H = H*RH + DO 645 I = 1,N + 645 Y(I) = YH(I,1) + CALL F (NEQ, TN, Y, SAVF) + NFE = NFE + 1 + DO 650 I = 1,N + 650 YH(I,2) = H*SAVF(I) + IPUP = MITER + IALTH = 5 + IF (NQ .EQ. 1) GO TO 200 + NQ = 1 + L = 2 + IRET = 3 + GO TO 150 +C----------------------------------------------------------------------- +C All returns are made through this section. H is saved in HOLD +C to allow the caller to change H on the next step. +C----------------------------------------------------------------------- + 660 KFLAG = -1 + GO TO 720 + 670 KFLAG = -2 + GO TO 720 + 680 KFLAG = -3 + GO TO 720 + 690 RMAX = 10.0E0 + 700 R = 1.0E0/TESCO(2,NQU) + DO 710 I = 1,N + 710 ACOR(I) = ACOR(I)*R + 720 HOLD = H + JSTART = 1 + RETURN +C----------------------- END OF SUBROUTINE SSTODE ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odepack/svnorm.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,34 @@ + REAL FUNCTION SVNORM (N, V, W) +C***BEGIN PROLOGUE SVNORM +C***SUBSIDIARY +C***PURPOSE Weighted root-mean-square vector norm. +C***TYPE SINGLE PRECISION (SVNORM-S, DVNORM-D) +C***AUTHOR Hindmarsh, Alan C., (LLNL) +C***DESCRIPTION +C +C This function routine computes the weighted root-mean-square norm +C of the vector of length N contained in the array V, with weights +C contained in the array W of length N: +C SVNORM = SQRT( (1/N) * SUM( V(i)*W(i) )**2 ) +C +C***SEE ALSO SLSODE +C***ROUTINES CALLED (NONE) +C***REVISION HISTORY (YYMMDD) +C 791129 DATE WRITTEN +C 890501 Modified prologue to SLATEC/LDOC format. (FNF) +C 890503 Minor cosmetic changes. (FNF) +C 930809 Renamed to allow single/double precision versions. (ACH) +C***END PROLOGUE SVNORM +C**End + INTEGER N, I + REAL V, W, SUM + DIMENSION V(N), W(N) +C +C***FIRST EXECUTABLE STATEMENT SVNORM + SUM = 0.0E0 + DO 10 I = 1,N + 10 SUM = SUM + (V(I)*W(I))**2 + SVNORM = SQRT(SUM/N) + RETURN +C----------------------- END OF FUNCTION SVNORM ------------------------ + END
--- a/libcruft/ordered-qz/Makefile.in Wed May 21 09:36:46 2008 -0400 +++ b/libcruft/ordered-qz/Makefile.in Sun May 11 22:51:50 2008 +0200 @@ -26,7 +26,7 @@ EXTERNAL_DISTFILES = $(DISTFILES) -FSRC = dsubsp.f exchqz.f +FSRC = dsubsp.f exchqz.f ssubsp.f sexchqz.f include $(TOPDIR)/Makeconf
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/ordered-qz/sexchqz.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,261 @@ + SUBROUTINE SEXCHQZ(NMAX, N, A, B, Z, L, LS1, LS2, EPS, FAIL) + INTEGER NMAX, N, L, LS1, LS2 + REAL A(NMAX,N), B(NMAX,N), Z(NMAX,N), EPS + LOGICAL FAIL +c modified july 9, 1998 a.s.hodel@eng.auburn.edu: +c calls to AMAX1 changed to call MAX instead. +c calls to giv changed to slartg (LAPACK); required new variable tempr +C* +C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A +C* WITH CONSECUTIVE LS1XLS1 AND LS2XLS2 DIAGONAL BLOCKS (LS1,LS2.LE.2) +C* STARTING AT ROW/COLUMN L, EXCHQZ PRODUCES EQUIVALENCE TRANSFORMA- +C* TIONS QT AND ZT THAT EXCHANGE THE BLOCKS ALONG WITH THEIR GENERALIZED +C* EIGENVALUES. EXCHQZ REQUIRES THE SUBROUTINES SROT (BLAS) AND GIV. +C* THE PARAMETERS IN THE CALLING SEQUENCE ARE (STARRED PARAMETERS ARE +C* ALTERED BY THE SUBROUTINE): +C* +C* NMAX THE FIRST DIMENSION OF A, B AND Z +C* N THE ORDER OF A, B AND Z +C* *A,*B THE MATRIX PAIR WHOSE BLOCKS ARE TO BE INTERCHANGED +C* *Z UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN +C* TRANSFORMATION ZT. +C* L THE POSITION OF THE BLOCKS +C* LS1 THE SIZE OF THE FIRST BLOCK +C* LS2 THE SIZE OF THE SECOND BLOCK +C* EPS THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT +C* *FAIL A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN, +C* TRUE OTHERWISE. +C* + INTEGER I, J, L1, L2, L3, LI, LJ, LL, IT1, IT2 + REAL U(3,3), D, E, F, G, SA, SB, A11B11, A21B11, + * A12B22, B12B22, + * A22B22, AMMBMM, ANMBMM, AMNBNN, BMNBNN, ANNBNN, TEMPR + LOGICAL ALTB + FAIL = .FALSE. + L1 = L + 1 + LL = LS1 + LS2 + IF (LL.GT.2) GO TO 10 +C*** INTERCHANGE 1X1 AND 1X1 BLOCKS VIA AN EQUIVALENCE +C*** TRANSFORMATION A:=Q*A*Z , B:=Q*B*Z +C*** WHERE Q AND Z ARE GIVENS ROTATIONS + F = MAX(ABS(A(L1,L1)),ABS(B(L1,L1))) + ALTB = .TRUE. + IF (ABS(A(L1,L1)).GE.F) ALTB = .FALSE. + SA = A(L1,L1)/F + SB = B(L1,L1)/F + F = SA*B(L,L) - SB*A(L,L) +C* CONSTRUCT THE COLUMN TRANSFORMATION Z + G = SA*B(L,L1) - SB*A(L,L1) + CALL SLARTG(F, G, D, E,TEMPR) + CALL SROT(L1, A(1,L), 1, A(1,L1), 1, E, -D) + CALL SROT(L1, B(1,L), 1, B(1,L1), 1, E, -D) + CALL SROT(N, Z(1,L), 1, Z(1,L1), 1, E, -D) +C* CONSTRUCT THE ROW TRANSFORMATION Q + IF (ALTB) CALL SLARTG(B(L,L), B(L1,L), D, E,TEMPR) + IF (.NOT.ALTB) CALL SLARTG(A(L,L), A(L1,L), D, E,TEMPR) + CALL SROT(N-L+1, A(L,L), NMAX, A(L1,L), NMAX, D, E) + CALL SROT(N-L+1, B(L,L), NMAX, B(L1,L), NMAX, D, E) + A(L1,L) = 0. + B(L1,L) = 0. + RETURN +C*** INTERCHANGE 1X1 AND 2X2 BLOCKS VIA AN EQUIVALENCE +C*** TRANSFORMATION A:=Q2*Q1*A*Z1*Z2 , B:=Q2*Q1*B*Z1*Z2 +C*** WHERE EACH QI AND ZI IS A GIVENS ROTATION + 10 L2 = L + 2 + IF (LS1.EQ.2) GO TO 60 + G = MAX(ABS(A(L,L)),ABS(B(L,L))) + ALTB = .TRUE. + IF (ABS(A(L,L)).LT.G) GO TO 20 + ALTB = .FALSE. + CALL SLARTG(A(L1,L1), A(L2,L1), D, E,TEMPR) + CALL SROT(N-L, A(L1,L1), NMAX, A(L2,L1), NMAX, D, E) + CALL SROT(N-L, B(L1,L1), NMAX, B(L2,L1), NMAX, D, E) +C** EVALUATE THE PENCIL AT THE EIGENVALUE CORRESPONDING +C** TO THE 1X1 BLOCK + 20 SA = A(L,L)/G + SB = B(L,L)/G + DO 40 J=1,2 + LJ = L + J + DO 30 I=1,3 + LI = L + I - 1 + U(I,J) = SA*B(LI,LJ) - SB*A(LI,LJ) + 30 CONTINUE + 40 CONTINUE + CALL SLARTG(U(3,1), U(3,2), D, E,TEMPR) + CALL SROT(3, U(1,1), 1, U(1,2), 1, E, -D) +C* PERFORM THE ROW TRANSFORMATION Q1 + CALL SLARTG(U(1,1), U(2,1), D, E,TEMPR) + U(2,2) = -U(1,2)*E + U(2,2)*D + CALL SROT(N-L+1, A(L,L), NMAX, A(L1,L), NMAX, D, E) + CALL SROT(N-L+1, B(L,L), NMAX, B(L1,L), NMAX, D, E) +C* PERFORM THE COLUMN TRANSFORMATION Z1 + IF (ALTB) CALL SLARTG(B(L1,L), B(L1,L1), D, E,TEMPR) + IF (.NOT.ALTB) CALL SLARTG(A(L1,L), A(L1,L1), D, E,TEMPR) + CALL SROT(L2, A(1,L), 1, A(1,L1), 1, E, -D) + CALL SROT(L2, B(1,L), 1, B(1,L1), 1, E, -D) + CALL SROT(N, Z(1,L), 1, Z(1,L1), 1, E, -D) +C* PERFORM THE ROW TRANSFORMATION Q2 + CALL SLARTG(U(2,2), U(3,2), D, E,TEMPR) + CALL SROT(N-L+1, A(L1,L), NMAX, A(L2,L), NMAX, D, E) + CALL SROT(N-L+1, B(L1,L), NMAX, B(L2,L), NMAX, D, E) +C* PERFORM THE COLUMN TRANSFORMATION Z2 + IF (ALTB) CALL SLARTG(B(L2,L1), B(L2,L2), D, E,TEMPR) + IF (.NOT.ALTB) CALL SLARTG(A(L2,L1), A(L2,L2), D, E,TEMPR) + CALL SROT(L2, A(1,L1), 1, A(1,L2), 1, E, -D) + CALL SROT(L2, B(1,L1), 1, B(1,L2), 1, E, -D) + CALL SROT(N, Z(1,L1), 1, Z(1,L2), 1, E, -D) + IF (ALTB) GO TO 50 + CALL SLARTG(B(L,L), B(L1,L), D, E,TEMPR) + CALL SROT(N-L+1, A(L,L), NMAX, A(L1,L), NMAX, D, E) + CALL SROT(N-L+1, B(L,L), NMAX, B(L1,L), NMAX, D, E) +C* PUT THE NEGLECTABLE ELEMENTS EQUAL TO ZERO + 50 A(L2,L) = 0. + A(L2,L1) = 0. + B(L1,L) = 0. + B(L2,L) = 0. + B(L2,L1) = 0. + RETURN +C*** INTERCHANGE 2X2 AND 1X1 BLOCKS VIA AN EQUIVALENCE +C*** TRANSFORMATION A:=Q2*Q1*A*Z1*Z2 , B:=Q2*Q1*B*Z1*Z2 +C*** WHERE EACH QI AND ZI IS A GIVENS ROTATION + 60 IF (LS2.EQ.2) GO TO 110 + G = MAX(ABS(A(L2,L2)),ABS(B(L2,L2))) + ALTB = .TRUE. + IF (ABS(A(L2,L2)).LT.G) GO TO 70 + ALTB = .FALSE. + CALL SLARTG(A(L,L), A(L1,L), D, E,TEMPR) + CALL SROT(N-L+1, A(L,L), NMAX, A(L1,L), NMAX, D, E) + CALL SROT(N-L+1, B(L,L), NMAX, B(L1,L), NMAX, D, E) +C** EVALUATE THE PENCIL AT THE EIGENVALUE CORRESPONDING +C** TO THE 1X1 BLOCK + 70 SA = A(L2,L2)/G + SB = B(L2,L2)/G + DO 90 I=1,2 + LI = L + I - 1 + DO 80 J=1,3 + LJ = L + J - 1 + U(I,J) = SA*B(LI,LJ) - SB*A(LI,LJ) + 80 CONTINUE + 90 CONTINUE + CALL SLARTG(U(1,1), U(2,1), D, E,TEMPR) + CALL SROT(3, U(1,1), 3, U(2,1), 3, D, E) +C* PERFORM THE COLUMN TRANSFORMATION Z1 + CALL SLARTG(U(2,2), U(2,3), D, E,TEMPR) + U(1,2) = U(1,2)*E - U(1,3)*D + CALL SROT(L2, A(1,L1), 1, A(1,L2), 1, E, -D) + CALL SROT(L2, B(1,L1), 1, B(1,L2), 1, E, -D) + CALL SROT(N, Z(1,L1), 1, Z(1,L2), 1, E, -D) +C* PERFORM THE ROW TRANSFORMATION Q1 + IF (ALTB) CALL SLARTG(B(L1,L1), B(L2,L1), D, E,TEMPR) + IF (.NOT.ALTB) CALL SLARTG(A(L1,L1), A(L2,L1), D, E,TEMPR) + CALL SROT(N-L+1, A(L1,L), NMAX, A(L2,L), NMAX, D, E) + CALL SROT(N-L+1, B(L1,L), NMAX, B(L2,L), NMAX, D, E) +C* PERFORM THE COLUMN TRANSFORMATION Z2 + CALL SLARTG(U(1,1), U(1,2), D, E,TEMPR) + CALL SROT(L2, A(1,L), 1, A(1,L1), 1, E, -D) + CALL SROT(L2, B(1,L), 1, B(1,L1), 1, E, -D) + CALL SROT(N, Z(1,L), 1, Z(1,L1), 1, E, -D) +C* PERFORM THE ROW TRANSFORMATION Q2 + IF (ALTB) CALL SLARTG(B(L,L), B(L1,L), D, E,TEMPR) + IF (.NOT.ALTB) CALL SLARTG(A(L,L), A(L1,L), D, E,TEMPR) + CALL SROT(N-L+1, A(L,L), NMAX, A(L1,L), NMAX, D, E) + CALL SROT(N-L+1, B(L,L), NMAX, B(L1,L), NMAX, D, E) + IF (ALTB) GO TO 100 + CALL SLARTG(B(L1,L1), B(L2,L1), D, E,TEMPR) + CALL SROT(N-L, A(L1,L1), NMAX, A(L2,L1), NMAX, D, E) + CALL SROT(N-L, B(L1,L1), NMAX, B(L2,L1), NMAX, D, E) +C* PUT THE NEGLECTABLE ELEMENTS EQUAL TO ZERO + 100 A(L1,L) = 0. + A(L2,L) = 0. + B(L1,L) = 0. + B(L2,1) = 0. + B(L2,L1) = 0. + RETURN +C*** INTERCHANGE 2X2 AND 2X2 BLOCKS VIA A SEQUENCE OF +C*** QZ-STEPS REALIZED BY THE EQUIVALENCE TRANSFORMATIONS +C*** A:=Q5*Q4*Q3*Q2*Q1*A*Z1*Z2*Z3*Z4*Z5 +C*** B:=Q5*Q4*Q3*Q2*Q1*B*Z1*Z2*Z3*Z4*Z5 +C*** WHERE EACH QI AND ZI IS A GIVENS ROTATION + 110 L3 = L + 3 +C* COMPUTE IMPLICIT SHIFT + AMMBMM = A(L,L)/B(L,L) + ANMBMM = A(L1,L)/B(L,L) + AMNBNN = A(L,L1)/B(L1,L1) + ANNBNN = A(L1,L1)/B(L1,L1) + BMNBNN = B(L,L1)/B(L1,L1) + DO 130 IT1=1,3 + U(1,1) = 1. + U(2,1) = 1. + U(3,1) = 1. + DO 120 IT2=1,10 +C* PERFORM ROW TRANSFORMATIONS Q1 AND Q2 + CALL SLARTG(U(2,1), U(3,1), D, E,TEMPR) + CALL SROT(N-L+1, A(L1,L), NMAX, A(L2,L), NMAX, D, E) + CALL SROT(N-L, B(L1,L1), NMAX, B(L2,L1), NMAX, D, E) + U(2,1) = D*U(2,1) + E*U(3,1) + CALL SLARTG(U(1,1), U(2,1), D, E,TEMPR) + CALL SROT(N-L+1, A(L,L), NMAX, A(L1,L), NMAX, D, E) + CALL SROT(N-L+1, B(L,L), NMAX, B(L1,L), NMAX, D, E) +C* PERFORM COLUMN TRANSFORMATIONS Z1 AND Z2 + CALL SLARTG(B(L2,L1), B(L2,L2), D, E,TEMPR) + CALL SROT(L3, A(1,L1), 1, A(1,L2), 1, E, -D) + CALL SROT(L2, B(1,L1), 1, B(1,L2), 1, E, -D) + CALL SROT(N, Z(1,L1), 1, Z(1,L2), 1, E, -D) + CALL SLARTG(B(L1,L), B(L1,L1), D, E,TEMPR) + CALL SROT(L3, A(1,L), 1, A(1,L1), 1, E, -D) + CALL SROT(L1, B(1,L), 1, B(1,L1), 1, E, -D) + CALL SROT(N, Z(1,L), 1, Z(1,L1), 1, E, -D) +C* PERFORM TRANSFORMATIONS Q3,Z3,Q4,Z4,Q5 AND Z5 IN +C* ORDER TO REDUCE THE PENCIL TO HESSENBERG FORM + CALL SLARTG(A(L2,L), A(L3,L), D, E,TEMPR) + CALL SROT(N-L+1, A(L2,L), NMAX, A(L3,L), NMAX, D, E) + CALL SROT(N-L1, B(L2,L2), NMAX, B(L3,L2), NMAX, D, E) + CALL SLARTG(B(L3,L2), B(L3,L3), D, E,TEMPR) + CALL SROT(L3, A(1,L2), 1, A(1,L3), 1, E, -D) + CALL SROT(L3, B(1,L2), 1, B(1,L3), 1, E, -D) + CALL SROT(N, Z(1,L2), 1, Z(1,L3), 1, E, -D) + CALL SLARTG(A(L1,L), A(L2,L), D, E,TEMPR) + CALL SROT(N-L+1, A(L1,L), NMAX, A(L2,L), NMAX, D, E) + CALL SROT(N-L, B(L1,L1), NMAX, B(L2,L1), NMAX, D, E) + CALL SLARTG(B(L2,L1), B(L2,L2), D, E,TEMPR) + CALL SROT(L3, A(1,L1), 1, A(1,L2), 1, E, -D) + CALL SROT(L2, B(1,L1), 1, B(1,L2), 1, E, -D) + CALL SROT(N, Z(1,L1), 1, Z(1,L2), 1, E, -D) + CALL SLARTG(A(L2,L1), A(L3,L1), D, E,TEMPR) + CALL SROT(N-L, A(L2,L1), NMAX, A(L3,L1), NMAX, D, E) + CALL SROT(N-L1, B(L2,L2), NMAX, B(L3,L2), NMAX, D, E) + CALL SLARTG(B(L3,L2), B(L3,L3), D, E,TEMPR) + CALL SROT(L3, A(1,L2), 1, A(1,L3), 1, E, -D) + CALL SROT(L3, B(1,L2), 1, B(1,L3), 1, E, -D) + CALL SROT(N, Z(1,L2), 1, Z(1,L3), 1, E, -D) +C* TEST OF CONVERGENCE ON THE ELEMENT SEPARATING THE BLOCKS + IF (ABS(A(L2,L1)).LE.EPS) GO TO 140 +C* COMPUTE A NEW SHIFT IN CASE OF NO CONVERGENCE + A11B11 = A(L,L)/B(L,L) + A12B22 = A(L,L1)/B(L1,L1) + A21B11 = A(L1,L)/B(L,L) + A22B22 = A(L1,L1)/B(L1,L1) + B12B22 = B(L,L1)/B(L1,L1) + U(1,1) = ((AMMBMM-A11B11)*(ANNBNN-A11B11)-AMNBNN* + * ANMBMM+ANMBMM*BMNBNN*A11B11)/A21B11 + A12B22 - A11B11*B12B22 + U(2,1) = (A22B22-A11B11) - A21B11*B12B22 - (AMMBMM-A11B11) - + * (ANNBNN-A11B11) + ANMBMM*BMNBNN + U(3,1) = A(L2,L1)/B(L1,L1) + 120 CONTINUE + 130 CONTINUE + FAIL = .TRUE. + RETURN +C* PUT THE NEGLECTABLE ELEMENTS EQUAL TO ZERO IN +C* CASE OF CONVERGENCE + 140 A(L2,L) = 0. + A(L2,L1) = 0. + A(L3,L) = 0. + A(L3,L1) = 0. + B(L1,L) = 0. + B(L2,L) = 0. + B(L2,L1) = 0. + B(L3,L) = 0. + B(L3,L1) = 0. + B(L3,L2) = 0. + RETURN + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/ordered-qz/ssubsp.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,104 @@ + SUBROUTINE SSUBSP(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND) + INTEGER NMAX, N, FTEST, NDIM, IND(N) + LOGICAL FAIL + REAL A(NMAX,N), B(NMAX,N), Z(NMAX,N), EPS +C* +C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A +C* WITH 1X1 OR 2X2 DIAGONAL BLOCKS, THIS ROUTINE REORDERS THE DIAGONAL +C* BLOCKS ALONG WITH THEIR GENERALIZED EIGENVALUES BY CONSTRUCTING EQUI- +C* VALENCE TRANSFORMATIONS QT AND ZT. THE ROW TRANSFORMATION ZT IS ALSO +C* PERFORMED ON THE GIVEN (INITIAL) TRANSFORMATION Z (RESULTING FROM A +C* POSSIBLE PREVIOUS STEP OR INITIALIZED WITH THE IDENTITY MATRIX). +C* AFTER REORDERING, THE EIGENVALUES INSIDE THE REGION SPECIFIED BY THE +C* FUNCTION FTEST APPEAR AT THE TOP. IF NDIM IS THEIR NUMBER THEN THE +C* NDIM FIRST COLUMNS OF Z SPAN THE REQUESTED SUBSPACE. DSUBSP REQUIRES +C* THE SUBROUTINE EXCHQZ AND THE INTEGER FUNCTION FTEST WHICH HAS TO BE +C* PROVIDED BY THE USER. THE PARAMETERS IN THE CALLING SEQUENCE ARE : +C* (STARRED PARAMETERS ARE ALTERED BY THE SUBROUTINE) +C* +C* NMAX THE FIRST DIMENSION OF A, B AND Z +C* N THE ORDER OF A, B AND Z +C* *A,*B THE MATRIX PAIR WHOSE BLOCKS ARE TO BE REORDERED. +C* *Z UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN +C* TRANSFORMATION ZT. +C* FTEST(LS,ALPHA,BETA,S,P) AN INTEGER FUNCTION DESCRIBING THE +C* SPECTRUM OF THE DEFLATING SUBSPACE TO BE COMPUTED: +C* WHEN LS=1 FTEST CHECKS IF ALPHA/BETA IS IN THAT SPECTRUM +C* WHEN LS=2 FTEST CHECKS IF THE TWO COMPLEX CONJUGATE +C* ROOTS WITH SUM S AND PRODUCT P ARE IN THAT SPECTRUM +C* IF THE ANSWER IS POSITIVE, FTEST=1, OTHERWISE FTEST=-1 +C* EPS THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT +C* *NDIM AN INTEGER GIVING THE DIMENSION OF THE COMPUTED +C* DEFLATING SUBSPACE +C* *FAIL A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN, +C* TRUE OTHERWISE (WHEN SEXCHQZ FAILS) +C* *IND AN INTEGER WORKING ARRAY OF DIMENSION AT LEAST N +C* + INTEGER L, LS, LS1, LS2, L1, LL, NUM, IS, L2I, L2K, I, K, II, + * ISTEP, IFIRST + REAL S, P, D, ALPHA, BETA + FAIL = .TRUE. + NDIM = 0 + NUM = 0 + L = 0 + LS = 1 +C*** CONSTRUCT ARRAY IND(I) WHERE : +C*** IABS(IND(I)) IS THE SIZE OF THE BLOCK I +C*** SIGN(IND(I)) INDICATES THE LOCATION OF ITS EIGENVALUES +C*** (AS DETERMINED BY FTEST). +C*** NUM IS THE NUMBER OF ELEMENTS IN THIS ARRAY + DO 30 LL=1,N + L = L + LS + IF (L.GT.N) GO TO 40 + L1 = L + 1 + IF (L1.GT.N) GO TO 10 + IF (A(L1,L).EQ.0.) GO TO 10 +C* HERE A 2X2 BLOCK IS CHECKED * + LS = 2 + D = B(L,L)*B(L1,L1) + S = (A(L,L)*B(L1,L1)+A(L1,L1)*B(L,L)-A(L1,L)*B(L,L1))/D + P = (A(L,L)*A(L1,L1)-A(L,L1)*A(L1,L))/D + IS = FTEST(LS,ALPHA,BETA,S,P) + GO TO 20 +C* HERE A 1X1 BLOCK IS CHECKED * + 10 LS = 1 + IS = FTEST(LS,A(L,L),B(L,L),S,P) + 20 NUM = NUM + 1 + IF (IS.EQ.1) NDIM = NDIM + LS + IND(NUM) = LS*IS + 30 CONTINUE +C*** REORDER BLOCKS SUCH THAT THOSE WITH POSITIVE VALUE +C*** OF IND(.) APPEAR FIRST. + 40 L2I = 1 + DO 100 I=1,NUM + IF (IND(I).GT.0) GO TO 90 +C* IF A NEGATIVE IND(I) IS ENCOUNTERED, THEN SEARCH FOR THE FIRST +C* POSITIVE IND(K) FOLLOWING ON IT + L2K = L2I + DO 60 K=I,NUM + IF (IND(K).LT.0) GO TO 50 + GO TO 70 + 50 L2K = L2K - IND(K) + 60 CONTINUE +C* IF THERE ARE NO POSITIVE INDICES FOLLOWING ON A NEGATIVE ONE +C* THEN STOP + GO TO 110 +C* IF A POSITIVE IND(K) FOLLOWS ON A NEGATIVE IND(I) THEN +C* INTERCHANGE BLOCK K BEFORE BLOCK I BY PERFORMING K-I SWAPS + 70 ISTEP = K - I + LS2 = IND(K) + L = L2K + DO 80 II=1,ISTEP + IFIRST = K - II + LS1 = -IND(IFIRST) + L = L - LS1 + CALL SEXCHQZ(NMAX, N, A, B, Z, L, LS1, LS2, EPS, FAIL) + IF (FAIL) RETURN + IND(IFIRST+1) = IND(IFIRST) + 80 CONTINUE + IND(I) = LS2 + 90 L2I = L2I + IND(I) + 100 CONTINUE + 110 FAIL = .FALSE. + RETURN + END
--- a/libcruft/quadpack/Makefile.in Wed May 21 09:36:46 2008 -0400 +++ b/libcruft/quadpack/Makefile.in Sun May 11 22:51:50 2008 +0200 @@ -27,7 +27,8 @@ EXTERNAL_DISTFILES = $(DISTFILES) FSRC = dqagi.f dqagie.f dqagp.f dqagpe.f dqelg.f dqk15i.f \ - dqk21.f dqpsrt.f xerror.f + dqk21.f dqpsrt.f qagie.f qagi.f qagpe.f qagp.f qelg.f \ + qk15i.f qk21.f qpsrt.f xerror.f include $(TOPDIR)/Makeconf
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qagi.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,190 @@ + subroutine qagi(f,bound,inf,epsabs,epsrel,result,abserr,neval, + * ier,limit,lenw,last,iwork,work) +c***begin prologue qagi +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a3a1,h2a4a1 +c***keywords automatic integrator, infinite intervals, +c general-purpose, transformation, extrapolation, +c globally adaptive +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. -k.u.leuven +c***purpose the routine calculates an approximation result to a given +c integral i = integral of f over (bound,+infinity) +c or i = integral of f over (-infinity,bound) +c or i = integral of f over (-infinity,+infinity) +c hopefully satisfying following claim for accuracy +c abs(i-result).le.max(epsabs,epsrel*abs(i)). +c***description +c +c integration over infinite intervals +c standard fortran subroutine +c +c parameters +c on entry +c f - subroutine f(x,result) defining the integrand +c function f(x). the actual name for f needs to be +c declared e x t e r n a l in the driver program. +c +c bound - real +c finite bound of integration range +c (has no meaning if interval is doubly-infinite) +c +c inf - integer +c indicating the kind of integration range involved +c inf = 1 corresponds to (bound,+infinity), +c inf = -1 to (-infinity,bound), +c inf = 2 to (-infinity,+infinity). +c +c epsabs - real +c absolute accuracy requested +c epsrel - real +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c +c on return +c result - real +c approximation to the integral +c +c abserr - real +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c - ier.gt.0 abnormal termination of the routine. the +c estimates for result and error are less +c reliable. it is assumed that the requested +c accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more +c subdivisions by increasing the value of +c limit (and taking the according dimension +c adjustments into account). however, if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. if +c the position of a local difficulty can be +c determined (e.g. singularity, +c discontinuity within the interval) one +c will probably gain from splitting up the +c interval at this point and calling the +c integrator on the subranges. if possible, +c an appropriate special-purpose integrator +c should be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is +c detected, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour occurs +c at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. +c it is assumed that the requested tolerance +c cannot be achieved, and that the returned +c result is the best which can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier. +c = 6 the input is invalid, because +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) +c or limit.lt.1 or leniw.lt.limit*4. +c result, abserr, neval, last are set to +c zero. exept when limit or leniw is +c invalid, iwork(1), work(limit*2+1) and +c work(limit*3+1) are set to zero, work(1) +c is set to a and work(limit+1) to b. +c +c dimensioning parameters +c limit - integer +c dimensioning parameter for iwork +c limit determines the maximum number of subintervals +c in the partition of the given integration interval +c (a,b), limit.ge.1. +c if limit.lt.1, the routine will end with ier = 6. +c +c lenw - integer +c dimensioning parameter for work +c lenw must be at least limit*4. +c if lenw.lt.limit*4, the routine will end +c with ier = 6. +c +c last - integer +c on return, last equals the number of subintervals +c produced in the subdivision process, which +c determines the number of significant elements +c actually in the work arrays. +c +c work arrays +c iwork - integer +c vector of dimension at least limit, the first +c k elements of which contain pointers +c to the error estimates over the subintervals, +c such that work(limit*3+iwork(1)),... , +c work(limit*3+iwork(k)) form a decreasing +c sequence, with k = last if last.le.(limit/2+2), and +c k = limit+1-last otherwise +c +c work - real +c vector of dimension at least lenw +c on return +c work(1), ..., work(last) contain the left +c end points of the subintervals in the +c partition of (a,b), +c work(limit+1), ..., work(limit+last) contain +c the right end points, +c work(limit*2+1), ...,work(limit*2+last) contain the +c integral approximations over the subintervals, +c work(limit*3+1), ..., work(limit*3) +c contain the error estimates. +c***references (none) +c***routines called qagie,xerror +c***end prologue qagi +c + real abserr, epsabs,epsrel,result,work + integer ier,iwork, lenw,limit,lvl,l1,l2,l3,neval +c + dimension iwork(limit),work(lenw) +c + external f +c +c check validity of limit and lenw. +c +c***first executable statement qagi + ier = 6 + neval = 0 + last = 0 + result = 0.0e+00 + abserr = 0.0e+00 + if(limit.lt.1.or.lenw.lt.limit*4) go to 10 +c +c prepare call for qagie. +c + l1 = limit+1 + l2 = limit+l1 + l3 = limit+l2 +c + call qagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, + * neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) +c +c call error handler if necessary. +c + lvl = 0 +10 if(ier.eq.6) lvl = 1 + if(ier.ne.0) call xerror('abnormal return from qagi',26,ier,lvl) + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qagie.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,460 @@ + subroutine qagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, + * neval,ier,alist,blist,rlist,elist,iord,last) +c***begin prologue qagie +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a3a1,h2a4a1 +c***keywords automatic integrator, infinite intervals, +c general-purpose, transformation, extrapolation, +c globally adaptive +c***author piessens,robert,appl. math & progr. div - k.u.leuven +c de doncker,elise,appl. math & progr. div - k.u.leuven +c***purpose the routine calculates an approximation result to a given +c integral i = integral of f over (bound,+infinity) +c or i = integral of f over (-infinity,bound) +c or i = integral of f over (-infinity,+infinity), +c hopefully satisfying following claim for accuracy +c abs(i-result).le.max(epsabs,epsrel*abs(i)) +c***description +c +c integration over infinite intervals +c standard fortran subroutine +c +c f - subroutine f(x,ierr,result) defining the integrand +c function f(x). the actual name for f needs to be +c declared e x t e r n a l in the driver program. +c +c bound - real +c finite bound of integration range +c (has no meaning if interval is doubly-infinite) +c +c inf - real +c indicating the kind of integration range involved +c inf = 1 corresponds to (bound,+infinity), +c inf = -1 to (-infinity,bound), +c inf = 2 to (-infinity,+infinity). +c +c epsabs - real +c absolute accuracy requested +c epsrel - real +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c limit - integer +c gives an upper bound on the number of subintervals +c in the partition of (a,b), limit.ge.1 +c +c on return +c result - real +c approximation to the integral +c +c abserr - real +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c - ier.gt.0 abnormal termination of the routine. the +c estimates for result and error are less +c reliable. it is assumed that the requested +c accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more +c subdivisions by increasing the value of +c limit (and taking the according dimension +c adjustments into account). however,if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. +c if the position of a local difficulty can +c be determined (e.g. singularity, +c discontinuity within the interval) one +c will probably gain from splitting up the +c interval at this point and calling the +c integrator on the subranges. if possible, +c an appropriate special-purpose integrator +c should be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is +c detected, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour occurs +c at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. +c it is assumed that the requested tolerance +c cannot be achieved, and that the returned +c result is the best which can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier. +c = 6 the input is invalid, because +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c result, abserr, neval, last, rlist(1), +c elist(1) and iord(1) are set to zero. +c alist(1) and blist(1) are set to 0 +c and 1 respectively. +c +c alist - real +c vector of dimension at least limit, the first +c last elements of which are the left +c end points of the subintervals in the partition +c of the transformed integration range (0,1). +c +c blist - real +c vector of dimension at least limit, the first +c last elements of which are the right +c end points of the subintervals in the partition +c of the transformed integration range (0,1). +c +c rlist - real +c vector of dimension at least limit, the first +c last elements of which are the integral +c approximations on the subintervals +c +c elist - real +c vector of dimension at least limit, the first +c last elements of which are the moduli of the +c absolute error estimates on the subintervals +c +c iord - integer +c vector of dimension limit, the first k +c elements of which are pointers to the +c error estimates over the subintervals, +c such that elist(iord(1)), ..., elist(iord(k)) +c form a decreasing sequence, with k = last +c if last.le.(limit/2+2), and k = limit+1-last +c otherwise +c +c last - integer +c number of subintervals actually produced +c in the subdivision process +c +c***references (none) +c***routines called qelg,qk15i,qpsrt,r1mach +c***end prologue qagie +c + real abseps,abserr,alist,area,area1,area12,area2,a1, + * a2,blist,boun,bound,b1,b2,correc,defabs,defab1,defab2, + * dres,r1mach,elist,epmach,epsabs,epsrel,erlarg,erlast, + * errbnd,errmax,error1,error2,erro12,errsum,ertest,oflow,resabs, + * reseps,result,res3la,rlist,rlist2,small,uflow + integer id,ier,ierro,inf,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn, + * ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2 + logical extrap,noext +c + dimension alist(limit),blist(limit),elist(limit),iord(limit), + * res3la(3),rlist(limit),rlist2(52) +c + external f +c +c the dimension of rlist2 is determined by the value of +c limexp in subroutine qelg. +c +c +c list of major variables +c ----------------------- +c +c alist - list of left end points of all subintervals +c considered up to now +c blist - list of right end points of all subintervals +c considered up to now +c rlist(i) - approximation to the integral over +c (alist(i),blist(i)) +c rlist2 - array of dimension at least (limexp+2), +c containing the part of the epsilon table +c wich is still needed for further computations +c elist(i) - error estimate applying to rlist(i) +c maxerr - pointer to the interval with largest error +c estimate +c errmax - elist(maxerr) +c erlast - error on the interval currently subdivided +c (before that subdivision has taken place) +c area - sum of the integrals over the subintervals +c errsum - sum of the errors over the subintervals +c errbnd - requested accuracy max(epsabs,epsrel* +c abs(result)) +c *****1 - variable for the left subinterval +c *****2 - variable for the right subinterval +c last - index for subdivision +c nres - number of calls to the extrapolation routine +c numrl2 - number of elements currently in rlist2. if an +c appropriate approximation to the compounded +c integral has been obtained, it is put in +c rlist2(numrl2) after numrl2 has been increased +c by one. +c small - length of the smallest interval considered up +c to now, multiplied by 1.5 +c erlarg - sum of the errors over the intervals larger +c than the smallest interval considered up to now +c extrap - logical variable denoting that the routine +c is attempting to perform extrapolation. i.e. +c before subdividing the smallest interval we +c try to decrease the value of erlarg. +c noext - logical variable denoting that extrapolation +c is no longer allowed (true-value) +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c oflow is the largest positive magnitude. +c + epmach = r1mach(4) +c +c test on validity of parameters +c ----------------------------- +c +c***first executable statement qagie + ier = 0 + neval = 0 + last = 0 + result = 0.0e+00 + abserr = 0.0e+00 + alist(1) = 0.0e+00 + blist(1) = 0.1e+01 + rlist(1) = 0.0e+00 + elist(1) = 0.0e+00 + iord(1) = 0 + if(epsabs.le.0.0e+00.and.epsrel.lt.amax1(0.5e+02*epmach,0.5e-14)) + * ier = 6 + if(ier.eq.6) go to 999 +c +c +c first approximation to the integral +c ----------------------------------- +c +c determine the interval to be mapped onto (0,1). +c if inf = 2 the integral is computed as i = i1+i2, where +c i1 = integral of f over (-infinity,0), +c i2 = integral of f over (0,+infinity). +c + boun = bound + if(inf.eq.2) boun = 0.0e+00 + call qk15i(f,boun,inf,0.0e+00,0.1e+01,result,abserr, + * defabs,resabs,ier) + if (ier.lt.0) return +c +c test on accuracy +c + last = 1 + rlist(1) = result + elist(1) = abserr + iord(1) = 1 + dres = abs(result) + errbnd = amax1(epsabs,epsrel*dres) + if(abserr.le.1.0e+02*epmach*defabs.and.abserr.gt. + * errbnd) ier = 2 + if(limit.eq.1) ier = 1 + if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or. + * abserr.eq.0.0e+00) go to 130 +c +c initialization +c -------------- +c + uflow = r1mach(1) + oflow = r1mach(2) + rlist2(1) = result + errmax = abserr + maxerr = 1 + area = result + errsum = abserr + abserr = oflow + nrmax = 1 + nres = 0 + ktmin = 0 + numrl2 = 2 + extrap = .false. + noext = .false. + ierro = 0 + iroff1 = 0 + iroff2 = 0 + iroff3 = 0 + ksgn = -1 + if(dres.ge.(0.1e+01-0.5e+02*epmach)*defabs) ksgn = 1 +c +c main do-loop +c ------------ +c + do 90 last = 2,limit +c +c bisect the subinterval with nrmax-th largest +c error estimate. +c + a1 = alist(maxerr) + b1 = 0.5e+00*(alist(maxerr)+blist(maxerr)) + a2 = b1 + b2 = blist(maxerr) + erlast = errmax + call qk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1,ier) + if (ier.lt.0) return + call qk15i(f,boun,inf,a2,b2,area2,error2,resabs,defab2,ier) + if (ier.lt.0) return +c +c improve previous approximations to integral +c and error and test for accuracy. +c + area12 = area1+area2 + erro12 = error1+error2 + errsum = errsum+erro12-errmax + area = area+area12-rlist(maxerr) + if(defab1.eq.error1.or.defab2.eq.error2)go to 15 + if(abs(rlist(maxerr)-area12).gt.0.1e-04*abs(area12) + * .or.erro12.lt.0.99e+00*errmax) go to 10 + if(extrap) iroff2 = iroff2+1 + if(.not.extrap) iroff1 = iroff1+1 + 10 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 + 15 rlist(maxerr) = area1 + rlist(last) = area2 + errbnd = amax1(epsabs,epsrel*abs(area)) +c +c test for roundoff error and eventually +c set error flag. +c + if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 + if(iroff2.ge.5) ierro = 3 +c +c set error flag in the case that the number of +c subintervals equals limit. +c + if(last.eq.limit) ier = 1 +c +c set error flag in the case of bad integrand behaviour +c at some points of the integration range. +c + if(amax1(abs(a1),abs(b2)).le.(0.1e+01+0.1e+03*epmach)* + * (abs(a2)+0.1e+04*uflow)) ier = 4 +c +c append the newly-created intervals to the list. +c + if(error2.gt.error1) go to 20 + alist(last) = a2 + blist(maxerr) = b1 + blist(last) = b2 + elist(maxerr) = error1 + elist(last) = error2 + go to 30 + 20 alist(maxerr) = a2 + alist(last) = a1 + blist(last) = b1 + rlist(maxerr) = area2 + rlist(last) = area1 + elist(maxerr) = error2 + elist(last) = error1 +c +c call subroutine qpsrt to maintain the descending ordering +c in the list of error estimates and select the +c subinterval with nrmax-th largest error estimate (to be +c bisected next). +c + 30 call qpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) + if(errsum.le.errbnd) go to 115 + if(ier.ne.0) go to 100 + if(last.eq.2) go to 80 + if(noext) go to 90 + erlarg = erlarg-erlast + if(abs(b1-a1).gt.small) erlarg = erlarg+erro12 + if(extrap) go to 40 +c +c test whether the interval to be bisected next is the +c smallest interval. +c + if(abs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 + extrap = .true. + nrmax = 2 + 40 if(ierro.eq.3.or.erlarg.le.ertest) go to 60 +c +c the smallest interval has the largest error. +c before bisecting decrease the sum of the errors +c over the larger intervals (erlarg) and perform +c extrapolation. +c + id = nrmax + jupbnd = last + if(last.gt.(2+limit/2)) jupbnd = limit+3-last + do 50 k = id,jupbnd + maxerr = iord(nrmax) + errmax = elist(maxerr) + if(abs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 + nrmax = nrmax+1 + 50 continue +c +c perform extrapolation. +c + 60 numrl2 = numrl2+1 + rlist2(numrl2) = area + call qelg(numrl2,rlist2,reseps,abseps,res3la,nres) + ktmin = ktmin+1 + if(ktmin.gt.5.and.abserr.lt.0.1e-02*errsum) ier = 5 + if(abseps.ge.abserr) go to 70 + ktmin = 0 + abserr = abseps + result = reseps + correc = erlarg + ertest = amax1(epsabs,epsrel*abs(reseps)) + if(abserr.le.ertest) go to 100 +c +c prepare bisection of the smallest interval. +c + 70 if(numrl2.eq.1) noext = .true. + if(ier.eq.5) go to 100 + maxerr = iord(1) + errmax = elist(maxerr) + nrmax = 1 + extrap = .false. + small = small*0.5e+00 + erlarg = errsum + go to 90 + 80 small = 0.375e+00 + erlarg = errsum + ertest = errbnd + rlist2(2) = area + 90 continue +c +c set final result and error estimate. +c ------------------------------------ +c + 100 if(abserr.eq.oflow) go to 115 + if((ier+ierro).eq.0) go to 110 + if(ierro.eq.3) abserr = abserr+correc + if(ier.eq.0) ier = 3 + if(result.ne.0.0e+00.and.area.ne.0.0e+00)go to 105 + if(abserr.gt.errsum)go to 115 + if(area.eq.0.0e+00) go to 130 + go to 110 + 105 if(abserr/abs(result).gt.errsum/abs(area))go to 115 +c +c test on divergence +c + 110 if(ksgn.eq.(-1).and.amax1(abs(result),abs(area)).le. + * defabs*0.1e-01) go to 130 + if(0.1e-01.gt.(result/area).or.(result/area).gt.0.1e+03. + *or.errsum.gt.abs(area)) ier = 6 + go to 130 +c +c compute global integral sum. +c + 115 result = 0.0e+00 + do 120 k = 1,last + result = result+rlist(k) + 120 continue + abserr = errsum + 130 neval = 30*last-15 + if(inf.eq.2) neval = 2*neval + if(ier.gt.2) ier=ier-1 + 999 return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qagp.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,223 @@ + subroutine qagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr, + * neval,ier,leniw,lenw,last,iwork,work) +c***begin prologue qagp +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a2a1 +c***keywords automatic integrator, general-purpose, +c singularities at user specified points, +c extrapolation, globally adaptive +c***author piessens,robert,appl. math. & progr. div - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose the routine calculates an approximation result to a given +c definite integral i = integral of f over (a,b), +c hopefully satisfying following claim for accuracy +c break points of the integration interval, where local +c difficulties of the integrand may occur(e.g. singularities, +c discontinuities), are provided by the user. +c***description +c +c computation of a definite integral +c standard fortran subroutine +c real version +c +c parameters +c on entry +c f - subroutine f(x,ierr,result) defining the integrand +c function f(x). the actual name for f needs to be +c declared e x t e r n a l in the driver program. +c +c a - real +c lower limit of integration +c +c b - real +c upper limit of integration +c +c npts2 - integer +c number equal to two more than the number of +c user-supplied break points within the integration +c range, npts.ge.2. +c if npts2.lt.2, the routine will end with ier = 6. +c +c points - real +c vector of dimension npts2, the first (npts2-2) +c elements of which are the user provided break +c points. if these points do not constitute an +c ascending sequence there will be an automatic +c sorting. +c +c epsabs - real +c absolute accuracy requested +c epsrel - real +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c on return +c result - real +c approximation to the integral +c +c abserr - real +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c ier.gt.0 abnormal termination of the routine. +c the estimates for integral and error are +c less reliable. it is assumed that the +c requested accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more +c subdivisions by increasing the value of +c limit (and taking the according dimension +c adjustments into account). however, if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. if +c the position of a local difficulty can be +c determined (i.e. singularity, +c discontinuity within the interval), it +c should be supplied to the routine as an +c element of the vector points. if necessary +c an appropriate special-purpose integrator +c must be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is +c detected, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour occurs +c at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. +c it is presumed that the requested +c tolerance cannot be achieved, and that +c the returned result is the best which +c can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier.gt.0. +c = 6 the input is invalid because +c npts2.lt.2 or +c break points are specified outside +c the integration range or +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) +c result, abserr, neval, last are set to +c zero. exept when leniw or lenw or npts2 is +c invalid, iwork(1), iwork(limit+1), +c work(limit*2+1) and work(limit*3+1) +c are set to zero. +c work(1) is set to a and work(limit+1) +c to b (where limit = (leniw-npts2)/2). +c +c dimensioning parameters +c leniw - integer +c dimensioning parameter for iwork +c leniw determines limit = (leniw-npts2)/2, +c which is the maximum number of subintervals in the +c partition of the given integration interval (a,b), +c leniw.ge.(3*npts2-2). +c if leniw.lt.(3*npts2-2), the routine will end with +c ier = 6. +c +c lenw - integer +c dimensioning parameter for work +c lenw must be at least leniw*2-npts2. +c if lenw.lt.leniw*2-npts2, the routine will end +c with ier = 6. +c +c last - integer +c on return, last equals the number of subintervals +c produced in the subdivision process, which +c determines the number of significant elements +c actually in the work arrays. +c +c work arrays +c iwork - integer +c vector of dimension at least leniw. on return, +c the first k elements of which contain +c pointers to the error estimates over the +c subintervals, such that work(limit*3+iwork(1)),..., +c work(limit*3+iwork(k)) form a decreasing +c sequence, with k = last if last.le.(limit/2+2), and +c k = limit+1-last otherwise +c iwork(limit+1), ...,iwork(limit+last) contain the +c subdivision levels of the subintervals, i.e. +c if (aa,bb) is a subinterval of (p1,p2) +c where p1 as well as p2 is a user-provided +c break point or integration limit, then (aa,bb) has +c level l if abs(bb-aa) = abs(p2-p1)*2**(-l), +c iwork(limit*2+1), ..., iwork(limit*2+npts2) have +c no significance for the user, +c note that limit = (leniw-npts2)/2. +c +c work - real +c vector of dimension at least lenw +c on return +c work(1), ..., work(last) contain the left +c end points of the subintervals in the +c partition of (a,b), +c work(limit+1), ..., work(limit+last) contain +c the right end points, +c work(limit*2+1), ..., work(limit*2+last) contain +c the integral approximations over the subintervals, +c work(limit*3+1), ..., work(limit*3+last) +c contain the corresponding error estimates, +c work(limit*4+1), ..., work(limit*4+npts2) +c contain the integration limits and the +c break points sorted in an ascending sequence. +c note that limit = (leniw-npts2)/2. +c +c***references (none) +c***routines called qagpe,xerror +c***end prologue qagp +c + real a,abserr,b,epsabs,epsrel,points,result,work + integer ier,iwork,leniw,lenw,limit,lvl,l1,l2,l3,neval,npts2 +c + dimension iwork(leniw),points(npts2),work(lenw) +c + external f +c +c check validity of limit and lenw. +c +c***first executable statement qagp + ier = 6 + neval = 0 + last = 0 + result = 0.0e+00 + abserr = 0.0e+00 + if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2) + * go to 10 +c +c prepare call for qagpe. +c + limit = (leniw-npts2)/2 + l1 = limit+1 + l2 = limit+l1 + l3 = limit+l2 + l4 = limit+l3 +c + call qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr, + * neval,ier,work(1),work(l1),work(l2),work(l3),work(l4), + * iwork(1),iwork(l1),iwork(l2),last) +c +c call error handler if necessary. +c + lvl = 0 +10 if(ier.eq.6) lvl = 1 + if(ier.ne.0) call xerror('abnormal return from qagp',26,ier,lvl) + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qagpe.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,560 @@ + subroutine qagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result, + * abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin, + * last) +c***begin prologue qagpe +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a2a1 +c***keywords automatic integrator, general-purpose, +c singularities at user specified points, +c extrapolation, globally adaptive. +c***author piessens,robert ,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose the routine calculates an approximation result to a given +c definite integral i = integral of f over (a,b),hopefully +c satisfying following claim for accuracy abs(i-result).le. +c max(epsabs,epsrel*abs(i)). break points of the integration +c interval, where local difficulties of the integrand may +c occur(e.g. singularities,discontinuities),provided by user. +c***description +c +c computation of a definite integral +c standard fortran subroutine +c real version +c +c parameters +c on entry +c f - subroutine f(x,ierr,result) defining the integrand +c function f(x). the actual name for f needs to be +c declared e x t e r n a l in the driver program. +c +c a - real +c lower limit of integration +c +c b - real +c upper limit of integration +c +c npts2 - integer +c number equal to two more than the number of +c user-supplied break points within the integration +c range, npts2.ge.2. +c if npts2.lt.2, the routine will end with ier = 6. +c +c points - real +c vector of dimension npts2, the first (npts2-2) +c elements of which are the user provided break +c points. if these points do not constitute an +c ascending sequence there will be an automatic +c sorting. +c +c epsabs - real +c absolute accuracy requested +c epsrel - real +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c limit - integer +c gives an upper bound on the number of subintervals +c in the partition of (a,b), limit.ge.npts2 +c if limit.lt.npts2, the routine will end with +c ier = 6. +c +c on return +c result - real +c approximation to the integral +c +c abserr - real +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c ier.gt.0 abnormal termination of the routine. +c the estimates for integral and error are +c less reliable. it is assumed that the +c requested accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more +c subdivisions by increasing the value of +c limit (and taking the according dimension +c adjustments into account). however, if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. if +c the position of a local difficulty can be +c determined (i.e. singularity, +c discontinuity within the interval), it +c should be supplied to the routine as an +c element of the vector points. if necessary +c an appropriate special-purpose integrator +c must be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is +c detected, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour occurs +c at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. it is presumed that +c the requested tolerance cannot be +c achieved, and that the returned result is +c the best which can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier.gt.0. +c = 6 the input is invalid because +c npts2.lt.2 or +c break points are specified outside +c the integration range or +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) +c or limit.lt.npts2. +c result, abserr, neval, last, rlist(1), +c and elist(1) are set to zero. alist(1) and +c blist(1) are set to a and b respectively. +c +c alist - real +c vector of dimension at least limit, the first +c last elements of which are the left end points +c of the subintervals in the partition of the given +c integration range (a,b) +c +c blist - real +c vector of dimension at least limit, the first +c last elements of which are the right end points +c of the subintervals in the partition of the given +c integration range (a,b) +c +c rlist - real +c vector of dimension at least limit, the first +c last elements of which are the integral +c approximations on the subintervals +c +c elist - real +c vector of dimension at least limit, the first +c last elements of which are the moduli of the +c absolute error estimates on the subintervals +c +c pts - real +c vector of dimension at least npts2, containing the +c integration limits and the break points of the +c interval in ascending sequence. +c +c level - integer +c vector of dimension at least limit, containing the +c subdivision levels of the subinterval, i.e. if +c (aa,bb) is a subinterval of (p1,p2) where p1 as +c well as p2 is a user-provided break point or +c integration limit, then (aa,bb) has level l if +c abs(bb-aa) = abs(p2-p1)*2**(-l). +c +c ndin - integer +c vector of dimension at least npts2, after first +c integration over the intervals (pts(i)),pts(i+1), +c i = 0,1, ..., npts2-2, the error estimates over +c some of the intervals may have been increased +c artificially, in order to put their subdivision +c forward. if this happens for the subinterval +c numbered k, ndin(k) is put to 1, otherwise +c ndin(k) = 0. +c +c iord - integer +c vector of dimension at least limit, the first k +c elements of which are pointers to the +c error estimates over the subintervals, +c such that elist(iord(1)), ..., elist(iord(k)) +c form a decreasing sequence, with k = last +c if last.le.(limit/2+2), and k = limit+1-last +c otherwise +c +c last - integer +c number of subintervals actually produced in the +c subdivisions process +c +c***references (none) +c***routines called qelg,qk21,qpsrt,r1mach +c***end prologue qagpe + real a,abseps,abserr,alist,area,area1,area12,area2,a1, + * a2,b,blist,b1,b2,correc,defabs,defab1,defab2, + * dres,r1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd, + * errmax,error1,erro12,error2,errsum,ertest,oflow,points,pts, + * resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp, + * uflow + integer i,id,ier,ierro,ind1,ind2,iord,ip1,iroff1,iroff2, + * iroff3,j,jlow,jupbnd,k,ksgn,ktmin,last,levcur,level,levmax, + * limit,maxerr,ndin,neval,nint,nintp1,npts,npts2,nres, + * nrmax,numrl2 + logical extrap,noext +c +c + dimension alist(limit),blist(limit),elist(limit),iord(limit), + * level(limit),ndin(npts2),points(npts2),pts(npts2),res3la(3), + * rlist(limit),rlist2(52) +c + external f +c +c the dimension of rlist2 is determined by the value of +c limexp in subroutine epsalg (rlist2 should be of dimension +c (limexp+2) at least). +c +c +c list of major variables +c ----------------------- +c +c alist - list of left end points of all subintervals +c considered up to now +c blist - list of right end points of all subintervals +c considered up to now +c rlist(i) - approximation to the integral over +c (alist(i),blist(i)) +c rlist2 - array of dimension at least limexp+2 +c containing the part of the epsilon table which +c is still needed for further computations +c elist(i) - error estimate applying to rlist(i) +c maxerr - pointer to the interval with largest error +c estimate +c errmax - elist(maxerr) +c erlast - error on the interval currently subdivided +c (before that subdivision has taken place) +c area - sum of the integrals over the subintervals +c errsum - sum of the errors over the subintervals +c errbnd - requested accuracy max(epsabs,epsrel* +c abs(result)) +c *****1 - variable for the left subinterval +c *****2 - variable for the right subinterval +c last - index for subdivision +c nres - number of calls to the extrapolation routine +c numrl2 - number of elements in rlist2. if an +c appropriate approximation to the compounded +c integral has been obtained, it is put in +c rlist2(numrl2) after numrl2 has been increased +c by one. +c erlarg - sum of the errors over the intervals larger +c than the smallest interval considered up to now +c extrap - logical variable denoting that the routine +c is attempting to perform extrapolation. i.e. +c before subdividing the smallest interval we +c try to decrease the value of erlarg. +c noext - logical variable denoting that extrapolation is +c no longer allowed (true-value) +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c oflow is the largest positive magnitude. +c +c***first executable statement qagpe + epmach = r1mach(4) +c +c test on validity of parameters +c ----------------------------- +c + ier = 0 + neval = 0 + last = 0 + result = 0.0e+00 + abserr = 0.0e+00 + alist(1) = a + blist(1) = b + rlist(1) = 0.0e+00 + elist(1) = 0.0e+00 + iord(1) = 0 + level(1) = 0 + npts = npts2-2 + if(npts2.lt.2.or.limit.le.npts.or.(epsabs.le.0.0e+00.and. + * epsrel.lt.amax1(0.5e+02*epmach,0.5e-14))) ier = 6 + if(ier.eq.6) go to 210 +c +c if any break points are provided, sort them into an +c ascending sequence. +c + sign = 1.0e+00 + if(a.gt.b) sign = -1.0e+00 + pts(1) = amin1(a,b) + if(npts.eq.0) go to 15 + do 10 i = 1,npts + pts(i+1) = points(i) + 10 continue + 15 pts(npts+2) = amax1(a,b) + nint = npts+1 + a1 = pts(1) + if(npts.eq.0) go to 40 + nintp1 = nint+1 + do 20 i = 1,nint + ip1 = i+1 + do 20 j = ip1,nintp1 + if(pts(i).le.pts(j)) go to 20 + temp = pts(i) + pts(i) = pts(j) + pts(j) = temp + 20 continue + if(pts(1).ne.amin1(a,b).or.pts(nintp1).ne. + * amax1(a,b)) ier = 6 + if(ier.eq.6) go to 999 +c +c compute first integral and error approximations. +c ------------------------------------------------ +c + 40 resabs = 0.0e+00 + do 50 i = 1,nint + b1 = pts(i+1) + call qk21(f,a1,b1,area1,error1,defabs,resa,ier) + if (ier.lt.0) return + abserr = abserr+error1 + result = result+area1 + ndin(i) = 0 + if(error1.eq.resa.and.error1.ne.0.0e+00) ndin(i) = 1 + resabs = resabs+defabs + level(i) = 0 + elist(i) = error1 + alist(i) = a1 + blist(i) = b1 + rlist(i) = area1 + iord(i) = i + a1 = b1 + 50 continue + errsum = 0.0e+00 + do 55 i = 1,nint + if(ndin(i).eq.1) elist(i) = abserr + errsum = errsum+elist(i) + 55 continue +c +c test on accuracy. +c + last = nint + neval = 21*nint + dres = abs(result) + errbnd = amax1(epsabs,epsrel*dres) + if(abserr.le.0.1e+03*epmach*resabs.and.abserr.gt. + * errbnd) ier = 2 + if(nint.eq.1) go to 80 + do 70 i = 1,npts + jlow = i+1 + ind1 = iord(i) + do 60 j = jlow,nint + ind2 = iord(j) + if(elist(ind1).gt.elist(ind2)) go to 60 + ind1 = ind2 + k = j + 60 continue + if(ind1.eq.iord(i)) go to 70 + iord(k) = iord(i) + iord(i) = ind1 + 70 continue + if(limit.lt.npts2) ier = 1 + 80 if(ier.ne.0.or.abserr.le.errbnd) go to 999 +c +c initialization +c -------------- +c + rlist2(1) = result + maxerr = iord(1) + errmax = elist(maxerr) + area = result + nrmax = 1 + nres = 0 + numrl2 = 1 + ktmin = 0 + extrap = .false. + noext = .false. + erlarg = errsum + ertest = errbnd + levmax = 1 + iroff1 = 0 + iroff2 = 0 + iroff3 = 0 + ierro = 0 + uflow = r1mach(1) + oflow = r1mach(2) + abserr = oflow + ksgn = -1 + if(dres.ge.(0.1e+01-0.5e+02*epmach)*resabs) ksgn = 1 +c +c main do-loop +c ------------ +c + do 160 last = npts2,limit +c +c bisect the subinterval with the nrmax-th largest +c error estimate. +c + levcur = level(maxerr)+1 + a1 = alist(maxerr) + b1 = 0.5e+00*(alist(maxerr)+blist(maxerr)) + a2 = b1 + b2 = blist(maxerr) + erlast = errmax + call qk21(f,a1,b1,area1,error1,resa,defab1,ier) + if (ier.lt.0) return + call qk21(f,a2,b2,area2,error2,resa,defab2,ier) + if (ier.lt.0) return +c +c improve previous approximations to integral +c and error and test for accuracy. +c + neval = neval+42 + area12 = area1+area2 + erro12 = error1+error2 + errsum = errsum+erro12-errmax + area = area+area12-rlist(maxerr) + if(defab1.eq.error1.or.defab2.eq.error2) go to 95 + if(abs(rlist(maxerr)-area12).gt.0.1e-04*abs(area12) + * .or.erro12.lt.0.99e+00*errmax) go to 90 + if(extrap) iroff2 = iroff2+1 + if(.not.extrap) iroff1 = iroff1+1 + 90 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 + 95 level(maxerr) = levcur + level(last) = levcur + rlist(maxerr) = area1 + rlist(last) = area2 + errbnd = amax1(epsabs,epsrel*abs(area)) +c +c test for roundoff error and eventually +c set error flag. +c + if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 + if(iroff2.ge.5) ierro = 3 +c +c set error flag in the case that the number of +c subintervals equals limit. +c + if(last.eq.limit) ier = 1 +c +c set error flag in the case of bad integrand behaviour +c at a point of the integration range +c + if(amax1(abs(a1),abs(b2)).le.(0.1e+01+0.1e+03*epmach)* + * (abs(a2)+0.1e+04*uflow)) ier = 4 +c +c append the newly-created intervals to the list. +c + if(error2.gt.error1) go to 100 + alist(last) = a2 + blist(maxerr) = b1 + blist(last) = b2 + elist(maxerr) = error1 + elist(last) = error2 + go to 110 + 100 alist(maxerr) = a2 + alist(last) = a1 + blist(last) = b1 + rlist(maxerr) = area2 + rlist(last) = area1 + elist(maxerr) = error2 + elist(last) = error1 +c +c call subroutine qpsrt to maintain the descending ordering +c in the list of error estimates and select the +c subinterval with nrmax-th largest error estimate (to be +c bisected next). +c + 110 call qpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) +c ***jump out of do-loop + if(errsum.le.errbnd) go to 190 +c ***jump out of do-loop + if(ier.ne.0) go to 170 + if(noext) go to 160 + erlarg = erlarg-erlast + if(levcur+1.le.levmax) erlarg = erlarg+erro12 + if(extrap) go to 120 +c +c test whether the interval to be bisected next is the +c smallest interval. +c + if(level(maxerr)+1.le.levmax) go to 160 + extrap = .true. + nrmax = 2 + 120 if(ierro.eq.3.or.erlarg.le.ertest) go to 140 +c +c the smallest interval has the largest error. +c before bisecting decrease the sum of the errors +c over the larger intervals (erlarg) and perform +c extrapolation. +c + id = nrmax + jupbnd = last + if(last.gt.(2+limit/2)) jupbnd = limit+3-last + do 130 k = id,jupbnd + maxerr = iord(nrmax) + errmax = elist(maxerr) +c ***jump out of do-loop + if(level(maxerr)+1.le.levmax) go to 160 + nrmax = nrmax+1 + 130 continue +c +c perform extrapolation. +c + 140 numrl2 = numrl2+1 + rlist2(numrl2) = area + if(numrl2.le.2) go to 155 + call qelg(numrl2,rlist2,reseps,abseps,res3la,nres) + ktmin = ktmin+1 + if(ktmin.gt.5.and.abserr.lt.0.1e-02*errsum) ier = 5 + if(abseps.ge.abserr) go to 150 + ktmin = 0 + abserr = abseps + result = reseps + correc = erlarg + ertest = amax1(epsabs,epsrel*abs(reseps)) +c ***jump out of do-loop + if(abserr.lt.ertest) go to 170 +c +c prepare bisection of the smallest interval. +c + 150 if(numrl2.eq.1) noext = .true. + if(ier.ge.5) go to 170 + 155 maxerr = iord(1) + errmax = elist(maxerr) + nrmax = 1 + extrap = .false. + levmax = levmax+1 + erlarg = errsum + 160 continue +c +c set the final result. +c --------------------- +c +c + 170 if(abserr.eq.oflow) go to 190 + if((ier+ierro).eq.0) go to 180 + if(ierro.eq.3) abserr = abserr+correc + if(ier.eq.0) ier = 3 + if(result.ne.0.0e+00.and.area.ne.0.0e+00)go to 175 + if(abserr.gt.errsum)go to 190 + if(area.eq.0.0e+00) go to 210 + go to 180 + 175 if(abserr/abs(result).gt.errsum/abs(area))go to 190 +c +c test on divergence. +c + 180 if(ksgn.eq.(-1).and.amax1(abs(result),abs(area)).le. + * resabs*0.1e-01) go to 210 + if(0.1e-01.gt.(result/area).or.(result/area).gt.0.1e+03.or. + * errsum.gt.abs(area)) ier = 6 + go to 210 +c +c compute global integral sum. +c + 190 result = 0.0e+00 + do 200 k = 1,last + result = result+rlist(k) + 200 continue + abserr = errsum + 210 if(ier.gt.2) ier = ier - 1 + result = result*sign + 999 return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qelg.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,184 @@ + subroutine qelg(n,epstab,result,abserr,res3la,nres) +c***begin prologue qelg +c***refer to qagie,qagoe,qagpe,qagse +c***routines called r1mach +c***revision date 830518 (yymmdd) +c***keywords epsilon algorithm, convergence acceleration, +c extrapolation +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math & progr. div. - k.u.leuven +c***purpose the routine determines the limit of a given sequence of +c approximations, by means of the epsilon algorithm of +c p. wynn. an estimate of the absolute error is also given. +c the condensed epsilon table is computed. only those +c elements needed for the computation of the next diagonal +c are preserved. +c***description +c +c epsilon algorithm +c standard fortran subroutine +c real version +c +c parameters +c n - integer +c epstab(n) contains the new element in the +c first column of the epsilon table. +c +c epstab - real +c vector of dimension 52 containing the elements +c of the two lower diagonals of the triangular +c epsilon table. the elements are numbered +c starting at the right-hand corner of the +c triangle. +c +c result - real +c resulting approximation to the integral +c +c abserr - real +c estimate of the absolute error computed from +c result and the 3 previous results +c +c res3la - real +c vector of dimension 3 containing the last 3 +c results +c +c nres - integer +c number of calls to the routine +c (should be zero at first call) +c +c***end prologue qelg +c + real abserr,delta1,delta2,delta3,r1mach, + * epmach,epsinf,epstab,error,err1,err2,err3,e0,e1,e1abs,e2,e3, + * oflow,res,result,res3la,ss,tol1,tol2,tol3 + integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,n,newelm,nres,num + dimension epstab(52),res3la(3) +c +c list of major variables +c ----------------------- +c +c e0 - the 4 elements on which the +c e1 computation of a new element in +c e2 the epsilon table is based +c e3 e0 +c e3 e1 new +c e2 +c newelm - number of elements to be computed in the new +c diagonal +c error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2) +c result - the element in the new diagonal with least value +c of error +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c oflow is the largest positive magnitude. +c limexp is the maximum number of elements the epsilon +c table can contain. if this number is reached, the upper +c diagonal of the epsilon table is deleted. +c +c***first executable statement qelg + epmach = r1mach(4) + oflow = r1mach(2) + nres = nres+1 + abserr = oflow + result = epstab(n) + if(n.lt.3) go to 100 + limexp = 50 + epstab(n+2) = epstab(n) + newelm = (n-1)/2 + epstab(n) = oflow + num = n + k1 = n + do 40 i = 1,newelm + k2 = k1-1 + k3 = k1-2 + res = epstab(k1+2) + e0 = epstab(k3) + e1 = epstab(k2) + e2 = res + e1abs = abs(e1) + delta2 = e2-e1 + err2 = abs(delta2) + tol2 = amax1(abs(e2),e1abs)*epmach + delta3 = e1-e0 + err3 = abs(delta3) + tol3 = amax1(e1abs,abs(e0))*epmach + if(err2.gt.tol2.or.err3.gt.tol3) go to 10 +c +c if e0, e1 and e2 are equal to within machine +c accuracy, convergence is assumed. +c result = e2 +c abserr = abs(e1-e0)+abs(e2-e1) +c + result = res + abserr = err2+err3 +c ***jump out of do-loop + go to 100 + 10 e3 = epstab(k1) + epstab(k1) = e1 + delta1 = e1-e3 + err1 = abs(delta1) + tol1 = amax1(e1abs,abs(e3))*epmach +c +c if two elements are very close to each other, omit +c a part of the table by adjusting the value of n +c + if(err1.le.tol1.or.err2.le.tol2.or.err3.le.tol3) go to 20 + ss = 0.1e+01/delta1+0.1e+01/delta2-0.1e+01/delta3 + epsinf = abs(ss*e1) +c +c test to detect irregular behaviour in the table, and +c eventually omit a part of the table adjusting the value +c of n. +c + if(epsinf.gt.0.1e-03) go to 30 + 20 n = i+i-1 +c ***jump out of do-loop + go to 50 +c +c compute a new element and eventually adjust +c the value of result. +c + 30 res = e1+0.1e+01/ss + epstab(k1) = res + k1 = k1-2 + error = err2+abs(res-e2)+err3 + if(error.gt.abserr) go to 40 + abserr = error + result = res + 40 continue +c +c shift the table. +c + 50 if(n.eq.limexp) n = 2*(limexp/2)-1 + ib = 1 + if((num/2)*2.eq.num) ib = 2 + ie = newelm+1 + do 60 i=1,ie + ib2 = ib+2 + epstab(ib) = epstab(ib2) + ib = ib2 + 60 continue + if(num.eq.n) go to 80 + indx = num-n+1 + do 70 i = 1,n + epstab(i)= epstab(indx) + indx = indx+1 + 70 continue + 80 if(nres.ge.4) go to 90 + res3la(nres) = result + abserr = oflow + go to 100 +c +c compute error estimate +c + 90 abserr = abs(result-res3la(3))+abs(result-res3la(2)) + * +abs(result-res3la(1)) + res3la(1) = res3la(2) + res3la(2) = res3la(3) + res3la(3) = result + 100 abserr = amax1(abserr,0.5e+01*epmach*abs(result)) + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qk15i.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,202 @@ + subroutine qk15i(f,boun,inf,a,b,result,abserr,resabs,resasc,ierr) +c***begin prologue qk15i +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a3a2,h2a4a2 +c***keywords 15-point transformed gauss-kronrod rules +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose the original (infinite integration range is mapped +c onto the interval (0,1) and (a,b) is a part of (0,1). +c it is the purpose to compute +c i = integral of transformed integrand over (a,b), +c j = integral of abs(transformed integrand) over (a,b). +c***description +c +c integration rule +c standard fortran subroutine +c real version +c +c parameters +c on entry +c f - subroutine f(x,ierr,result) defining the integrand +c function f(x). the actual name for f needs to be +c declared e x t e r n a l in the calling program. +c +c boun - real +c finite bound of original integration +c range (set to zero if inf = +2) +c +c inf - integer +c if inf = -1, the original interval is +c (-infinity,bound), +c if inf = +1, the original interval is +c (bound,+infinity), +c if inf = +2, the original interval is +c (-infinity,+infinity) and +c the integral is computed as the sum of two +c integrals, one over (-infinity,0) and one over +c (0,+infinity). +c +c a - real +c lower limit for integration over subrange +c of (0,1) +c +c b - real +c upper limit for integration over subrange +c of (0,1) +c +c on return +c result - real +c approximation to the integral i +c result is computed by applying the 15-point +c kronrod rule(resk) obtained by optimal addition +c of abscissae to the 7-point gauss rule(resg). +c +c abserr - real +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c resabs - real +c approximation to the integral j +c +c resasc - real +c approximation to the integral of +c abs((transformed integrand)-i/(b-a)) over (a,b) +c +c***references (none) +c***routines called r1mach +c***end prologue qk15i +c + real a,absc,absc1,absc2,abserr,b,boun,centr, + * dinf,r1mach,epmach,fc,fsum,fval1,fval2,fvalt,fv1, + * fv2,hlgth,resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2, + * uflow,wg,wgk,xgk + integer inf,j,min0 + external f +c + dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8) +c +c the abscissae and weights are supplied for the interval +c (-1,1). because of symmetry only the positive abscissae and +c their corresponding weights are given. +c +c xgk - abscissae of the 15-point kronrod rule +c xgk(2), xgk(4), ... abscissae of the 7-point +c gauss rule +c xgk(1), xgk(3), ... abscissae which are optimally +c added to the 7-point gauss rule +c +c wgk - weights of the 15-point kronrod rule +c +c wg - weights of the 7-point gauss rule, corresponding +c to the abscissae xgk(2), xgk(4), ... +c wg(1), wg(3), ... are set to zero. +c + data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7), + * xgk(8)/ + * 0.9914553711208126e+00, 0.9491079123427585e+00, + * 0.8648644233597691e+00, 0.7415311855993944e+00, + * 0.5860872354676911e+00, 0.4058451513773972e+00, + * 0.2077849550078985e+00, 0.0000000000000000e+00/ +c + data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7), + * wgk(8)/ + * 0.2293532201052922e-01, 0.6309209262997855e-01, + * 0.1047900103222502e+00, 0.1406532597155259e+00, + * 0.1690047266392679e+00, 0.1903505780647854e+00, + * 0.2044329400752989e+00, 0.2094821410847278e+00/ +c + data wg(1),wg(2),wg(3),wg(4),wg(5),wg(6),wg(7),wg(8)/ + * 0.0000000000000000e+00, 0.1294849661688697e+00, + * 0.0000000000000000e+00, 0.2797053914892767e+00, + * 0.0000000000000000e+00, 0.3818300505051189e+00, + * 0.0000000000000000e+00, 0.4179591836734694e+00/ +c +c +c list of major variables +c ----------------------- +c +c centr - mid point of the interval +c hlgth - half-length of the interval +c absc* - abscissa +c tabsc* - transformed abscissa +c fval* - function value +c resg - result of the 7-point gauss formula +c resk - result of the 15-point kronrod formula +c reskh - approximation to the mean value of the transformed +c integrand over (a,b), i.e. to i/(b-a) +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c +c***first executable statement qk15i + epmach = r1mach(4) + uflow = r1mach(1) + dinf = min0(1,inf) +c + centr = 0.5e+00*(a+b) + hlgth = 0.5e+00*(b-a) + tabsc1 = boun+dinf*(0.1e+01-centr)/centr + call f(tabsc1, ierr, fval1) + if (ierr.lt.0) return + if(inf.eq.2) then + call f(-tabsc1, ierr, fval1) + if (ierr.lt.0) return + fval1 = fval1 + fvalt + endif + fc = (fval1/centr)/centr +c +c compute the 15-point kronrod approximation to +c the integral, and estimate the error. +c + resg = wg(8)*fc + resk = wgk(8)*fc + resabs = abs(resk) + do 10 j=1,7 + absc = hlgth*xgk(j) + absc1 = centr-absc + absc2 = centr+absc + tabsc1 = boun+dinf*(0.1e+01-absc1)/absc1 + tabsc2 = boun+dinf*(0.1e+01-absc2)/absc2 + call f(tabsc1, ierr, fval1) + if (ierr.lt.0) return + call f(tabsc2, ierr, fval2) + if (ierr.lt.0) return + if(inf.eq.2) then + call f(-tabsc1,ierr,fvalt) + if (ierr.lt.0) return + fval1 = fval1 + fvalt; + endif + if(inf.eq.2) then + call f(-tabsc2,ierr,fvalt) + if (ierr.lt.0) return + fval2 = fval2 + fvalt; + endif + fval1 = (fval1/absc1)/absc1 + fval2 = (fval2/absc2)/absc2 + fv1(j) = fval1 + fv2(j) = fval2 + fsum = fval1+fval2 + resg = resg+wg(j)*fsum + resk = resk+wgk(j)*fsum + resabs = resabs+wgk(j)*(abs(fval1)+abs(fval2)) + 10 continue + reskh = resk*0.5e+00 + resasc = wgk(8)*abs(fc-reskh) + do 20 j=1,7 + resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) + 20 continue + result = resk*hlgth + resasc = resasc*hlgth + resabs = resabs*hlgth + abserr = abs((resk-resg)*hlgth) + if(resasc.ne.0.0e+00.and.abserr.ne.0.e0) abserr = resasc* + * amin1(0.1e+01,(0.2e+03*abserr/resasc)**1.5e+00) + if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1 + * ((epmach*0.5e+02)*resabs,abserr) + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qk21.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,175 @@ + subroutine qk21(f,a,b,result,abserr,resabs,resasc,ierr) +c***begin prologue qk21 +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a1a2 +c***keywords 21-point gauss-kronrod rules +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose to compute i = integral of f over (a,b), with error +c estimate +c j = integral of abs(f) over (a,b) +c***description +c +c integration rules +c standard fortran subroutine +c real version +c +c parameters +c on entry +c f - subroutine f(x,ierr,result) defining the integrand +c function f(x). the actual name for f needs to be +c declared e x t e r n a l in the driver program. +c +c a - real +c lower limit of integration +c +c b - real +c upper limit of integration +c +c on return +c result - real +c approximation to the integral i +c result is computed by applying the 21-point +c kronrod rule (resk) obtained by optimal addition +c of abscissae to the 10-point gauss rule (resg). +c +c abserr - real +c estimate of the modulus of the absolute error, +c which should not exceed abs(i-result) +c +c resabs - real +c approximation to the integral j +c +c resasc - real +c approximation to the integral of abs(f-i/(b-a)) +c over (a,b) +c +c***references (none) +c***routines called r1mach +c***end prologue qk21 +c + real a,absc,abserr,b,centr,dhlgth,epmach,fc,fsum,fval1,fval2, + * fv1,fv2,hlgth,resabs,resg,resk,reskh,result,r1mach,uflow,wg,wgk, + * xgk + integer j,jtw,jtwm1 + external f +c + dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11) +c +c the abscissae and weights are given for the interval (-1,1). +c because of symmetry only the positive abscissae and their +c corresponding weights are given. +c +c xgk - abscissae of the 21-point kronrod rule +c xgk(2), xgk(4), ... abscissae of the 10-point +c gauss rule +c xgk(1), xgk(3), ... abscissae which are optimally +c added to the 10-point gauss rule +c +c wgk - weights of the 21-point kronrod rule +c +c wg - weights of the 10-point gauss rule +c + data xgk(1),xgk(2),xgk(3),xgk(4),xgk(5),xgk(6),xgk(7), + * xgk(8),xgk(9),xgk(10),xgk(11)/ + * 0.9956571630258081e+00, 0.9739065285171717e+00, + * 0.9301574913557082e+00, 0.8650633666889845e+00, + * 0.7808177265864169e+00, 0.6794095682990244e+00, + * 0.5627571346686047e+00, 0.4333953941292472e+00, + * 0.2943928627014602e+00, 0.1488743389816312e+00, + * 0.0000000000000000e+00/ +c + data wgk(1),wgk(2),wgk(3),wgk(4),wgk(5),wgk(6),wgk(7), + * wgk(8),wgk(9),wgk(10),wgk(11)/ + * 0.1169463886737187e-01, 0.3255816230796473e-01, + * 0.5475589657435200e-01, 0.7503967481091995e-01, + * 0.9312545458369761e-01, 0.1093871588022976e+00, + * 0.1234919762620659e+00, 0.1347092173114733e+00, + * 0.1427759385770601e+00, 0.1477391049013385e+00, + * 0.1494455540029169e+00/ +c + data wg(1),wg(2),wg(3),wg(4),wg(5)/ + * 0.6667134430868814e-01, 0.1494513491505806e+00, + * 0.2190863625159820e+00, 0.2692667193099964e+00, + * 0.2955242247147529e+00/ +c +c +c list of major variables +c ----------------------- +c +c centr - mid point of the interval +c hlgth - half-length of the interval +c absc - abscissa +c fval* - function value +c resg - result of the 10-point gauss formula +c resk - result of the 21-point kronrod formula +c reskh - approximation to the mean value of f over (a,b), +c i.e. to i/(b-a) +c +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c +c***first executable statement qk21 + epmach = r1mach(4) + uflow = r1mach(1) +c + centr = 0.5e+00*(a+b) + hlgth = 0.5e+00*(b-a) + dhlgth = abs(hlgth) +c +c compute the 21-point kronrod approximation to +c the integral, and estimate the absolute error. +c + resg = 0.0e+00 + call f(centr, ierr, fc) + if (ierr .lt. 0) return + resk = wgk(11)*fc + resabs = abs(resk) + do 10 j=1,5 + jtw = 2*j + absc = hlgth*xgk(jtw) + call f(centr-absc,ierr,fval1) + if (ierr .lt. 0) return + call f(centr+absc,ierr,fval2) + if (ierr .lt. 0) return + fv1(jtw) = fval1 + fv2(jtw) = fval2 + fsum = fval1+fval2 + resg = resg+wg(j)*fsum + resk = resk+wgk(jtw)*fsum + resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2)) + 10 continue + do 15 j = 1,5 + jtwm1 = 2*j-1 + absc = hlgth*xgk(jtwm1) + call f(centr-absc,ierr,fval1) + if (ierr .lt. 0) return + call f(centr+absc,ierr,fval2) + if (ierr .lt. 0) return + fv1(jtwm1) = fval1 + fv2(jtwm1) = fval2 + fsum = fval1+fval2 + resk = resk+wgk(jtwm1)*fsum + resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2)) + 15 continue + reskh = resk*0.5e+00 + resasc = wgk(11)*abs(fc-reskh) + do 20 j=1,10 + resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) + 20 continue + result = resk*hlgth + resabs = resabs*dhlgth + resasc = resasc*dhlgth + abserr = abs((resk-resg)*hlgth) + if(resasc.ne.0.0e+00.and.abserr.ne.0.0e+00) + * abserr = resasc*amin1(0.1e+01, + * (0.2e+03*abserr/resasc)**1.5e+00) + if(resabs.gt.uflow/(0.5e+02*epmach)) abserr = amax1 + * ((epmach*0.5e+02)*resabs,abserr) + return + end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/quadpack/qpsrt.f Sun May 11 22:51:50 2008 +0200 @@ -0,0 +1,136 @@ + subroutine qpsrt(limit,last,maxerr,ermax,elist,iord,nrmax) +c***begin prologue qpsrt +c***refer to qage,qagie,qagpe,qagse,qawce,qawse,qawoe +c***routines called (none) +c***keywords sequential sorting +c***description +c +c 1. qpsrt +c ordering routine +c standard fortran subroutine +c real version +c +c 2. purpose +c this routine maintains the descending ordering +c in the list of the local error estimates resulting from +c the interval subdivision process. at each call two error +c estimates are inserted using the sequential search +c method, top-down for the largest error estimate +c and bottom-up for the smallest error estimate. +c +c 3. calling sequence +c call qpsrt(limit,last,maxerr,ermax,elist,iord,nrmax) +c +c parameters (meaning at output) +c limit - integer +c maximum number of error estimates the list +c can contain +c +c last - integer +c number of error estimates currently +c in the list +c +c maxerr - integer +c maxerr points to the nrmax-th largest error +c estimate currently in the list +c +c ermax - real +c nrmax-th largest error estimate +c ermax = elist(maxerr) +c +c elist - real +c vector of dimension last containing +c the error estimates +c +c iord - integer +c vector of dimension last, the first k +c elements of which contain pointers +c to the error estimates, such that +c elist(iord(1)),... , elist(iord(k)) +c form a decreasing sequence, with +c k = last if last.le.(limit/2+2), and +c k = limit+1-last otherwise +c +c nrmax - integer +c maxerr = iord(nrmax) +c +c 4. no subroutines or functions needed +c***end prologue qpsrt +c + real elist,ermax,errmax,errmin + integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr, + * nrmax + dimension elist(last),iord(last) +c +c check whether the list contains more than +c two error estimates. +c +c***first executable statement qpsrt + if(last.gt.2) go to 10 + iord(1) = 1 + iord(2) = 2 + go to 90 +c +c this part of the routine is only executed +c if, due to a difficult integrand, subdivision +c increased the error estimate. in the normal case +c the insert procedure should start after the +c nrmax-th largest error estimate. +c + 10 errmax = elist(maxerr) + if(nrmax.eq.1) go to 30 + ido = nrmax-1 + do 20 i = 1,ido + isucc = iord(nrmax-1) +c ***jump out of do-loop + if(errmax.le.elist(isucc)) go to 30 + iord(nrmax) = isucc + nrmax = nrmax-1 + 20 continue +c +c compute the number of elements in the list to +c be maintained in descending order. this number +c depends on the number of subdivisions still +c allowed. +c + 30 jupbn = last + if(last.gt.(limit/2+2)) jupbn = limit+3-last + errmin = elist(last) +c +c insert errmax by traversing the list top-down, +c starting comparison from the element elist(iord(nrmax+1)). +c + jbnd = jupbn-1 + ibeg = nrmax+1 + if(ibeg.gt.jbnd) go to 50 + do 40 i=ibeg,jbnd + isucc = iord(i) +c ***jump out of do-loop + if(errmax.ge.elist(isucc)) go to 60 + iord(i-1) = isucc + 40 continue + 50 iord(jbnd) = maxerr + iord(jupbn) = last + go to 90 +c +c insert errmin by traversing the list bottom-up. +c + 60 iord(i-1) = maxerr + k = jbnd + do 70 j=i,jbnd + isucc = iord(k) +c ***jump out of do-loop + if(errmin.lt.elist(isucc)) go to 80 + iord(k+1) = isucc + k = k-1 + 70 continue + iord(i) = last + go to 90 + 80 iord(k+1) = last +c +c set maxerr and ermax. +c + 90 maxerr = iord(nrmax) + ermax = elist(maxerr) + return + end
--- a/liboctave/ChangeLog Wed May 21 09:36:46 2008 -0400 +++ b/liboctave/ChangeLog Sun May 11 22:51:50 2008 +0200 @@ -1,3 +1,48 @@ +2008-05-21 David Bateman <dbateman@free.fr> + + * CmplxGEBAL.cc (ComplexGEPBALANCE), dbleGEPBAL.cc (GEPBALANCE), + fCmplxGEPBAL.cc (FloatComplexGEPBALANCE), floatGEPBAL.cc + (FloatGEPBALANCE): New class for generalized eigenvalue balancing. + * CmplxGEBAL.h (ComplexGEPBALANCE), dbleGEPBAL.h (GEPBALANCE), + fCmplxGEPBAL.h (FloatComplexGEPBALANCE), floatGEPBAL.h + (FloatGEPBALANCE): Declare them. + * Makefile.in (MATRIX_INC): Include them here. + (MATRIX_SRC): and here. + + * floatAEPBAL.cc (FloatAEPBALANCE), fCmplxAEPBAL.cc + (FloatComplexAEPBALANCE): New classes for single precision + Algebraic eignvalue balancing. + * floatAEPBAL.h (FloatAEPBALANCE), fCmplxAEPBAL.h + (FloatComplexAEPBALANCE): Declare them. + * Makefile.in (MATRIX_INC): Include them here. + (MATRIX_SRC): and here. + + * floatHESS.cc (FloatHESS), fCmplxHESS.cc (FloatComplexHESS): New + classes for single precision Hessenberg decomposition. + * floatHESS.h (FloatHESS), fCmplxHESS.h (FloatComplexHESS): + Declare them. + * Makefile.in (MATRIX_INC): Include them here. + (MATRIX_SRC): and here. + + * floatQR.cc (FloatQR), fCmplxQR.cc (FloatComplexQR): New + classes for single precision QR decomposition. + * floatQR.h (FloatQR), fCmplxQR.h (FloatComplexQR): + Declare them. + * Makefile.in (MATRIX_INC): Include them here. + (MATRIX_SRC): and here. + + * floatQRP.cc (FloatQRP), fCmplxQRP.cc (FloatComplexQRP): New + classes for single precision permuted QR decomposition. + * floatQRP.h (FloatQRP), fCmplxQRP.h (FloatComplexQRP): + Declare them. + * Makefile.in (MATRIX_INC): Include them here. + (MATRIX_SRC): and here. + + * mx-defs (FloatAEPBALANCE, FloatComplexAEPBALANCE, + ComplexGEPBALANCE, FloatGEPBALANCE,FloatComplexGEPBALANCE, + FloatHESS, FloatComplexHESS, FloatQR, FloatComplexQR, QRP, + ComplexQRP, FloatQRP, FloatComplexQRP): Declare classes. + 2008-05-20 David Bateman <dbateman@free.fr> * Array.cc (Array<T> Array<T>::transpose () const): Modify for tiled