Mercurial > octave-nkf
changeset 3983:7a37caf6ed43
[project @ 2002-07-11 02:36:25 by jwe]
| author | jwe |
|---|---|
| date | Thu, 11 Jul 2002 02:37:44 +0000 |
| parents | 334b3ee09f4b |
| children | addebffd4961 |
| files | ChangeLog configure.in libcruft/ChangeLog libcruft/Makefile.in libcruft/odessa/Makefile.in libcruft/odessa/odessa.f libcruft/odessa/odessa_addx.f libcruft/odessa/odessa_cfode.f libcruft/odessa/odessa_ewset.f libcruft/odessa/odessa_intdy.f libcruft/odessa/odessa_prepd.f libcruft/odessa/odessa_prepj.f libcruft/odessa/odessa_rscom.f libcruft/odessa/odessa_solsy.f libcruft/odessa/odessa_sprim.f libcruft/odessa/odessa_stesa.f libcruft/odessa/odessa_stode.f libcruft/odessa/odessa_svcom.f libcruft/odessa/odessa_vnorm.f libcruft/odessa/odessa_xsetf.f libcruft/odessa/xerr.f libcruft/odessa/xsetun.f |
| diffstat | 22 files changed, 3445 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog Wed Jul 10 19:23:58 2002 +0000 +++ b/ChangeLog Thu Jul 11 02:37:44 2002 +0000 @@ -1,3 +1,8 @@ +2002-07-10 John W. Eaton <jwe@bevo.che.wisc.edu> + + * configure.in (AC_CONFIG_FILES): Add libcruft/odessa/Makefile to + the list. + 2002-05-24 John W. Eaton <jwe@bevo.che.wisc.edu> * configure.in: Maybe add -fno-coalesce-templates to XTRA_CXXFLAGS
--- a/configure.in Wed Jul 10 19:23:58 2002 +0000 +++ b/configure.in Thu Jul 11 02:37:44 2002 +0000 @@ -22,7 +22,7 @@ ### 02111-1307, USA. AC_INIT -AC_REVISION($Revision: 1.361 $) +AC_REVISION($Revision: 1.362 $) AC_PREREQ(2.52) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -1263,7 +1263,8 @@ libcruft/dassl/Makefile libcruft/fftpack/Makefile \ libcruft/lapack/Makefile libcruft/linpack/Makefile \ libcruft/minpack/Makefile libcruft/misc/Makefile \ - libcruft/odepack/Makefile libcruft/ordered-qz/Makefile \ + libcruft/odepack/Makefile libcruft/odessa/Makefile \ + libcruft/ordered-qz/Makefile \ libcruft/quadpack/Makefile libcruft/ranlib/Makefile \ libcruft/slatec-fn/Makefile libcruft/slatec-err/Makefile \ libcruft/villad/Makefile \
--- a/libcruft/ChangeLog Wed Jul 10 19:23:58 2002 +0000 +++ b/libcruft/ChangeLog Thu Jul 11 02:37:44 2002 +0000 @@ -1,3 +1,8 @@ +2002-07-10 John W. Eaton <jwe@bevo.che.wisc.edu> + + * odessa: New subdirectory. + * Makefile.in (CRUFT_DIRS): Add it to the list. + 2002-06-27 John W. Eaton <jwe@bevo.che.wisc.edu> * slatec-err/xermsg.f (XERMSG): If MAXMES .LT. 0, messages may be
--- a/libcruft/Makefile.in Wed Jul 10 19:23:58 2002 +0000 +++ b/libcruft/Makefile.in Thu Jul 11 02:37:44 2002 +0000 @@ -29,7 +29,7 @@ # dirname is substituted by configure and may be the empty string. CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dassl @FFT_DIR@ @LAPACK_DIR@ \ - lapack-xtra linpack minpack misc odepack ordered-qz quadpack \ + lapack-xtra linpack minpack misc odepack odessa ordered-qz quadpack \ ranlib slatec-err slatec-fn villad SUBDIRS = $(CRUFT_DIRS)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/Makefile.in Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,19 @@ +# +# Makefile for octave's libcruft/odessa directory +# +# John W. Eaton +# jwe@bevo.che.wisc.edu +# University of Wisconsin-Madison +# Department of Chemical Engineering + +TOPDIR = ../.. + +srcdir = @srcdir@ +top_srcdir = @top_srcdir@ +VPATH = @srcdir@ + +EXTERNAL_DISTFILES = $(DISTFILES) + +include $(TOPDIR)/Makeconf + +include ../Makerules
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,1926 @@ +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- +C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA.. +C AN ORDINARY DIFFERENTIAL EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS +C SENSITIVITY ANALYSIS. +C +C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF +C LSODE.. LIVERMORE SOLVER FOR ORDINARY DIFFERENTIAL EQUATIONS. +C THIS VERSION IS IN DOUBLE PRECISION. +C +C ODESSA SOLVES FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS.. +C DY(I)/DP, FOR A SINGLE PARAMETER, OR, +C DY(I)/DP(J), FOR MULTIPLE PARAMETERS, +C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.. +C DY/DT = F(Y,T;P). +C----------------------------------------------------------------------- +C REFERENCES... +C +C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND +C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY +C DIFFERENTIAL EQUATIONS. SUBMITTED TO ACM TRANS. MATH. SOFTWARE, +C (1985). +C +C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY DIFFERENTIA +C EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS SENSITIVITY ANALYSIS. +C SUBMITTED TO ACM TRANS. MATH. SOFTWARE, (1985). +C +C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE +C ORDINARY DIFFERENTIAL EQUATION SOLVERS, ACM-SIGNUM NEWSLETTER, +C VOL. 15, NO. 4 (1980), PP. 10-11. +C----------------------------------------------------------------------- +C PROBLEM STATEMENT.. +C +C THE ODESSA MODIFICATION OF THE LSODE PACKAGE PROVIDES THE OPTION TO +C CALCULATE FIRST-ORDER SENSITIVITY COEFFICIENTS FOR A SYSTEM OF STIFF +C OR NON-STIFF EXPLICIT ORDINARY DIFFERENTIAL EQUATIONS OF THE GENERAL +C FORM : +C +C DY/DT = F(Y,T;P) (1) +C +C WHERE Y IS AN N-DIMENSIONAL DEPENDENT VARIABLE VECTOR, T IS THE +C INDEPENDENT INTEGRATION VARIABLE, AND P IS AN NPAR-DIMENSIONAL +C CONSTANT VECTOR. THE GOVERNING EQUATIONS FOR THE FIRST-ORDER +C SENSITIVITY COEFFICIENTS ARE GIVEN BY : +C +C S'(T) = J(T)*S(T) + DF/DP (2) +C +C WHERE +C +C S(T) = DY(T)/DP (= SENSITIVITY FUNCTIONS) +C S'(T) = D(DY(T)/DP)/DT +C J(T) = DF(Y,T;P)/DY(T) (= JACOBIAN MATRIX) +C AND DF/DP = DF(Y,T;P)/DP (= INHOMOGENEITY MATRIX) +C +C SOLUTION OF EQUATIONS (1) AND (2) PROCEEDS SIMULTANEOUSLY VIA AN +C EXTENSION OF THE LSODE PACKAGE AS DESCRIBED IN [1]. +C---------------------------------------------------------------------- +C ACKNOWLEDGEMENT : THE FOLLOWING ODESSA PACKAGE DOCUMENTATION IS A +C MODIFICATION OF THE LSODE DOCUMENTATION WHICH +C ACCOMPANIES THE LSODE PACKAGE CODE. +C---------------------------------------------------------------------- +C SUMMARY OF USAGE. +C +C COMMUNICATION BETWEEN THE USER AND THE ODESSA PACKAGE, FOR NORMAL +C SITUATIONS, IS SUMMARIZED HERE. THIS SUMMARY DESCRIBES ONLY A SUBSET +C OF THE FULL SET OF OPTIONS AVAILABLE. SEE THE FULL DESCRIPTION FOR +C DETAILS, INCLUDING OPTIONAL COMMUNICATION, NONSTANDARD OPTIONS, +C AND INSTRUCTIONS FOR SPECIAL SITUATIONS. SEE ALSO THE EXAMPLE +C PROBLEM (WITH PROGRAM AND OUTPUT) FOLLOWING THIS SUMMARY. +C +C A. FIRST PROVIDE A SUBROUTINE OF THE FORM.. +C SUBROUTINE F (N, T, Y, PAR, YDOT) +C DOUBLE PRECISION T, Y, PAR, YDOT +C DIMENSION Y(N), YDOT(N), PAR(NPAR) +C WHICH SUPPLIES THE VECTOR FUNCTION F BY LOADING YDOT(I) WITH F(I). +C N IS THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS IN THE +C ABOVE MODEL. NPAR IS THE NUMBER OF MODEL PARAMETERS FOR WHICH +C VECTOR SENSITIVITY FUNCTIONS ARE DESIRED. YOU ARE ALSO ENCOURAGED +C TO PROVIDE A SUBROUTINE OF THE FORM.. +C SUBROUTINE DF (N, T, Y, PAR, DFDP, JPAR) +C DOUBLE PRECISION T, Y, PAR, DFDP +C DIMENSION Y(N), PAR(NPAR), DFDP(N) +C GO TO (1,...,NPAR) JPAR +C 1 DFDP(1) = DF(1)/DP(1) +C . +C DFDP(I) = DF(I)/DP(1) +C . +C DFDP(N) = DF(N)/DP(1) +C RETURN +C 2 DFDP(1) = DF(1)/DP(2) +C . +C DFDP(I) = DF(I)/DP(2) +C . +C DFDP(N) = DF(N)/DP(2) +C RETURN +C . . +C . . +C RETURN +C NPAR DFDP(1) = DF(1)/DP(NPAR) +C . +C DFDP(I) = DF(I)/DP(NPAR) +C . +C DFDP(N) = DF(N)/DP(NPAR) +C RETURN +C END +C ONLY NONZERO ELEMENTS NEED BE LOADED. IF THIS IS NOT FEASIBLE, +C ODESSA WILL GENERATE THIS MATRIX INTERNALLY BY DIFFERENCE QUOTIENTS. +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 EIGENVALUE +C WHOSE REAL PART IS NEGATIVE AND LARGE IN MAGNITUDE, COMPARED TO THE +C RECIPROCAL OF THE T SPAN OF INTEREST. IF THE PROBLEM IS NONSTIFF, +C USE METH = 10. IF IT IS STIFF, USE METH = 20. THE USER IS REQUIRED +C TO INPUT THE METHOD FLAG MF = 10*METH + MITER. THERE ARE FOUR +C STANDARD CHOICES FOR MITER WHEN A SENSITIVITY ANALYSIS IS DESIRED, +C AND ODESSA REQUIRES THE JACOBIAN MATRIX IN SOME FORM. +C THIS MATRIX IS REGARDED EITHER AS FULL (MITER = 1 OR 2), +C OR BANDED (MITER = 4 OR 5). IN THE BANDED CASE, ODESSA REQUIRES TWO +C HALF-BANDWIDTH PARAMETERS ML AND MU. THESE ARE, RESPECTIVELY, THE +C WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, EXCLUDING THE MAIN +C DIAGONAL. THUS THE BAND CONSISTS OF THE LOCATIONS (I,J) WITH +C I-ML .LE. J .LE. I+MU, AND THE FULL BANDWIDTH IS ML+MU+1. +C +C C. YOU ARE ENCOURAGED TO SUPPLY THE JACOBIAN DIRECTLY (MF = 11, 14, +C 21, OR 24), BUT IF THIS IS NOT FEASIBLE, ODESSA WILL COMPUTE IT +C INTERNALLY BY DIFFERENCE QUOTIENTS (MF = 12, 15, 22, OR 25). IF YOU +C ARE SUPPLYING THE JACOBIAN, PROVIDE A SUBROUTINE OF THE FORM.. +C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD) +C DOUBLE PRECISION T, Y, PAR, PD +C DIMENSION Y(N), PD(NROWPD,N), PAR(NPAR) +C WHICH SUPPLIES DF/DY BY LOADING PD AS FOLLOWS.. +C FOR A FULL JACOBIAN (MF = 11, OR 21), LOAD PD(I,J) WITH DF(I)/DY(J), +C THE PARTIAL DERIVATIVE OF F(I) WITH RESPECT TO Y(J). (IGNORE THE +C ML AND MU ARGUMENTS IN THIS CASE.) +C FOR A BANDED JACOBIAN (MF = 14, OR 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 ROWS OF +C PD FROM THE TOP DOWN. +C IN EITHER CASE, ONLY NONZERO ELEMENTS NEED BE LOADED. +C +C D. WRITE A MAIN PROGRAM WHICH CALLS SUBROUTINE ODESSA ONCE FOR +C EACH POINT AT WHICH ANSWERS ARE DESIRED. THIS SHOULD ALSO PROVIDE +C FOR POSSIBLE USE OF LOGICAL UNIT 6 FOR OUTPUT OF ERROR MESSAGES BY +C ODESSA. ON THE FIRST CALL TO ODESSA, SUPPLY ARGUMENTS AS FOLLOWS.. +C F = NAME OF SUBROUTINE FOR RIGHT-HAND SIDE VECTOR F (MODEL). +C THIS NAME MUST BE DECLARED EXTERNAL IN CALLING PROGRAM. +C DF = NAME OF SUBROUTINE FOR INHOMOGENEITY MATRIX DF/DP. +C IF USED (IDF = 1), THIS NAME MUST BE DECLARED EXTERNAL IN +C CALLING PROGRAM. IF NOT USED (IDF = 0), PASS A DUMMY NAME. +C N = NUMBER OF FIRST ORDER ODE-S IN MODEL; LOAD INTO NEQ(1). +C NPAR = NUMBER OF MODEL PARAMETERS OF INTEREST; LOAD INTO NEQ(2). +C Y = AN (N) BY (NPAR+1) REAL ARRAY OF INITIAL VALUES.. +C Y(I,1) , I = 1,N , CONTAIN THE STATE, OR MODEL, DEPENDENT +C VARIABLES, +C Y(I,J) , J = 2,NPAR , CONTAIN THE DEPENDENT SENSITIVITY +C COEFFICIENTS. +C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING MODEL PARAMETERS +C OF INTEREST. +C T = THE INITIAL VALUE OF THE INDEPENDENT VARIABLE. +C TOUT = FIRST POINT WHERE OUTPUT IS DESIRED (.NE. T). +C ITOL = 1, 2, 3, OR 4 ACCORDING AS RTOL, ATOL (BELOW) ARE SCALARS +C OR ARRAYS. +C RTOL = RELATIVE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1) +C ARRAY). +C ATOL = ABSOLUTE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1) +C ARRAY). +C THE ESTIMATED LOCAL ERROR IN Y(I,J) WILL BE CONTROLLED SO AS +C TO BE ROUGHLY LESS (IN MAGNITUDE) THAN +C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL IF ITOL = 1, +C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL(I,J) IF ITOL = 2, +C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL IF ITOL = 3, OR +C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL(I,J) IF ITOL = 4. +C THUS THE LOCAL ERROR TEST PASSES IF, IN EACH COMPONENT, +C EITHER THE ABSOLUTE ERROR IS LESS THAN ATOL (OR ATOL(I,J)), +C OR THE RELATIVE ERROR IS LESS THAN RTOL (OR RTOL(I,J)). +C USE RTOL = 0.0 FOR PURE ABSOLUTE ERROR CONTROL, AND +C USE ATOL = 0.0 FOR PURE RELATIVE ERROR CONTROL. +C CAUTION.. ACTUAL (GLOBAL) ERRORS MAY EXCEED THESE LOCAL +C TOLERANCES, SO CHOOSE THEM CONSERVATIVELY. +C ITASK = 1 FOR NORMAL COMPUTATION OF OUTPUT VALUES OF Y AT T = TOUT. +C ISTATE = INTEGER FLAG (INPUT AND OUTPUT). SET ISTATE = 1. +C IOPT = 0, TO INDICATE NO OPTIONAL INPUTS FOR INTEGRATION; +C LOAD INTO IOPT(1). +C ISOPT = 1, TO INDICATE SENSITIVITY ANALYSIS, = 0, TO INDICATE +C NO SENSITIVITY ANALYSIS; LOAD INTO IOPT(2). +C IDF = 1, IF SUBROUTINE DF (ABOVE) IS SUPPLIED BY THE USER, +C = 0, OTHERWISE; LOAD INTO IOPT(3). +C RWORK = REAL WORK ARRAY OF LENGTH AT LEAST.. +C 22 + 16*N + N**2 FOR MF = 11 OR 12, +C 22 + 17*N + (2*ML + MU)*N FOR MF = 14 OR 15, +C 22 + 9*N + N**2 FOR MF = 21 OR 22, +C 22 + 10*N + (2*ML + MU)*N FOR MF = 24 OR 25, +C IF ISOPT = 0, OR.. +C 22 + 15*(NPAR+1)*N + N**2 + N FOR MF = 11 OR 12, +C 24 + 15*(NPAR+1)*N + (2*ML+MU+2)*N + N FOR MF = 14 OR 15, +C 22 + 8*(NPAR+1)*N + N**2 + N FOR MF = 21 OR 22, +C 24 + 8*(NPAR+1)*N + (2*ML+MU+2)*N + N FOR MF = 24 OR 25, +C IF ISOPT = 1. +C LRW = DECLARED LENGTH OF RWORK (IN USER-S DIMENSION STATEMENT). +C IWORK = INTEGER WORK ARRAY OF LENGTH AT LEAST.. +C 20 + N IF ISOPT = 0, +C 21 + N + NPAR IF ISOPT = 1. +C IF MITER = 4 OR 5, INPUT IN IWORK(1),IWORK(2) THE LOWER +C AND UPPER HALF-BANDWIDTHS ML,MU (EXCLUDING MAIN DIAGONAL). +C LIW = DECLARED LENGTH OF IWORK (IN USER-S DIMENSION STATEMENT). +C JAC = NAME OF SUBROUTINE FOR JACOBIAN MATRIX. +C IF USED, THIS NAME MUST BE DECLARED EXTERNAL IN CALLING +C PROGRAM. IF NOT USED, PASS A DUMMY NAME. +C MF = METHOD FLAG. STANDARD VALUES FOR ISOPT = 0 ARE.. +C 10 FOR NONSTIFF (ADAMS) METHOD, NO JACOBIAN USED. +C 21 FOR STIFF (BDF) METHOD, USER-SUPPLIED FULL JACOBIAN. +C 22 FOR STIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN. +C 24 FOR STIFF METHOD, USER-SUPPLIED BANDED JACOBIAN. +C 25 FOR STIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN. +C IF ISOPT = 1, MF = 10 IS ILLEGAL AND CAN BE REPLACED BY.. +C 11 FOR NONSTIFF METHOD, USER-SUPPLIED FULL JACOBIAN. +C 12 FOR NONSTIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN. +C 14 FOR NONSTIFF METHOD, USER-SUPPLIED BANDED JACOBIAN. +C 15 FOR NONSTIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN. +C NOTE THAT THE MAIN PROGRAM MUST DECLARE ARRAYS Y, RWORK, IWORK, AND +C POSSIBLY ATOL AND RTOL, AS WELL AS NEQ, IOPT, AND PAR IF ISOPT = 1. +C +C E. THE OUTPUT FROM THE FIRST CALL (OR ANY CALL) IS.. +C Y = ARRAY OF COMPUTED VALUES OF Y(T) VECTOR. +C T = CORRESPONDING VALUE OF INDEPENDENT VARIABLE (NORMALLY TOUT). +C ISTATE = 2 IF ODESSA WAS SUCCESSFUL, NEGATIVE OTHERWISE. +C -1 MEANS EXCESS WORK DONE ON THIS CALL (PERHAPS WRONG MF). +C -2 MEANS EXCESS ACCURACY REQUESTED (TOLERANCES TOO SMALL). +C -3 MEANS ILLEGAL INPUT DETECTED (SEE PRINTED MESSAGE). +C -4 MEANS REPEATED ERROR TEST FAILURES (CHECK ALL INPUTS). +C -5 MEANS REPEATED CONVERGENCE FAILURES (PERHAPS BAD JACOBIAN +C SUPPLIED OR WRONG CHOICE OF MF OR TOLERANCES). +C -6 MEANS ERROR WEIGHT BECAME ZERO DURING PROBLEM. (SOLUTION +C COMPONENT I,J VANISHED, AND ATOL OR ATOL(I,J) = 0.0) +C +C F. TO CONTINUE THE INTEGRATION AFTER A SUCCESSFUL RETURN, SIMPLY +C RESET TOUT AND CALL ODESSA AGAIN. NO OTHER PARAMETERS NEED BE RESET. +C---------------------------------------------------------------------- +C EXAMPLE PROBLEM. +C +C THE FOLLOWING IS A SIMPLE EXAMPLE PROBLEM, WITH THE CODING +C NEEDED FOR ITS SOLUTION BY ODESSA. THE PROBLEM IS FROM CHEMICAL +C KINETICS, AND CONSISTS OF THE FOLLOWING THREE RATE EQUATIONS.. +C DY1/DT = -PAR(1)*Y1 + PAR(2)*Y2*Y3 ; PAR(1) = .04, PAR(2) = 1.E4 +C DY2/DT = PAR(1)*Y1 - PAR(2)*Y2*Y3 - PAR(3)*Y2**2 ; PAR(3) = 3.E7 +C DY3/DT = PAR(3)*Y2**2 +C ON THE INTERVAL FROM T = 0.0 TO T = 4.E10, WITH INITIAL CONDITIONS +C Y1 = 1.0, Y2 = Y3 = 0, AND S(I,J) = 0, I = 1,3, J = 1,3. +C THE PROBLEM IS STIFF. +C +C THE FOLLOWING CODING SOLVES THIS PROBLEM WITH ODESSA, USING +C MF = 21 AND PRINTING RESULTS AT T = .4, 4., ..., 4.E10. +C IT USES ITOL = 4 AND ATOL MUCH SMALLER FOR Y2 THAN Y1 OR Y3, +C BECAUSE Y2 HAS MUCH SMALLER VALUES. LESS STRINGENT TOLERANCES +C ARE ASSIGNED FOR THE SENSITIVITIES TO ACHIEVE GREATER EFFICIENCY. +C AT THE END OF THE RUN, STATISTICAL QUANTITIES OF INTEREST ARE +C PRINTED (SEE OPTIONAL OUTPUTS IN THE FULL DESCRIPTION BELOW). +C +C DOUBLE PRECISION ATOL, RWORK, RTOL, T, TOUT, Y, PAR +C EXTERNAL FEX, JEX, DFEX +C DIMENSION Y(3,4), PAR(3), ATOL(3,4), RTOL(3,4), RWORK(130), +C 1 IWORK(27), NEQ(2), IOPT(3) +C N = 3 +C NPAR = 3 +C NEQ(1) = N +C NEQ(2) = NPAR +C NSV = NPAR+1 +C DO 10 I = 1,N +C DO 10 J = 1,NSV +C 10 Y(I,J) = 0.0D0 +C Y(1,1) = 1.0D0 +C PAR(1) = 0.04D0 +C PAR(2) = 1.0D4 +C PAR(3) = 3.0D7 +C T = 0.D0 +C TOUT = .4D0 +C ITOL = 4 +C ATOL(1,1) = 1.D-6 +C ATOL(2,1) = 1.D-10 +C ATOL(3,1) = 1.D-6 +C DO 20 I = 1,N +C RTOL(I,1) = 1.D-4 +C DO 15 J = 2,NSV +C RTOL(I,J) = 1.D-3 +C 15 ATOL(I,J) = 1.D2 * ATOL(I,1) +C 20 CONTINUE +C ITASK = 1 +C ISTATE = 1 +C IOPT(1) = 0 +C IOPT(2) = 1 +C IOPT(3) = 1 +C LRW = 130 +C LIW = 27 +C MF = 21 +C DO 60 IOUT = 1,12 +C CALL ODESSA(FEX,DFEX,NEQ,Y,PAR,T,TOUT,ITOL,RTOL,ATOL, +C 1 ITASK,ISTATE, IOPT,RWORK,LRW,IWORK,LIW,JEX,MF) +C WRITE(6,30)T,Y(1,1),Y(2,1),Y(3,1) +C 30 FORMAT(1X,7H AT T =,E12.4,6H Y =,3E14.6) +C DO 50 J = 2,NSV +C JPAR = J-1 +C WRITE(6,40)JPAR,Y(1,J),Y(2,J),Y(3,J) +C 40 FORMAT(20X,2HS(,I1,3H) =,3E14.6) +C 50 CONTINUE +C IF (ISTATE .LT. 0) GO TO 80 +C 60 TOUT = TOUT*10.D0 +C WRITE(6,70)IWORK(11),IWORK(12),IWORK(13),IWORK(19) +C 70 FORMAT(1X,/,12H NO. STEPS =,I4,11H NO. F-S =,I4,11H NO. J-S =, +C 1 I4,12H NO. DF-S =,I4) +C STOP +C 80 WRITE(6,90)ISTATE +C 90 FORMAT(///22H ERROR HALT.. ISTATE =,I3) +C STOP +C END +C +C SUBROUTINE FEX (NEQ, T, Y, PAR, YDOT) +C DOUBLE PRECISION T, Y, YDOT, PAR +C DIMENSION Y(3), YDOT(3), PAR(3) +C YDOT(1) = -PAR(1)*Y(1) + PAR(2)*Y(2)*Y(3) +C YDOT(3) = PAR(3)*Y(2)*Y(2) +C YDOT(2) = -YDOT(1) - YDOT(3) +C RETURN +C END +C +C SUBROUTINE JEX (NEQ, T, Y, PAR, ML, MU, PD, NRPD) +C DOUBLE PRECISION PD, T, Y, PAR +C DIMENSION Y(3), PD(NRPD,3), PAR(3) +C PD(1,1) = -PAR(1) +C PD(1,2) = PAR(2)*Y(3) +C PD(1,3) = PAR(2)*Y(2) +C PD(2,1) = PAR(1) +C PD(2,3) = -PD(1,3) +C PD(3,2) = 2.D0*PAR(3)*Y(2) +C PD(2,2) = -PD(1,2) - PD(3,2) +C RETURN +C END +C +C SUBROUTINE DFEX (NEQ, T, Y, PAR, DFDP, JPAR) +C DOUBLE PRECISION T, Y, PAR, DFDP +C DIMENSION Y(3), PAR(3), DFDP(3) +C GO TO (1,2,3), JPAR +C 1 DFDP(1) = -Y(1) +C DFDP(2) = Y(1) +C RETURN +C 2 DFDP(1) = Y(2)*Y(3) +C DFDP(2) = -Y(2)*Y(3) +C RETURN +C 3 DFDP(2) = -Y(2)*Y(2) +C DFDP(3) = Y(2)*Y(2) +C RETURN +C END +C +C THE OUTPUT OF THIS PROGRAM (ON A DATA GENERAL MV-8000 IN +C DOUBLE PRECISION IS AS FOLLOWS: +C +C AT T = .4000E+00 Y = .985173E+00 .338641E-04 .147930E-01 +C S(1) = -.355914E+00 .390261E-03 .355524E+00 +C S(2) = .955150E-07 -.213065E-09 -.953019E-07 +C S(3) = -.158466E-10 -.529012E-12 .163756E-10 +C AT T = .4000E+01 Y = .905516E+00 .224044E-04 .944615E-01 +C S(1) = -.187621E+01 .179197E-03 .187603E+01 +C S(2) = .296093E-05 -.583104E-09 -.296034E-05 +C S(3) = -.493267E-09 -.276246E-12 .493544E-09 +C AT T = .4000E+02 Y = .715848E+00 .918628E-05 .284143E+00 +C S(1) = -.424730E+01 .459360E-04 .424726E+01 +C S(2) = .137294E-04 -.235815E-09 -.137291E-04 +C S(3) = -.228818E-08 -.113803E-12 .228829E-08 +C AT T = .4000E+03 Y = .450526E+00 .322299E-05 .549471E+00 +C S(1) = -.595837E+01 .354310E-05 .595836E+01 +C S(2) = .227380E-04 -.226041E-10 -.227380E-04 +C S(3) = -.378971E-08 -.499501E-13 .378976E-08 +C AT T = .4000E+04 Y = .183185E+00 .894131E-06 .816814E+00 +C S(1) = -.475006E+01 -.599504E-05 .475007E+01 +C S(2) = .188089E-04 .231330E-10 -.188089E-04 +C S(3) = -.313478E-08 -.187575E-13 .313480E-08 +C AT T = .4000E+05 Y = .389733E-01 .162133E-06 .961027E+00 +C S(1) = -.157477E+01 -.276199E-05 .157477E+01 +C S(2) = .628668E-05 .110026E-10 -.628670E-05 +C S(3) = -.104776E-08 -.453588E-14 .104776E-08 +C AT T = .4000E+06 Y = .493609E-02 .198411E-07 .995064E+00 +C S(1) = -.236244E+00 -.458262E-06 .236244E+00 +C S(2) = .944669E-06 .183193E-11 -.944671E-06 +C S(3) = -.157441E-09 -.635990E-15 .157442E-09 +C AT T = .4000E+07 Y = .516087E-03 .206540E-08 .999484E+00 +C S(1) = -.256277E-01 -.509808E-07 .256278E-01 +C S(2) = .102506E-06 .203905E-12 -.102506E-06 +C S(3) = -.170825E-10 -.684002E-16 .170826E-10 +C AT T = .4000E+08 Y = .519314E-04 .207736E-09 .999948E+00 +C S(1) = -.259316E-02 -.518029E-08 .259316E-02 +C S(2) = .103726E-07 .207209E-13 -.103726E-07 +C S(3) = -.172845E-11 -.691450E-17 .172845E-11 +C AT T = .4000E+09 Y = .544710E-05 .217885E-10 .999995E+00 +C S(1) = -.271637E-03 -.541849E-09 .271638E-03 +C S(2) = .108655E-08 .216739E-14 -.108655E-08 +C S(3) = -.180902E-12 -.723615E-18 .180902E-12 +C AT T = .4000E+10 Y = .446748E-06 .178699E-11 .100000E+01 +C S(1) = -.322322E-04 -.842541E-10 .322323E-04 +C S(2) = .128929E-09 .337016E-15 -.128929E-09 +C S(3) = -.209715E-13 -.838859E-19 .209715E-13 +C AT T = .4000E+11 Y = -.363960E-07 -.145584E-12 .100000E+01 +C S(1) = -.164109E-06 -.429604E-11 .164113E-06 +C S(2) = .656436E-12 .171842E-16 -.656451E-12 +C S(3) = -.689361E-15 -.275745E-20 .689363E-15 +C +C NO. STEPS = 340 NO. F-S = 412 NO. J-S = 343 NO. DF-S =1023 +C---------------------------------------------------------------------- +C FULL DESCRIPTION OF USER INTERFACE TO ODESSA. +C +C THE USER INTERFACE TO ODESSA CONSISTS OF THE FOLLOWING PARTS. +C +C I. THE CALL SEQUENCE TO SUBROUTINE ODESSA, WHICH IS A DRIVER +C ROUTINE FOR THE SOLVER. THIS INCLUDES DESCRIPTIONS OF BOTH +C THE CALL SEQUENCE ARGUMENTS AND OF USER-SUPPLIED ROUTINES. +C FOLLOWING THESE DESCRIPTIONS IS A DESCRIPTION OF +C OPTIONAL INPUTS AVAILABLE THROUGH THE CALL SEQUENCE, AND THEN +C A DESCRIPTION OF OPTIONAL OUTPUTS (IN THE WORK ARRAYS). +C +C II. DESCRIPTIONS OF OTHER ROUTINES IN THE ODESSA 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 III. DESCRIPTIONS OF COMMON BLOCKS TO BE DECLARED IN OVERLAY +C OR SIMILAR ENVIRONMENTS, OR TO BE SAVED WHEN DOING AN INTERRUPT +C OF THE PROBLEM AND CONTINUED SOLUTION LATER. +C +C IV. DESCRIPTION OF TWO SUBROUTINES IN THE ODESSA 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 V. GENERAL REMARKS WHICH HIGHLIGHT DIFFERENCES BETWEEN THE LSODE +C PACKAGE AND THE ODESSA PACKAGE. +C---------------------------------------------------------------------- +C PART I. CALL SEQUENCE. +C +C THE CALL SEQUENCE PARAMETERS USED FOR INPUT ONLY ARE.. +C F, DF, NEQ, PAR, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, +C JAC, MF, +C AND THOSE USED FOR BOTH INPUT AND OUTPUT ARE +C Y, T, ISTATE. +C THE WORK ARRAYS RWORK AND IWORK ARE ALSO USED FOR CONDITIONAL AND +C OPTIONAL INPUTS AND OPTIONAL OUTPUTS. (THE TERM OUTPUT HERE REFERS +C TO THE RETURN FROM SUBROUTINE ODESSA TO THE USER-S CALLING 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 +C ODE MODEL. THIS SYSTEM MUST BE PUT IN THE FIRST-ORDER +C FORM DY/DT = F(Y,T;P), WHERE F IS A VECTOR-VALUED FUNCTION +C OF THE SCALAR T AND VECTORS Y, AND PAR. SUBROUTINE F IS TO +C COMPUTE THE FUNCTION F. IT IS TO HAVE THE FORM.. +C SUBROUTINE F (NEQ, T, Y, PAR, YDOT) +C DOUBLE PRECISION T, Y, PAR, YDOT +C DIMENSION Y(1), PAR(1), YDOT(1) +C WHERE NEQ, T, Y, AND PAR ARE INPUT, AND YDOT = F(Y,T;P) +C IS OUTPUT. Y AND YDOT ARE ARRAYS OF LENGTH N (= NEQ(1)). +C (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY +C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.) +C F SHOULD NOT ALTER ARRAY Y, OR PAR(1),...,PAR(NPAR). +C F MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. +C +C SUBROUTINE F MAY ACCESS USER-DEFINED QUANTITIES IN +C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY +C (DIMENSIONED IN F) AND PAR HAS LENGTH EXCEEDING NPAR. +C SEE THE DESCRIPTIONS OF NEQ AND PAR BELOW. +C +C DF = THE NAME OF THE USER-SUPPLIED ROUTINE (IDF = 1) TO COMPUTE +C THE INHOMOGENEITY MATRIX, DF/DP, AS A FUNCTION OF THE SCALAR +C T, AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM +C SUBROUTINE DF (NEQ, T, Y, PAR, DFDP, JPAR) +C DOUBLE PRECISION T, Y, PAR, DFDP +C DIMENSION Y(1), PAR(1), DFDP(1) +C GO TO (1,2,...,NPAR) JPAR +C 1 DFDP(1) = DF(1)/DP(1) +C . +C DFDP(I) = DF(I)/DP(1) +C . +C DFDP(N) = DF(N)/DP(1) +C RETURN +C 2 DFDP(1) = DF(1)/DP(2) +C . +C DFDP(I) = DF(I)/DP(2) +C . +C DFDP(N) = DF(N)/DP(2) +C . +C RETURN +C . . +C . . +C NPAR DFDP(1) = DF(1)/DP(NPAR) +C . +C DFDP(I) = DF(I)/DP(NPAR) +C . +C DFDP(N) = DF(N)/DP(NPAR) +C RETURN +C END +C WHERE NEQ, T, Y, PAR, AND JPAR ARE INPUT AND THE VECTOR +C DFDP(*,JPAR) IS TO BE LOADED WITH THE PARTIAL DERIVATIVES +C DF(Y,T;PAR)/DP(JPAR) ON OUTPUT. ONLY NONZERO ELEMENTS NEED +C BE LOADED. T, Y, AND PAR HAVE THE SAME MEANING AS IN +C SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY +C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE). +C +C DFDP(*,JPAR) IS PRESET TO ZERO BY THE SOLVER, SO THAT ONLY +C THE NONZERO ELEMENTS NEED BE LOADED BY DF. SUBROUTINE DF +C MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM IF USED. +C IF IDF = 0 (OR ISOPT = 0), A DUMMY ARGUMENT CAN BE USED. +C +C SUBROUTINE DF MAY ACCESS USER-DEFINED QUANTITIES IN +C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY +C (DIMENSIONED IN DF) AND PAR HAS A LENGTH EXCEEDING NPAR. +C SEE THE DESCRIPTIONS OF NEQ AND PAR (BELOW). +C +C NEQ = THE SIZE OF THE ODE SYSTEM (NUMBER OF FIRST ORDER ORDINARY +C DIFFERENTIAL EQUATIONS (N) IN THE MODEL). USED ONLY FOR +C INPUT. NEQ MAY NOT BE CHANGED DURING THE PROBLEM. +C +C FOR ISOPT = 0, NEQ IS NORMALLY A SCALAR. HOWEVER, NEQ MAY +C BE AN ARRAY, WITH NEQ(1) SET TO THE SYSTEM SIZE (N), IN WHICH +C CASE THE ODESSA PACKAGE ACCESSES ONLY NEQ(1). HOWEVER, +C THIS PARAMETER IS PASSED AS THE NEQ ARGUMENT IN ALL CALLS +C TO F, DF, AND JAC. HENCE, IF IT IS AN ARRAY, LOCATIONS +C NEQ(2),... MAY BE USED TO STORE OTHER INTEGER DATA AND PASS +C IT TO F, DF, AND/OR JAC. FOR ISOPT = 1, NPAR MUST BE LOADED +C INTO NEQ(2), AND IS NOT ALLOWED TO CHANGE DURING THE PROBLEM. +C IN THESE CASES, SUBROUTINES F, DF, AND/OR JAC MUST INCLUDE +C NEQ IN A DIMENSION STATEMENT. +C +C Y = A REAL ARRAY FOR THE VECTOR OF DEPENDENT VARIABLES, OF +C DIMENSION (N) BY (NPAR+1). USED FOR BOTH INPUT AND +C OUTPUT ON THE FIRST CALL (ISTATE = 1), AND ONLY FOR +C OUTPUT ON OTHER CALLS. ON THE FIRST CALL, Y MUST CONTAIN +C THE VECTORS OF INITIAL VALUES. ON OUTPUT, Y CONTAINS THE +C COMPUTED SOLUTION VECTORS, EVALUATED AT T. +C +C PAR = A REAL ARRAY FOR THE VECTOR OF CONSTANT MODEL PARAMETERS +C OF INTEREST IN THE SENSITIVITY ANALYSIS, OF LENGTH NPAR +C OR MORE. PAR IS PASSED AS AN ARGUMENT IN ALL CALLS TO F, +C DF, AND JAC. HENCE LOCATIONS PAR(NPAR+1),... MAY BE USED +C TO STORE OTHER REAL DATA AND PASS IT TO F, DF, AND/OR JAC. +C LOCATIONS PAR(1),...,PAR(NPAR) ARE USED AS INPUT ONLY, +C AND MUST NOT BE CHANGED DURING THE PROBLEM. +C +C T = THE INDEPENDENT VARIABLE. ON INPUT, T IS USED ONLY ON THE +C 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 TOUT). +C ON AN ERROR RETURN, T IS THE FARTHEST POINT REACHED. +C +C TOUT = THE NEXT VALUE OF T AT WHICH A COMPUTED SOLUTION IS DESIRED. +C USED ONLY FOR INPUT. +C +C WHEN STARTING THE PROBLEM (ISTATE = 1), TOUT MAY BE EQUAL +C TO T FOR ONE CALL, THEN SHOULD .NE. T FOR THE NEXT CALL. +C FOR THE INITIAL T, AN INPUT VALUE OF TOUT .NE. T IS USED +C IN ORDER TO DETERMINE THE DIRECTION OF THE INTEGRATION +C (I.E. THE ALGEBRAIC SIGN OF THE STEP SIZES) AND THE ROUGH +C SCALE OF THE PROBLEM. INTEGRATION IN EITHER DIRECTION +C (FORWARD OR BACKWARD IN T) IS PERMITTED. +C +C IF ITASK = 2 OR 5 (ONE-STEP MODES), TOUT IS IGNORED AFTER +C THE FIRST CALL (I.E. THE FIRST CALL WITH TOUT .NE. T). +C 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 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 SPACE (N) BY (NPAR+1). SEE DESCRIPTION BELOW +C UNDER ATOL. INPUT ONLY. +C +C ATOL = AN ABSOLUTE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR +C AN ARRAY OF SPACE (N) BY (NPAR+1). INPUT ONLY. +C +C THE INPUT PARAMETERS ITOL, RTOL, AND ATOL DETERMINE +C THE ERROR CONTROL PERFORMED BY THE SOLVER. THE SOLVER WILL +C CONTROL THE VECTOR E = (E(I,J)) OF ESTIMATED LOCAL ERRORS +C IN Y, ACCORDING TO AN INEQUALITY OF THE FORM +C RMS-NORM OF ( E(I,J)/EWT(I,J) ) .LE. 1, +C WHERE EWT(I,J) = RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J), +C AND THE RMS-NORM (ROOT-MEAN-SQUARE NORM) HERE IS +C RMS-NORM(V) = SQRT ( (1/N) * SUM (V(I,J)**2) ); I =1,...,N. +C HERE EWT = (EWT(I,J)) IS A VECTOR OF WEIGHTS WHICH MUST +C ALWAYS BE POSITIVE, AND THE VALUES OF RTOL AND ATOL SHOULD +C ALL BE NON-NEGATIVE. THE FOLLOWING TABLE GIVES THE TYPES +C (SCALAR/ARRAY) OF RTOL AND ATOL, AND THE CORRESPONDING FORM +C OF EWT(I,J). +C +C ITOL RTOL ATOL EWT(I,J) +C 1 SCALAR SCALAR RTOL*ABS(Y(I,J)) + ATOL +C 2 SCALAR ARRAY RTOL*ABS(Y(I,J)) + ATOL(I,J) +C 3 ARRAY SCALAR RTOL(I,J)*ABS(Y(I,J)) + ATOL +C 4 ARRAY ARRAY RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J) +C +C WHEN EITHER OF THESE PARAMETERS IS A SCALAR, IT NEED NOT +C BE DIMENSIONED IN THE USER-S CALLING PROGRAM. +C +C THE TOTAL NUMBER OF ERROR TEST FAILURES DUE TO THE SENSITIVITY +C ANALYSIS, AND WHICH REQUIRE AN INTEGRATION STEP TO BE +C REPEATED, ARE ACCUMULATED IN THE LAST NPAR+1 LOCATIONS OF THE +C INTEGER WORK ARRAY IWORK (SEE OPTIONAL OUTPUTS BELOW). +C THIS INFORMATION MAY BE OF VALUE IN DETERMINING APPROPRIATE +C ERROR TOLERANCES TO BE APPLIED TO THE SENSITIVITY FUNCTIONS. +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 IV 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 SCALED +C DOWN UNIFORMLY. +C +C ITASK = AN INDEX SPECIFYING THE TASK TO BE PERFORMED. +C INPUT ONLY. ITASK HAS THE FOLLOWING VALUES AND MEANINGS. +C 1 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT +C T = TOUT (BY OVERSHOOTING AND INTERPOLATING). +C 2 MEANS TAKE ONE STEP ONLY AND RETURN. +C 3 MEANS STOP AT THE FIRST INTERNAL MESH POINT AT OR +C BEYOND T = TOUT AND RETURN. +C 4 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT +C T = TOUT BUT WITHOUT OVERSHOOTING T = TCRIT. +C TCRIT MUST BE INPUT AS RWORK(1). TCRIT MAY BE EQUAL TO +C OR 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 MEANS 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 TCRIT, +C IN WHICH CASE ANSWERS AT T = TOUT ARE RETURNED FIRST). +C +C ISTATE = AN INDEX USED FOR INPUT AND OUTPUT TO SPECIFY THE +C THE STATE OF THE CALCULATION. +C +C ON INPUT, THE VALUES OF ISTATE ARE AS FOLLOWS. +C 1 MEANS THIS IS THE FIRST CALL FOR THE PROBLEM +C (INITIALIZATIONS WILL BE DONE). SEE NOTE BELOW. +C 2 MEANS THIS IS NOT THE FIRST CALL, AND THE CALCULATION +C IS TO CONTINUE NORMALLY, WITH NO CHANGE IN ANY INPUT +C PARAMETERS EXCEPT POSSIBLY TOUT AND ITASK. +C (IF ITOL, RTOL, AND/OR ATOL ARE CHANGED BETWEEN CALLS +C WITH ISTATE = 2, THE NEW VALUES WILL BE USED BUT NOT +C TESTED FOR LEGALITY.) +C 3 MEANS THIS IS NOT THE FIRST CALL, AND THE +C CALCULATION IS TO CONTINUE NORMALLY, BUT WITH +C A CHANGE IN INPUT PARAMETERS OTHER THAN +C TOUT AND ITASK. CHANGES ARE ALLOWED IN +C ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU, +C AND ANY OF THE OPTIONAL INPUTS EXCEPT H0. +C (SEE IWORK DESCRIPTION FOR ML AND MU.) +C NOTE.. A PRELIMINARY CALL WITH TOUT = T IS NOT COUNTED +C AS 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.) +C THUS THE FIRST CALL FOR WHICH TOUT .NE. T REQUIRES +C ISTATE = 1 ON INPUT. +C +C ON OUTPUT, ISTATE HAS THE FOLLOWING VALUES AND MEANINGS. +C 1 MEANS NOTHING WAS DONE, AS TOUT WAS EQUAL TO T WITH +C ISTATE = 1 ON INPUT. (HOWEVER, AN INTERNAL COUNTER WAS +C SET TO DETECT AND PREVENT REPEATED CALLS OF THIS TYPE.) +C 2 MEANS THE INTEGRATION WAS PERFORMED SUCCESSFULLY. +C -1 MEANS AN EXCESSIVE AMOUNT OF WORK (MORE THAN MXSTEP +C STEPS) 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 .GT. 1 AND CALL AGAIN +C (THE EXCESS WORK STEP COUNTER WILL BE RESET TO 0). +C IN ADDITION, THE USER MAY INCREASE MXSTEP TO AVOID +C THIS ERROR RETURN (SEE BELOW ON OPTIONAL INPUTS). +C -2 MEANS TOO MUCH ACCURACY WAS REQUESTED FOR THE PRECISION +C OF 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 TOLERANCE +C PARAMETERS MUST BE RESET, AND ISTATE MUST BE SET +C TO 3. THE OPTIONAL OUTPUT TOLSF MAY BE USED FOR THIS +C PURPOSE. (NOTE.. IF THIS CONDITION IS DETECTED BEFORE +C TAKING ANY STEPS, THEN AN ILLEGAL INPUT RETURN +C (ISTATE = -3) OCCURS INSTEAD.) +C -3 MEANS 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 CALLS +C TO THE SOLVER WITH ILLEGAL INPUT, IT WILL CAUSE +C THE RUN TO STOP. +C -4 MEANS THERE WERE REPEATED ERROR TEST FAILURES ON +C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED +C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T. +C THE PROBLEM MAY HAVE A SINGULARITY, OR THE INPUT +C MAY BE INAPPROPRIATE. +C -5 MEANS THERE WERE REPEATED CONVERGENCE TEST FAILURES ON +C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED +C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T. +C THIS MAY BE CAUSED BY AN INACCURATE JACOBIAN MATRIX, +C IF ONE IS BEING USED. +C -6 MEANS EWT(I,J) BECAME ZERO FOR SOME I,J DURING THE +C INTEGRATION. PURE RELATIVE ERROR CONTROL (ATOL(I,J)=0.0) +C WAS REQUESTED ON A VARIABLE WHICH HAS NOW VANISHED. +C THE INTEGRATION WAS SUCCESSFUL AS FAR AS T. +C +C NOTE.. SINCE THE NORMAL OUTPUT VALUE OF ISTATE IS 2, +C IT DOES NOT NEED TO BE RESET FOR NORMAL CONTINUATION. +C ALSO, SINCE A NEGATIVE INPUT VALUE OF ISTATE WILL BE +C REGARDED AS ILLEGAL, A NEGATIVE OUTPUT VALUE REQUIRES THE +C USER TO CHANGE IT, AND POSSIBLY OTHER INPUTS, BEFORE +C CALLING THE SOLVER AGAIN. +C +C IOPT = AN INTEGER ARRAY FLAG TO SPECIFY WHETHER OR NOT ANY OPTIONAL +C INPUTS ARE BEING USED ON THIS CALL. INPUT ONLY. +C THE OPTIONAL INPUTS ARE LISTED SEPARATELY BELOW. +C IOPT(1) = 0 MEANS NO OPTIONAL INPUTS FOR THE SOLVER WILL BE +C USED. DEFAULT VALUES WILL BE USED IN ALL CASES. +C = 1 MEANS ONE OR MORE OPTIONAL INPUTS FOR THE +C SOLVER ARE BEING USED. +C NOTE : IOPT(1) IS INDEPENDENT OF ISOPT AND IDF. +C IOPT(2) = 0 MEANS NO SENSITIVITY ANALYSIS WILL BE PERFORMED. +C = 1 MEANS A SENSITIVITY ANALYSIS WILL BE PERFORMED. +C NOTE : IOPT(2) IS RENAMED TO ISOPT IN ODESSA. +C = 0 MEANS DF/DP WILL BE CALCULATED BY FINITE +C DIFFERENCE WITHIN ODESSA. +C IOPT(3) = 1 MEANS DF/DP WILL BE CALCULATED BY A USER-SUPPLIED +C ROUTINE. +C NOTE : IOPT(3) IS RENAMED TO IDF IN ODESSA. +C IF IDF = 1, THE USER MUST SUPPLY A +C SUBROUTINE DF (THE NAME IS ARBITRARY) AS +C DESCRIBED BELOW UNDER DF. FOR IDF = 0, +C A DUMMY ARGUMENT CAN BE USED. +C +C RWORK = A REAL WORKING ARRAY (DOUBLE PRECISION). +C FOR ISOPT = 0, THE LENGTH OF RWORK MUST BE AT LEAST.. +C 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM +C FOR ISOPT = 1, THE LENGTH OF RWORK MUST BE AT LEAST.. +C 20 + NYH*(MAXORD + 1) + 2*NYH + LWM + N +C WHERE.. +C NYH = THE TOTAL NUMBER OF DEPENDENT VARIABLES; +C (= N IF ISOPT = 0, AND N*(NPAR+1) IF ISOPT = 1). +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 = N**2 + 2 IF MITER IS 1 OR 2, +C LWM = N + 2 IF MITER = 3, AND +C LWM = (2*ML+MU+1)*N + 2 IF MITER IS 4 OR 5. +C (SEE THE MF DESCRIPTION FOR METH AND MITER.) +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 = CRITICAL VALUE OF T WHICH THE SOLVER +C IS NOT TO OVERSHOOT. REQUIRED IF ITASK IS +C 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. THE LENGTH MUST BE AT LEAST.. +C 20 IF MITER = 0 OR 3 (MF = 10, 13, 20, 23), OR +C 20 + N OTHERWISE (MF = 11, 12, 14, 15, 21, 22, 24, 25). +C FOR ISOPT = 0, OR.. +C 21 + N + NPAR +C FOR ISOPT = 1. +C THE FIRST FEW WORDS OF IWORK ARE USED FOR CONDITIONAL AND +C OPTIONAL INPUTS AND OPTIONAL OUTPUTS. +C +C THE FOLLOWING 2 WORDS IN IWORK ARE CONDITIONAL INPUTS.. +C IWORK(1) = ML THESE ARE THE LOWER AND UPPER +C IWORK(2) = MU HALF-BANDWIDTHS, RESPECTIVELY, OF THE +C BANDED JACOBIAN, EXCLUDING THE MAIN DIAGONAL. +C THE BAND IS DEFINED BY THE MATRIX LOCATIONS +C (I,J) WITH I-ML .LE. J .LE. I+MU. ML AND MU +C MUST SATISFY 0 .LE. ML,MU .LE. NEQ-1. +C THESE ARE REQUIRED IF MITER IS 4 OR 5, AND +C IGNORED OTHERWISE. ML AND MU MAY IN FACT BE +C THE BAND PARAMETERS FOR A MATRIX TO WHICH +C DF/DY IS ONLY 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 ODESSA +C FOR THE SAME PROBLEM, EXCEPT POSSIBLY FOR THE CONDITIONAL AND +C OPTIONAL INPUTS, AND EXCEPT FOR THE LAST 2*NYH + N WORDS OF RWORK. +C THE LATTER SPACE IS USED FOR INTERNAL SCRATCH SPACE, AND SO IS +C AVAILABLE FOR USE BY THE USER OUTSIDE ODESSA BETWEEN CALLS, IF +C DESIRED (BUT NOT FOR USE BY F, DF, 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 VECTORS Y, AND PAR. IT IS TO HAVE THE FORM +C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD) +C DOUBLE PRECISION T, Y, PAR, PD +C DIMENSION Y(1), PAR(1), PD(NROWPD,1) +C WHERE NEQ, T, Y, PAR, ML, MU, AND NROWPD ARE INPUT AND THE +C ARRAY PD IS TO BE LOADED WITH PARTIAL DERIVATIVES (ELEMENTS +C OF THE JACOBIAN MATRIX) ON OUTPUT. PD MUST BE GIVEN A FIRST +C DIMENSION OF NROWPD. T, Y, AND PAR HAVE THE SAME MEANING AS +C IN SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A +C DUMMY DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.) +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 IN THE BAND MATRIX CASE (MITER = 4), THE ELEMENTS +C WITHIN THE BAND ARE TO BE LOADED INTO PD IN COLUMNWISE +C MANNER, WITH DIAGONAL LINES OF DF/DY LOADED INTO THE ROWS +C OF PD. THUS DF(I)/DY(J) IS TO BE LOADED INTO PD(I-J+MU+1,J). +C ML 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 ODESSA. +C PD IS PRESET TO ZERO BY THE SOLVER, SO THAT ONLY THE +C NONZERO ELEMENTS NEED BE LOADED BY JAC. EACH CALL TO JAC IS +C PRECEDED BY A CALL TO F WITH THE SAME ARGUMENTS NEQ, T, Y, +C AND PAR. THUS TO GAIN SOME EFFICIENCY, INTERMEDIATE +C QUANTITIES SHARED BY BOTH CALCULATIONS MAY BE SAVED IN A +C USER COMMON BLOCK BY F AND NOT RECOMPUTED BY JAC, IF +C DESIRED. ALSO, JAC MAY ALTER THE Y ARRAY, IF DESIRED. +C JAC MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM. +C SUBROUTINE JAC MAY ACCESS USER-DEFINED QUANTITIES IN +C NEQ(2),... AND PAR(NPAR+1),.... SEE THE DESCRIPTIONS OF +C NEQ (ABOVE) AND PAR (BELOW). +C +C MF = THE METHOD FLAG. USED ONLY FOR INPUT. THE LEGAL VALUES OF +C MF ARE 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, AND 25. +C MF HAS DECIMAL DIGITS METH AND MITER.. MF = 10*METH + MITER. +C METH INDICATES THE BASIC LINEAR MULTISTEP METHOD.. +C METH = 1 MEANS THE IMPLICIT ADAMS METHOD. +* +C METH = 2 MEANS THE METHOD BASED ON BACKWARD +C DIFFERENTIATION FORMULAS (BDF-S). +C MITER INDICATES THE CORRECTOR ITERATION METHOD.. +C MITER = 0 MEANS FUNCTIONAL ITERATION (NO JACOBIAN MATRIX +C IS INVOLVED). +C MITER = 1 MEANS CHORD ITERATION WITH A USER-SUPPLIED +C FULL (NEQ BY NEQ) JACOBIAN. +C MITER = 2 MEANS CHORD ITERATION WITH AN INTERNALLY +C GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN +C (USING NEQ EXTRA CALLS TO F PER DF/DY VALUE). +C MITER = 3 MEANS CHORD ITERATION WITH AN INTERNALLY +C GENERATED DIAGONAL JACOBIAN APPROXIMATION. +C (USING 1 EXTRA CALL TO F PER DF/DY EVALUATION). +C MITER = 4 MEANS CHORD ITERATION WITH A USER-SUPPLIED +C BANDED JACOBIAN. +C MITER = 5 MEANS CHORD ITERATION WITH AN INTERNALLY +C GENERATED BANDED JACOBIAN (USING ML+MU+1 EXTRA +C CALLS TO F PER DF/DY EVALUATION). +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 IF A SENSITIVITY ANLYSIS IS DESIRED (ISOPT = 1), MITER = 0 +C AND 3 ARE DISALLOWED. IN THESE CASES, THE USER IS RECOMMENDED +C TO SUPPLY AN ANALYTICAL JACOBIAN (MITER = 1 OR 4) AND AN +C ANALYTICAL INHOMOGENEITY MATRIX (IDF = 1). +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 II.) 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) = 1, AND IN THAT +C CASE ALL OF THESE INPUTS ARE EXAMINED. A VALUE OF ZERO FOR ANY +C OF THESE OPTIONAL INPUTS WILL CAUSE THE DEFAULT VALUE TO BE USED. +C THUS TO USE A SUBSET OF THE OPTIONAL INPUTS, SIMPLY PRELOAD +C LOCATIONS 5 TO 10 IN RWORK AND IWORK TO 0.0 AND 0 RESPECTIVELY, AND +C THEN SET THOSE OF INTEREST TO NONZERO VALUES. +C +C NAME LOCATION MEANING AND DEFAULT VALUE +C +C H0 RWORK(5) THE STEP SIZE TO BE ATTEMPTED ON THE FIRST STEP. +C THE DEFAULT VALUE IS DETERMINED BY THE SOLVER. +C +C HMAX RWORK(6) THE MAXIMUM ABSOLUTE STEP SIZE ALLOWED. +C THE DEFAULT VALUE IS INFINITE. +C +C HMIN RWORK(7) THE MINIMUM ABSOLUTE STEP SIZE ALLOWED. +C THE DEFAULT VALUE IS 0. (THIS LOWER BOUND IS NOT +C ENFORCED ON THE FINAL STEP BEFORE REACHING TCRIT +C WHEN ITASK = 4 OR 5.) +C +C MAXORD IWORK(5) THE MAXIMUM ORDER TO BE ALLOWED. THE DEFAULT +C VALUE IS 12 IF METH = 1, AND 5 IF METH = 2. +C IF MAXORD EXCEEDS THE DEFAULT VALUE, IT WILL +C BE REDUCED TO THE DEFAULT VALUE. +C IF MAXORD IS CHANGED DURING THE PROBLEM, IT MAY +C CAUSE THE CURRENT ORDER TO BE REDUCED. +C +C MXSTEP IWORK(6) MAXIMUM NUMBER OF (INTERNALLY DEFINED) STEPS +C ALLOWED DURING ONE CALL TO THE SOLVER. +C THE DEFAULT VALUE IS 500. +C +C MXHNIL IWORK(7) MAXIMUM NUMBER OF MESSAGES PRINTED (PER PROBLEM) +C WARNING THAT T + H = T ON A STEP (H = STEP SIZE). +C THIS MUST BE POSITIVE TO RESULT IN A NON-DEFAULT +C VALUE. THE DEFAULT VALUE IS 10. +C---------------------------------------------------------------------- +C OPTIONAL OUTPUTS. +C +C AS OPTIONAL ADDITIONAL OUTPUT FROM ODESSA, THE VARIABLES LISTED +C BELOW ARE QUANTITIES RELATED TO THE PERFORMANCE OF ODESSA +C WHICH ARE AVAILABLE TO THE USER. THESE ARE COMMUNICATED BY WAY OF +C THE WORK ARRAYS, BUT ALSO HAVE INTERNAL MNEMONIC NAMES AS SHOWN. +C EXCEPT WHERE STATED OTHERWISE, ALL OF THESE OUTPUTS ARE DEFINED +C ON ANY SUCCESSFUL RETURN FROM ODESSA, AND ON ANY RETURN WITH +C ISTATE = -1, -2, -4, -5, OR -6. ON AN ILLEGAL INPUT RETURN +C (ISTATE = -3), THEY WILL BE UNCHANGED FROM THEIR EXISTING VALUES +C (IF ANY), EXCEPT POSSIBLY FOR TOLSF, LENRW, AND LENIW. +C ON ANY ERROR RETURN, OUTPUTS RELEVANT TO THE ERROR WILL BE DEFINED, +C AS NOTED BELOW. +C +C NAME LOCATION MEANING +C +C HU RWORK(11) THE STEP SIZE IN T LAST USED (SUCCESSFULLY). +C +C HCUR RWORK(12) THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. +C +C TCUR RWORK(13) THE CURRENT VALUE OF THE INDEPENDENT VARIABLE +C WHICH THE SOLVER HAS ACTUALLY REACHED, I.E. THE +C CURRENT INTERNAL MESH POINT IN T. ON OUTPUT, TCUR +C WILL ALWAYS BE AT LEAST AS FAR AS THE ARGUMENT +C T, BUT MAY BE FARTHER (IF INTERPOLATION WAS DONE). +C +C TOLSF RWORK(14) A TOLERANCE SCALE FACTOR, GREATER THAN 1.0, +C COMPUTED WHEN A REQUEST FOR TOO MUCH ACCURACY WAS +C DETECTED (ISTATE = -3 IF DETECTED AT THE START OF +C THE PROBLEM, ISTATE = -2 OTHERWISE). IF ITOL IS +C LEFT UNALTERED BUT RTOL AND ATOL ARE UNIFORMLY +C SCALED UP BY A FACTOR OF TOLSF FOR THE NEXT CALL, +C THEN THE SOLVER IS DEEMED LIKELY TO SUCCEED. +C (THE USER MAY ALSO IGNORE TOLSF AND ALTER THE +C TOLERANCE PARAMETERS IN ANY OTHER WAY APPROPRIATE.) +C +C NST IWORK(11) THE NUMBER OF STEPS TAKEN FOR THE PROBLEM SO FAR. +C +C NFE IWORK(12) THE NUMBER OF F EVALUATIONS FOR THE PROBLEM SO FAR. +C +C NJE IWORK(13) THE NUMBER OF JACOBIAN EVALUATIONS (AND OF MATRIX +C LU DECOMPOSITIONS IF ISOPT = 0) FOR THE PROBLEM SO +C FAR. IF ISOPT = 1, THE NUMBER OF LU DECOMPOSITIONS +C IS EQUAL TO NJE - NSPE (SEE BELOW). +C +C NQU IWORK(14) THE METHOD ORDER LAST USED (SUCCESSFULLY). +C +C NQCUR IWORK(15) THE ORDER TO BE ATTEMPTED ON THE NEXT STEP. +C +C IMXER IWORK(16) THE INDEX OF THE COMPONENT OF LARGEST MAGNITUDE IN +C THE WEIGHTED LOCAL ERROR VECTOR (E(I,J)/EWT(I,J)), +C ON AN ERROR RETURN WITH ISTATE = -4 OR -5. +C +C LENRW IWORK(17) THE LENGTH OF RWORK ACTUALLY REQUIRED. +C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL +C INPUT RETURN FOR INSUFFICIENT STORAGE. +C +C LENIW IWORK(18) THE LENGTH OF IWORK ACTUALLY REQUIRED. +C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL +C INPUT RETURN FOR INSUFFICIENT STORAGE. +C +C NDFE IWORK(19) THE NUMBER OF DF/DP (VECTOR) EVALUATIONS. +C +C NSPE IWORK(20) THE NUMBER OF CALLS TO SUBROUTINE ODESSA_SPRIME. EACH CALL +C TO ODESSA_SPRIME REQUIRES A JACOBIAN EVALUATION, BUT NOT +C AN LU DECOMPOSITION. +C +C THE FOLLOWING ARRAYS ARE SEGMENTS OF THE RWORK AND IWORK ARRAYS +C WHICH MAY ALSO BE OF INTEREST TO THE USER AS OPTIONAL OUTPUTS. +C FOR EACH ARRAY, THE TABLE BELOW GIVES ITS INTERNAL NAME, ITS BASE +C ADDRESS IN RWORK OR IWORK, AND ITS DESCRIPTION. +C +C NAME BASE ADDRESS DESCRIPTION +C +C YH 21 IN RWORK THE NORDSIECK HISTORY ARRAY, OF SIZE NYH BY +C (NQCUR + 1). FOR J = 0,1,...,NQCUR, COLUMN J+1 +C OF YH CONTAINS HCUR**J/FACTORIAL(J) TIMES +C THE J-TH DERIVATIVE OF THE INTERPOLATING +C POLYNOMIAL CURRENTLY REPRESENTING THE SOLUTION, +C EVALUATED AT T = TCUR. +C +C ACOR LENRW-NYH+1 ARRAY OF SIZE NYH USED FOR THE ACCUMULATED +C IN RWORK CORRECTIONS ON EACH STEP, SCALED ON OUTPUT +C TO REPRESENT THE ESTIMATED LOCAL ERROR IN Y +C ON THE LAST STEP. THIS IS THE VECTOR E IN +C THE DESCRIPTION OF THE ERROR CONTROL. +C IT IS DEFINED ONLY ON A SUCCESSFUL RETURN +C FROM ODESSA. +C NRS LENIW-NPAR ARRAY OF SIZE NPAR+1, USED TO STORE THE +C IN IWORK ACCUMULATED NUMBER OF REPEATED STEPS DUE TO +C THE SENSITIVITY ANALYSIS.. +C NRS(1) = TOTAL NUMBER OF REPEATED STEPS, +C NRS(2),... = NUMBER OF REPEATED STEPS DUE TO +C MODEL PARAMETER 1,... +C +C---------------------------------------------------------------------- +C PART II. OTHER ROUTINES CALLABLE. +C +C THE FOLLOWING ARE OPTIONAL CALLS WHICH THE USER MAY MAKE TO +C GAIN ADDITIONAL CAPABILITIES IN CONJUNCTION WITH ODESSA. +C (THE ROUTINES XSETUN AND ODESSA_XSETF ARE DESIGNED TO CONFORM TO THE +C SLATEC ERROR HANDLING PACKAGE.) +C +C FORM OF CALL FUNCTION +C CALL XSETUN(LUN) SET THE LOGICAL UNIT NUMBER, LUN, FOR +C OUTPUT OF MESSAGES FROM ODESSA, IF +C THE DEFAULT IS NOT DESIRED. +C THE DEFAULT VALUE OF LUN IS 6. +C +C CALL ODESSA_XSETF(MFLAG) SET A FLAG TO CONTROL THE PRINTING OF +C MESSAGES BY ODESSA.. +C MFLAG = 0 MEANS DO NOT PRINT. (DANGER.. +C THIS RISKS LOSING VALUABLE INFORMATION.) +C MFLAG = 1 MEANS PRINT (THE DEFAULT). +C +C EITHER OF THE ABOVE CALLS MAY BE MADE AT +C ANY TIME AND WILL TAKE EFFECT IMMEDIATELY. +C +C CALL ODESSA_SVCOM (RSAV, ISAV) STORE IN RSAV AND ISAV THE CONTENTS +C OF THE INTERNAL COMMON BLOCKS USED BY +C ODESSA (SEE PART III BELOW). +C RSAV MUST BE A REAL ARRAY OF LENGTH 222 +C OR MORE, AND ISAV MUST BE AN INTEGER +C ARRAY OF LENGTH 54 OR MORE. +C +C CALL ODESSA_RSCOM (RSAV, ISAV) RESTORE, FROM RSAV AND ISAV, THE CONTENTS +C OF THE INTERNAL COMMON BLOCKS USED BY +C ODESSA. PRESUMES A PRIOR CALL TO ODESSA_SVCOM +C WITH THE SAME ARGUMENTS. +C +C ODESSA_SVCOM AND ODESSA_RSCOM ARE USEFUL IF +C INTERRUPTING A RUN AND RESTARTING +C LATER, OR ALTERNATING BETWEEN TWO OR +C MORE PROBLEMS SOLVED WITH ODESSA. +C +C CALL ODESSA_INTDY(,,,,,) PROVIDE DERIVATIVES OF Y, OF VARIOUS +C (SEE BELOW) ORDERS, AT A SPECIFIED POINT T, IF +C DESIRED. IT MAY BE CALLED ONLY AFTER +C A SUCCESSFUL RETURN FROM ODESSA. +C +C THE DETAILED INSTRUCTIONS FOR USING ODESSA_INTDY ARE AS FOLLOWS. +C THE FORM OF THE CALL IS.. +C +C CALL ODESSA_INTDY (T, K, RWORK(21), NYH, DKY, IFLAG) +C +C THE INPUT PARAMETERS ARE.. +C +C T = VALUE OF INDEPENDENT VARIABLE WHERE ANSWERS ARE DESIRED +C (NORMALLY THE SAME AS THE T LAST RETURNED BY ODESSA). +C FOR VALID RESULTS, T MUST LIE BETWEEN TCUR - HU AND TCUR. +C (SEE OPTIONAL OUTPUTS FOR TCUR AND HU.) +C K = INTEGER ORDER OF THE DERIVATIVE DESIRED. K MUST SATISFY +C 0 .LE. K .LE. NQCUR, WHERE NQCUR IS THE CURRENT ORDER +C (SEE OPTIONAL OUTPUTS). THE CAPABILITY CORRESPONDING +C TO K = 0, I.E. COMPUTING Y(T), IS ALREADY PROVIDED +C BY ODESSA DIRECTLY. SINCE NQCUR .GE. 1, THE FIRST +C DERIVATIVE DY/DT IS ALWAYS AVAILABLE WITH ODESSA_INTDY. +C RWORK(21) = THE BASE ADDRESS OF THE HISTORY ARRAY YH. +C NYH = COLUMN LENGTH OF YH, EQUAL TO THE TOTAL NUMBER OF +C DEPENDENT VARIABLES. IF ISOPT = 0, NYH = N. IF ISOPT = 1, +C NYH = N * (NPAR + 1). +C +C THE OUTPUT PARAMETERS ARE.. +C +C DKY = A REAL ARRAY OF LENGTH NYH CONTAINING THE COMPUTED VALUE +C OF THE K-TH 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 PART III. COMMON BLOCKS. +C +C IF ODESSA IS TO BE USED IN AN OVERLAY SITUATION, THE USER +C MUST DECLARE, IN THE PRIMARY OVERLAY, THE VARIABLES IN.. +C (1) THE CALL SEQUENCE TO ODESSA, +C (2) THE THREE INTERNAL COMMON BLOCKS +C /ODE001/ OF LENGTH 258 (219 DOUBLE PRECISION WORDS +C FOLLOWED BY 39 INTEGER WORDS), +C /ODE002/ OF LENGTH 14 (3 DOUBLE PRECISION WORDS FOLLOWED +C BY 11 INTEGER WORDS), +C /EH0001/ OF LENGTH 2 (INTEGER WORDS). +C +C IF ODESSA 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 THREE COMMON BLOCKS IN HIS MAIN PROGRAM TO INSURE +C THAT THEIR CONTENTS ARE PRESERVED. +C +C IF THE SOLUTION OF A GIVEN PROBLEM BY ODESSA IS TO BE INTERRUPTED +C AND THEN LATER CONTINUED, SUCH AS WHEN RESTARTING AN INTERRUPTED RUN +C OR ALTERNATING BETWEEN TWO OR MORE PROBLEMS, THE USER SHOULD SAVE, +C FOLLOWING THE RETURN FROM THE LAST ODESSA CALL PRIOR TO THE +C INTERRUPTION, THE CONTENTS OF THE CALL SEQUENCE VARIABLES AND THE +C INTERNAL COMMON BLOCKS, AND LATER RESTORE THESE VALUES BEFORE THE +C NEXT ODESSA CALL FOR THAT PROBLEM. TO SAVE AND RESTORE THE COMMON +C BLOCKS, USE SUBROUTINES ODESSA_SVCOM AND ODESSA_RSCOM (SEE PART II ABOVE). +C +C---------------------------------------------------------------------- +C PART IV. OPTIONALLY REPLACEABLE SOLVER ROUTINES. +C +C BELOW ARE DESCRIPTIONS OF TWO ROUTINES IN THE ODESSA PACKAGE WHICH +C RELATE TO THE MEASUREMENT OF ERRORS. EITHER ROUTINE CAN BE +C REPLACED BY A USER-SUPPLIED VERSION, IF DESIRED. HOWEVER, SINCE SUCH +C A REPLACEMENT MAY HAVE A MAJOR IMPACT ON PERFORMANCE, IT SHOULD BE +C DONE ONLY WHEN ABSOLUTELY NECESSARY, AND ONLY WITH GREAT CAUTION. +C (NOTE.. THE MEANS BY WHICH THE PACKAGE VERSION OF A ROUTINE IS +C SUPERSEDED BY THE USER-S VERSION MAY BE SYSTEM-DEPENDENT.) +C +C (A) ODESSA_EWSET. +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 SUBROUTINE ODESSA_EWSET (NYH, ITOL, RTOL, ATOL, YCUR, EWT) +C WHERE NEQ, ITOL, RTOL, AND ATOL ARE AS IN THE ODESSA CALL SEQUENCE, +C YCUR CONTAINS THE CURRENT DEPENDENT VARIABLE VECTOR, AND +C EWT IS THE ARRAY OF WEIGHTS SET BY ODESSA_EWSET. +C +C IF THE USER SUPPLIES THIS SUBROUTINE, IT MUST RETURN IN EWT(I) +C (I = 1,...,NYH) A POSITIVE QUANTITY SUITABLE FOR COMPARING ERRORS +C IN Y(I) TO. THE EWT ARRAY RETURNED BY ODESSA_EWSET IS PASSED TO THE +C ODESSA_VNORM ROUTINE (SEE BELOW), AND ALSO USED BY ODESSA IN THE COMPUTATION +C OF THE OPTIONAL OUTPUT IMXER, THE DIAGONAL JACOBIAN APPROXIMATION, +C AND THE INCREMENTS FOR DIFFERENCE QUOTIENT JACOBIANS. +C +C IN THE USER-SUPPLIED VERSION OF ODESSA_EWSET, 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 ODESSA_EWSET, 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 THE QUANTITIES NQ, NYH, H, AND NST CAN BE OBTAINED BY INCLUDING +C IN ODESSA_EWSET THE STATEMENTS.. +C DOUBLE PRECISION H, RLS +C COMMON /ODE001/ RLS(219),ILS(39) +C NQ = ILS(35) +C NYH = ILS(14) +C NST = ILS(36) +C H = RLS(213) +C THUS, FOR EXAMPLE, THE CURRENT VALUE OF DY/DT CAN BE OBTAINED AS +C YCUR(NYH+I)/H (I=1,...,N) (AND THE DIVISION BY H IS +C UNNECESSARY WHEN NST = 0). +C +C (B) ODESSA_VNORM. +C THE FOLLOWING IS A REAL FUNCTION ROUTINE WHICH COMPUTES THE WEIGHTED +C ROOT-MEAN-SQUARE NORM OF A VECTOR V.. +C D = ODESSA_VNORM (LV, V, W) +C WHERE.. +C LV = 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 ODESSA_VNORM IS CALLED WITH LV = N AND WITH W(I) = 1.0/EWT(I), WHERE +C EWT IS AS SET BY SUBROUTINE ODESSA_EWSET. +C +C IF THE USER SUPPLIES THIS FUNCTION, IT SHOULD RETURN A NON-NEGATIVE +C VALUE OF ODESSA_VNORM SUITABLE FOR USE IN THE ERROR CONTROL IN ODESSA. +C NONE OF THE ARGUMENTS SHOULD BE ALTERED BY ODESSA_VNORM. +C FOR EXAMPLE, A USER-SUPPLIED ODESSA_VNORM 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 OTHER ROUTINES IN THE ODESSA PACKAGE. +C +C IN ADDITION TO SUBROUTINE ODESSA, THE ODESSA PACKAGE INCLUDES THE +C FOLLOWING SUBROUTINES AND FUNCTION ROUTINES.. +C ODESSA_INTDY COMPUTES AN INTERPOLATED VALUE OF THE Y VECTOR AT T = TOUT. +C ODESSA_STODE IS THE CORE INTEGRATOR, WHICH DOES ONE STEP OF THE +C INTEGRATION AND THE ASSOCIATED ERROR CONTROL. +C ODESSA_STESA MANAGES THE SOLUTION OF THE SENSITIVITY FUNCTIONS. +C ODESSA_CFODE SETS ALL METHOD COEFFICIENTS AND TEST CONSTANTS. +C ODESSA_PREPJ COMPUTES AND PREPROCESSES THE JACOBIAN MATRIX J = DF/DY +C AND THE NEWTON ITERATION MATRIX P = I - H*L0*J. +C IT IS ALSO CALLED BY ODESSA_SPRIME (WITH JOPT = 1) TO JUST +C COMPUTE THE JACOBIAN MATRIX. +C ODESSA_PREPDF COMPUTES THE INHOMOGENEITY MATRIX DF/DP. +C ODESSA_SPRIME DEFINES THE SYSTEM OF SENSITIVITY EQUATIONS. +C ODESSA_SOLSY MANAGES SOLUTION OF LINEAR SYSTEM IN CHORD ITERATION. +C ODESSA_EWSET SETS THE ERROR WEIGHT VECTOR EWT BEFORE EACH STEP. +C ODESSA_VNORM COMPUTES THE WEIGHTED R.M.S. NORM OF A VECTOR. +C ODESSA_SVCOM AND ODESSA_RSCOM ARE USER-CALLABLE ROUTINES TO SAVE AND RESTORE, +C RESPECTIVELY, THE CONTENTS OF THE INTERNAL COMMON BLOCKS. +C DGEFA AND DGESL ARE ROUTINES FROM LINPACK FOR SOLVING FULL +C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS. +C DGBFA AND DGBSL ARE ROUTINES FROM LINPACK FOR SOLVING BANDED +C LINEAR SYSTEMS. +C DAXPY, DSCAL, IDAMAX, AND DDOT ARE BASIC LINEAR ALGEBRA MODULES +C (BLAS) USED BY THE ABOVE LINPACK ROUTINES. +C D1MACH COMPUTES THE UNIT ROUNDOFF IN A MACHINE-INDEPENDENT MANNER. +C XERR, XSETUN, AND ODESSA_XSETF HANDLE THE PRINTING OF ALL ERROR +C MESSAGES AND WARNINGS. +C NOTE.. ODESSA_VNORM, IDAMAX, DDOT, AND D1MACH ARE FUNCTION ROUTINES. +C ALL THE OTHERS ARE SUBROUTINES. +C +C THE FORTRAN GENERIC INTRINSIC FUNCTIONS USED BY ODESSA ARE.. +C ABS, MAX, MIN, REAL, MOD, SIGN, SQRT, AND WRITE +C +C A BLOCK DATA SUBPROGRAM IS ALSO INCLUDED WITH THE PACKAGE, +C FOR LOADING SOME OF THE VARIABLES IN INTERNAL COMMON. +C +C---------------------------------------------------------------------- +C PART V. GENERAL REMARKS +C +C THIS SECTION HIGHLIGHTS THE BASIC DIFFERENCES BETWEEN THE ORIGINAL +C LSODE PACKAGE AND THE ODESSA MODIFICATION. THIS IS PROVIDED AS A +C SERVICE TO EXPERIENCED LSODE USERS TO EXPEDITE FAMILIARIZATION WITH +C ODESSA. +C +C (A). ORIGINAL SUBROUTINES AND FUNCTIONS. +C +C OF THE ORIGINAL 22 SUBROUTINES AND FUNCTIONS USED IN THE LSODE +C PACKAGE, ALL ARE USED BY ODESSA, WITH THE FOLLOWING HAVING BEEN +C MODIFIED.. +C +C LSODE THE ORIGINAL DRIVER SUBROUTINE FOR THE LSODE PACKAGE IS +C EXTENSIVELY MODIFIED AND RENAMED ODESSA, WHICH NOW +C CONTAINS A CALL TO ODESSA_SPRIME TO ESTABLISH INITIAL CONDITIONS +C FOR THE SENSITIVITY CALCULATIONS. +C +C ODESSA_STODE THE ONE STEP INTEGRATOR IS SLIGHTLY MODIFIED AND RETAINS +C ITS ORIGINAL NAME. IT NOW CONTAINS THE CALL TO ODESSA_STESA, +C AND ALSO CALLS ODESSA_SPRIME IF KFLAG .LE. -3. +C +C ODESSA_PREPJ ALSO NAMED ODESSA_PREPJ IN ODESSA IS SLIGHTLY MODIFIED TO ALLOW +C FOR THE CALCULATION OF JACOBIAN WITH NO PREPROCESSING +C (JOPT = 1). +C +C (B). NEW SUBROUTINES. +C +C IN ADDITION TO THE CHANGES NOTED ABOVE, THREE NEW SUBROUTINES +C HAVE BEEN INTRODUCED (SEE ODESSA_STESA, ODESSA_SPRIME, AND ODESSA_PREPDF AS DESCRIBED +C IN PART IV. ABOVE). +C +C (C). COMMON BLOCKS. +C +C /LS0001/ RETAINS THE SAME LENGTH AND IS RENAMED /ODE001/; +C HOWEVER THE REAL ARRAY ROWNS(209) IS SHORTENED TO A +C LENGTH OF (173) REAL WORDS, ALLOWING THE REMOVAL OF +C TESCO(3,12) WHICH IS NOW PASSED FROM ODESSA_STODE TO ODESSA_STESA. +C IN ADDITION, THE INTEGER ARRAY IOWNS(6) IS SHORTENED +C TO A LENGTH OF (4) INTEGER WORDS, ALLOWING THE REMOVAL +C OF IALTH AND LMAX WHICH ARE NOW PASSED FROM ODESSA_STODE TO +C ODESSA_STESA. +C +C /ODE002/ ADDED COMMON BLOCK FOR VARIABLES IMPORTANT TO +C SENSITIVITY ANALYSIS (SEE PART III. ABOVE). A BLOCK +C DATA PROGRAM IS NOT REQUIRED FOR THIS COMMON BLOCK. +C +C ODESSA_SVCOM,ODESSA_RSCOM THESE TWO SUBROUTINES ARE MODIFIED TO HANDLE +C COMMON BLOCK /ODE002/ AS WELL. +C +C (D). OPTIONAL INPUTS. +C +C THE FULL SET OF OPTIONAL INPUTS AVAILABLE IN LSODE IS ALSO +C AVAILABLE IN ODESSA, WITH THE EXCEPTION THAT THE NUMBER OF ODE'S +C IN THE MODEL (NEQ(1)), MAY NOT BE CHANGED DURING THE PROBLEM. +C IN ODESSA, NYH NOW REFERS TO THE TOTAL NUMBER OF FIRST-ORDER +C ODE'S (MODEL AND SENSITIVITY EQUATIONS) WHICH IS EQUAL TO +C NEQ(1) IF ISOPT = 0, OR NEQ(1)*(NEQ(2)+1) IF ISOPT = 1. +C NEQ(1), NEQ(2), AND NYH ARE NOT ALLOWED TO CHANGE DURING +C THE COURSE OF AN INTEGRATION. +C +C (E). OPTIONAL OUTPUTS. +C +C THE FULL SET OF OPTIONAL OUTPUTS AVAILABLE IN LSODE IS ALSO +C AVAILABLE IN ODESSA. IN ADDITION, IWORK(19) AND IWORK(20) ARE +C LOADED WITH NDFE AND NSPE, RESPECTIVELY, UPON OUTPUT. THE TOTAL +C NUMBER OF LU DECOMPOSITIONS OF THE PROCESSED JACOBIAN IS EQUAL +C TO NJE - NSPE. +C----------------------------------------------------------------------- + SUBROUTINE ODESSA (F, DF, NEQ, Y, PAR, T, TOUT, ITOL, RTOL, ATOL, + 1 ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + LOGICAL IHIT + EXTERNAL F, DF, JAC, ODESSA_PREPJ, ODESSA_SOLSY, ODESSA_PREPDF + DIMENSION NEQ(*), Y(*), PAR(*), RTOL(*), ATOL(*), IOPT(*), + 1 RWORK(LRW), IWORK(LIW), MORD(2) +C----------------------------------------------------------------------- +C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA.. +C AN ORDINARY DIFFERENTIAL EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS +C SENSITIVITY ANALYSIS. +C +C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF +C LSODE.. LIVERMORE SOLVER FOR ORDINARY DIFFERENTIAL EQUATIONS. +C THIS VERSION IS IN DOUBLE PRECISION. +C +C ODESSA SOLVES FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS.. +C DY(I)/DP, FOR A SINGLE PARAMETER, OR, +C DY(I)/DP(J), FOR MULTIPLE PARAMETERS, +C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.. +C DY(T)/DT = F(Y,T;P). +C----------------------------------------------------------------------- +C REFERENCES... +C +C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND +C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY +C DIFFERENTIAL EQUATIONS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE, +C (1985). +C +C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY +C DIFFERENTIAL EQUATION SOLVER WITH EXPLICIT SIMULTANEOUS +C SENSITIVITY ANALYSIS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE. +C (1985). +C +C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE +C ORDINARY DIFFERENTIAL EQUATION SOLVERS, ACM-SIGNUM NEWSLETTER, +C VOL. 15, NO. 4 (1980), PP. 10-11. +C----------------------------------------------------------------------- +C THE FOLLOWING INTERNAL COMMON BLOCKS CONTAIN +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 STRUCTURE OF THE BLOCKS ARE AS FOLLOWS.. ALL REAL VARIABLES ARE +C LISTED FIRST, FOLLOWED BY ALL INTEGERS. WITHIN EACH TYPE, THE +C VARIABLES ARE GROUPED WITH THOSE LOCAL TO SUBROUTINE ODESSA FIRST, +C THEN THOSE LOCAL TO SUBROUTINE ODESSA_STODE, AND FINALLY THOSE USED +C FOR COMMUNICATION. THE BLOCKS ARE DECLARED IN SUBROUTINES ODESSA +C ODESSA_INTDY, ODESSA_STODE, ODESSA_STESA, ODESSA_PREPJ, ODESSA_PREPDF, +C AND ODESSA_SOLSY. GROUPS OF VARIABLES ARE REPLACED BY DUMMY ARRAYS IN +C THE COMMON DECLARATIONS IN ROUTINES WHERE THOSE VARIABLES ARE NOT USED. +C----------------------------------------------------------------------- + COMMON /ODE001/ TRET, ROWNS(173), + 1 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 2 ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, + 3 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, NYH, IOWNS(4), + 4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, + 5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + COMMON /ODE002/ DUPS, DSMS, DDNS, + 1 NPAR, LDFDP, LNRS, + 2 ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS + PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0) + 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 APPROPIATELY. +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, JUMP TO BLOCK G AND RETURN IMMEDIATELY. +C----------------------------------------------------------------------- + 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) GO TO 430 + 20 NTREP = 0 +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----------------------------------------------------------------------- + IF (NEQ(1) .LE. 0) GO TO 604 + IF (ISTATE .EQ. 1) GO TO 25 + IF (NEQ(1) .NE. N) GO TO 605 + 25 N = NEQ(1) + IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606 + DO 26 I = 1,3 + 26 IF (IOPT(I) .LT. 0 .OR. IOPT(I) .GT. 1) GO TO 607 + ISOPT = IOPT(2) + IDF = IOPT(3) + NYH = N + NSV = 1 + 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 IF (ISOPT .EQ. 0) GO TO 32 +C CHECK LEGALITY OF THE NON-OPTIONAL INPUTS ISOPT, NPAR. +C COMPUTE NUMBER OF SOLUTION VECTORS AND TOTAL NUMBER OF EQUATIONS. + IF (NEQ(2) .LE. 0) GO TO 628 + IF (ISTATE .EQ. 1) GO TO 31 + IF (NEQ(2) .NE. NPAR) GO TO 629 + 31 NPAR = NEQ(2) + NSV = NPAR + 1 + NYH = NSV * N + IF (MITER .EQ. 0 .OR. MITER .EQ. 3) GO TO 630 +C NEXT PROCESS AND CHECK THE OPTIONAL INPUTS. -------------------------- + 32 IF (IOPT(1) .EQ. 1) GO TO 40 + MAXORD = MORD(METH) + MXSTEP = MXSTP0 + MXHNIL = MXHNL0 + IF (ISTATE .EQ. 1) H0 = ZERO + HMXI = ZERO + HMIN = ZERO + 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. ZERO) GO TO 614 + 50 HMAX = RWORK(6) + IF (HMAX .LT. ZERO) GO TO 615 + HMXI = ZERO + IF (HMAX .GT. ZERO) HMXI = ONE/HMAX + HMIN = RWORK(7) + IF (HMIN .LT. ZERO) 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 WORK SPACE FOR DFDP IS CONTAINED IN ACOR. +C----------------------------------------------------------------------- + 60 LYH = 21 + 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 + NYH + LACOR = LSAVF + N + LDFDP = LACOR + N + LENRW = LACOR + NYH - 1 + IWORK(17) = LENRW + LIWM = 1 + LENIW = 20 + N + IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20 + LNRS = LENIW + LIWM + IF (ISOPT .EQ. 1) LENIW = LNRS + NPAR + 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,NYH + IF (ITOL .GE. 3) RTOLI = RTOL(I) + IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) + IF (RTOLI .LT. ZERO) GO TO 619 + IF (ATOLI .LT. ZERO) GO TO 620 + 70 CONTINUE + IF (ISTATE .EQ. 1) GO TO 100 +C IF ISTATE = 3, SET FLAG TO SIGNAL PARAMETER CHANGES TO ODESSA_STODE. - + 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) + 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 THE INITIAL CALL TO ODESSA_SPRIME IF ISOPT = 1, +C AND THE CALCULATION OF THE INITIAL STEP SIZE. +C THE ERROR WEIGHTS IN EWT ARE INVERTED AFTER BEING LOADED. +C----------------------------------------------------------------------- + 100 UROUND = D1MACH(4) + TN = T + IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 105 + TCRIT = RWORK(1) + IF ((TCRIT - TOUT)*(TOUT - T) .LT. ZERO) GO TO 625 + IF (H0 .NE. ZERO .AND. (T + H0 - TCRIT)*H0 .GT. ZERO) + 1 H0 = TCRIT - T + 105 JSTART = 0 + IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) + NHNIL = 0 + NST = 0 + NJE = 0 + NSLAST = 0 + HU = ZERO + NQU = 0 + CCMAX = 0.3D0 + MAXCOR = 3 + IF (ISOPT .EQ. 1) MAXCOR = 4 + MSBP = 20 + MXNCF = 10 +C INITIAL CALL TO F. (LF0 POINTS TO YH(1,2) AND LOADS IN VALUES). + LF0 = LYH + NYH + CALL F (NEQ, T, Y, PAR, RWORK(LF0)) + NFE = 1 + DUPS = ZERO + DSMS = ZERO + DDNS = ZERO + NDFE = 0 + NSPE = 0 + IF (ISOPT .EQ. 0) GO TO 114 +C INITIALIZE COUNTS FOR REPEATED STEPS DUE TO SENSITIVITY ANALYSIS. + DO 110 J = 1,NSV + 110 IWORK(J + LNRS - 1) = 0 +C LOAD THE INITIAL VALUE VECTOR IN YH. --------------------------------- + 114 DO 115 I = 1,NYH + 115 RWORK(I+LYH-1) = Y(I) +C LOAD AND INVERT THE EWT ARRAY. (H IS TEMPORARILY SET TO ONE.) ------- + NQ = 1 + H = ONE + CALL ODESSA_EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) + DO 120 I = 1,NYH + IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 621 + 120 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) + IF (ISOPT .EQ. 0) GO TO 125 +C CALL ODESSA_SPRIME TO LOAD FIRST-ORDER SENSITIVITY DERIVATIVES INTO +C REMAINING YH(*,2) POSITIONS. + CALL ODESSA_SPRIME (NEQ, Y, RWORK(LYH), NYH, N, NSV, RWORK(LWM), + 1 IWORK(LIWM), RWORK(LEWT), RWORK(LF0), RWORK(LACOR), + 2 RWORK(LDFDP), PAR, F, JAC, DF, ODESSA_PREPJ, ODESSA_PREPDF) + IF (IERSP .EQ. -1) GO TO 631 + IF (IERSP .EQ. -2) GO TO 632 +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. ONLY THE ORIGINAL +C SOLUTION VECTOR IS CONSIDERED IN THIS CALCULATION (ISOPT = 0 OR 1). +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----------------------------------------------------------------------- + 125 IF (H0 .NE. ZERO) GO TO 180 + TDIST = ABS(TOUT - T) + W0 = MAX(ABS(T),ABS(TOUT)) + IF (TDIST .LT. TWO*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. ZERO) 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. ZERO) TOL = MAX(TOL,ATOLI/AYI) + 150 CONTINUE + 160 TOL = MAX(TOL,100.0D0*UROUND) + TOL = MIN(TOL,0.001D0) + SUM = ODESSA_VNORM (N, RWORK(LF0), RWORK(LEWT)) + SUM = ONE/(TOL*W0*W0) + TOL*SUM**2 + H0 = ONE/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. ONE) H0 = H0/RH +C LOAD H WITH H0 AND SCALE YH(*,2) BY H0. ------------------------------ + H = H0 + DO 190 I = 1,NYH + 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. ZERO) GO TO 250 + CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) + IF (IFLAG .NE. 0) GO TO 627 + T = TOUT + GO TO 420 + 220 TP = TN - HU*(ONE + 100.0D0*UROUND) + IF ((TP - TOUT)*H .GT. ZERO) GO TO 623 + IF ((TN - TOUT)*H .LT. ZERO) GO TO 250 + GO TO 400 + 230 TCRIT = RWORK(1) + IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624 + IF ((TCRIT - TOUT)*H .LT. ZERO) GO TO 625 + IF ((TN - TOUT)*H .LT. ZERO) GO TO 245 + CALL ODESSA_INTDY (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. ZERO) GO TO 624 + 245 HMX = ABS(TN) + ABS(H) + IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX + IF (IHIT) GO TO 400 + TNEXT = TN + H*(ONE + FOUR*UROUND) + IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250 + H = (TCRIT - TN)*(ONE - FOUR*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 ODESSA_STODE. +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 TOLSF IS CALCULATED CONSIDERING ALL SOLUTION VECTORS. +C----------------------------------------------------------------------- + 250 CONTINUE + IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500 + CALL ODESSA_EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) + DO 260 I = 1,NYH + IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 510 + 260 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) + 270 TOLSF = UROUND*ODESSA_VNORM (NYH, RWORK(LYH), RWORK(LEWT)) + IF (TOLSF .LE. ONE) GO TO 280 + TOLSF = TOLSF*2.0D0 + IF (NST .EQ. 0) GO TO 626 + GO TO 520 + 280 IF (ODESSA_ADDX(TN,H) .NE. TN) GO TO 290 + NHNIL = NHNIL + 1 + IF (NHNIL .GT. MXHNIL) GO TO 290 + CALL XERR ('ODESSA - WARNING..INTERNAL T (=R1) AND H (=R2) ARE', + 1 101, 1, 0, 0, 0, 0, ZERO, ZERO) + CALL XERR ('SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP', + 1 101, 1, 0, 0, 0, 0, ZERO, ZERO) + CALL XERR ('(H = STEP SIZE). SOLVER WILL CONTINUE ANYWAY', + 1 101, 1, 0, 0, 0, 2, TN, H) + IF (NHNIL .LT. MXHNIL) GO TO 290 + CALL XERR ('ODESSA - ABOVE WARNING HAS BEEN ISSUED I1 TIMES.', + 1 102, 1, 0, 0, 0, 0, ZERO, ZERO) + CALL XERR ('IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM', + 1 102, 1, 1, MXHNIL, 0, 0, ZERO,ZERO) + 290 CONTINUE +C----------------------------------------------------------------------- +C CALL ODESSA_STODE(NEQ,Y,YH,NYH,YH,WM,IWM,EWT,SAVF,ACOR,PAR,NRS, +C 1 F,JAC,DF,ODESSA_PREPJ,ODESSA_PREPDF,ODESSA_SOLSY) +C----------------------------------------------------------------------- + CALL ODESSA_STODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), + 1 RWORK(LWM), IWORK(LIWM), RWORK(LEWT), RWORK(LSAVF), + 2 RWORK(LACOR), PAR, IWORK(LNRS), F, JAC, DF, ODESSA_PREPJ, + 3 ODESSA_PREPDF, ODESSA_SOLSY) + KGO = 1 - KFLAG + GO TO (300, 530, 540, 633), 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. ZERO) GO TO 250 + CALL ODESSA_INTDY (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. ZERO) 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. ZERO) GO TO 345 + CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) + T = TOUT + GO TO 420 + 345 HMX = ABS(TN) + ABS(H) + IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX + IF (IHIT) GO TO 400 + TNEXT = TN + H*(ONE + FOUR*UROUND) + IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250 + H = (TCRIT - TN)*(ONE - FOUR*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.0D0*UROUND*HMX +C----------------------------------------------------------------------- +C BLOCK G. +C THE FOLLOWING BLOCK HANDLES ALL SUCCESSFUL RETURNS FROM ODESSA. +C IF ITASK .NE. 1, Y IS LOADED FROM YH AND T IS SET ACCORDINGLY. +C ISTATE IS SET TO 2, THE ILLEGAL INPUT COUNTER IS ZEROED, AND THE +C OPTIONAL OUTPUTS ARE LOADED INTO THE WORK ARRAYS BEFORE RETURNING. +C IF ISTATE = 1 AND TOUT = T, THERE IS A RETURN WITH NO ACTION TAKEN, +C EXCEPT THAT IF THIS HAS HAPPENED REPEATEDLY, THE RUN IS TERMINATED. +C----------------------------------------------------------------------- + 400 DO 410 I = 1,NYH + 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 + ILLIN = 0 + RWORK(11) = HU + RWORK(12) = H + RWORK(13) = TN + IWORK(11) = NST + IWORK(12) = NFE + IWORK(13) = NJE + IWORK(14) = NQU + IWORK(15) = NQ + IF (ISOPT .EQ. 0) RETURN + IWORK(19) = NDFE + IWORK(20) = NSPE + RETURN + 430 NTREP = NTREP + 1 + IF (NTREP .LT. 5) RETURN + CALL XERR ('ODESSA -- REPEATED CALLS WITH ISTATE = 1 AND + 1TOUT = T (=R1)', 301, 1, 0, 0, 0, 1, T, ZERO) + GO TO 800 +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, T IS SET TO TN, AND THE ILLEGAL INPUT +C COUNTER ILLIN IS SET TO 0. THE OPTIONAL OUTPUTS ARE LOADED INTO +C THE WORK ARRAYS BEFORE RETURNING. +C----------------------------------------------------------------------- +C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE REACHING TOUT. ---------- + 500 CALL XERR ('ODESSA - AT CURRENT T (=R1), MXSTEP (=I1) STEPS', + 1 201, 1, 0, 0, 0, 0, ZERO,ZERO) + CALL XERR ('TAKEN ON THIS CALL BEFORE REACHING TOUT', + 1 201, 1, 1, MXSTEP, 0, 1, TN, ZERO) + 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) + CALL XERR ('ODESSA - AT T (=R1), EWT(I1) HAS BECOME R2 .LE. 0.', + 1 202, 1, 1, I, 0, 2, TN, EWTI) + ISTATE = -6 + GO TO 580 +C TOO MUCH ACCURACY REQUESTED FOR MACHINE PRECISION. ------------------- + 520 CALL XERR ('ODESSA - AT T (=R1), TOO MUCH ACCURACY REQUESTED', + 1 203, 1, 0, 0, 0, 0, ZERO,ZERO) + CALL XERR ('FOR PRECISION OF MACHINE.. SEE TOLSF (=R2)', + 1 203, 1, 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 CALL XERR ('ODESSA - AT T(=R1) AND STEP SIZE H(=R2), THE ERROR', + 1 204, 1, 0, 0, 0, 0, ZERO, ZERO) + CALL XERR ('TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN', + 1 204, 1, 0, 0, 0, 2, TN, H) + ISTATE = -4 + GO TO 560 +C KFLAG = -2. CONVERGENCE FAILED REPEATEDLY OR WITH ABS(H) = HMIN. ---- + 540 CALL XERR ('ODESSA - AT T (=R1) AND STEP SIZE H (=R2), THE', + 1 205, 1, 0, 0, 0, 0, ZERO,ZERO) + CALL XERR ('CORRECTOR CONVERGENCE FAILED REPEATEDLY', + 1 205, 1, 0, 0, 0, 0, ZERO,ZERO) + CALL XERR ('OR WITH ABS(H) = HMIN', + 1 205, 1, 0, 0, 0, 2, TN, H) + ISTATE = -5 +C COMPUTE IMXER IF RELEVANT. ------------------------------------------- + 560 BIG = ZERO + IMXER = 1 + DO 570 I = 1,NYH + 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, ILLIN, AND OPTIONAL OUTPUTS. ------------------------ + 580 DO 590 I = 1,NYH + 590 Y(I) = RWORK(I+LYH-1) + T = TN + ILLIN = 0 + RWORK(11) = HU + RWORK(12) = H + RWORK(13) = TN + IWORK(11) = NST + IWORK(12) = NFE + IWORK(13) = NJE + IWORK(14) = NQU + IWORK(15) = NQ + IF (ISOPT .EQ. 0) RETURN + IWORK(19) = NDFE + IWORK(20) = NSPE + 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. THEN IF THERE HAVE BEEN +C 5 CONSECUTIVE SUCH RETURNS JUST BEFORE THIS CALL TO THE SOLVER, +C THE RUN IS HALTED. +C----------------------------------------------------------------------- + 601 CALL XERR ('ODESSA - ISTATE (=I1) ILLEGAL', + 1 1, 1, 1, ISTATE, 0, 0, ZERO,ZERO) + GO TO 700 + 602 CALL XERR ('ODESSA - ITASK (=I1) ILLEGAL', + 1 2, 1, 1, ITASK, 0, 0, ZERO,ZERO) + GO TO 700 + 603 CALL XERR ('ODESSA - ISTATE .GT. 1 BUT ODESSA NOT INITIALIZED', + 1 3, 1, 0, 0, 0, 0, ZERO,ZERO) + GO TO 700 + 604 CALL XERR ('ODESSA - NEQ (=I1) .LT. 1', + 1 4, 1, 1, NEQ(1), 0, 0, ZERO,ZERO) + GO TO 700 + 605 CALL XERR ('ODESSA - ISTATE = 3 AND NEQ CHANGED. (I1 TO I2)', + 1 5, 1, 2, N, NEQ(1), 0, ZERO,ZERO) + GO TO 700 + 606 CALL XERR ('ODESSA - ITOL (=I1) ILLEGAL', + 1 6, 1, 1, ITOL, 0, 0, ZERO,ZERO) + GO TO 700 + 607 CALL XERR ('ODESSA - IOPT (=I1) ILLEGAL', + 1 7, 1, 1, IOPT, 0, 0, ZERO,ZERO) + GO TO 700 + 608 CALL XERR('ODESSA - MF (=I1) ILLEGAL', + 1 8, 1, 1, MF, 0, 0, ZERO,ZERO) + GO TO 700 + 609 CALL XERR('ODESSA - ML (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)', + 1 9, 1, 2, ML, NEQ(1), 0, ZERO,ZERO) + GO TO 700 + 610 CALL XERR('ODESSA - MU (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)', + 1 10, 1, 2, MU, NEQ(1), 0, ZERO,ZERO) + GO TO 700 + 611 CALL XERR('ODESSA - MAXORD (=I1) .LT. 0', + 1 11, 1, 1, MAXORD, 0, 0, ZERO,ZERO) + GO TO 700 + 612 CALL XERR('ODESSA - MXSTEP (=I1) .LT. 0', + 1 12, 1, 1, MXSTEP, 0, 0, ZERO,ZERO) + GO TO 700 + 613 CALL XERR('ODESSA - MXHNIL (=I1) .LT. 0', + 1 13, 1, 1, MXHNIL, 0, 0, ZERO,ZERO) + GO TO 700 + 614 CALL XERR('ODESSA - TOUT (=R1) BEHIND T (=R2)', + 1 14, 1, 0, 0, 0, 2, TOUT, T) + CALL XERR('INTEGRATION DIRECTION IS GIVEN BY H0 (=R1)', + 1 14, 1, 0, 0, 0, 1, H0, ZERO) + GO TO 700 + 615 CALL XERR('ODESSA - HMAX (=R1) .LT. 0.0', + 1 15, 1, 0, 0, 0, 1, HMAX, ZERO) + GO TO 700 + 616 CALL XERR('ODESSA - HMIN (=R1) .LT. 0.0', + 1 16, 1, 0, 0, 0, 1, HMIN, ZERO) + GO TO 700 + 617 CALL XERR('ODESSA - RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS + 1 LRW (=I2)', 17, 1, 2, LENRW, LRW, 0, ZERO,ZERO) + GO TO 700 + 618 CALL XERR('ODESSA - IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS + 1 LIW (=I2)', 18, 1, 2, LENIW, LIW, 0, ZERO,ZERO) + GO TO 700 + 619 CALL XERR('ODESSA - RTOL(I1) IS R1 .LT. 0.0', + 1 19, 1, 1, I, 0, 1, RTOLI, ZREO) + GO TO 700 + 620 CALL XERR('ODESSA - ATOL(I1) IS R1 .LT. 0.0', + 1 20, 1, 1, I, 0, 1, ATOLI, ZERO) + GO TO 700 +* + 621 EWTI = RWORK(LEWT+I-1) + CALL XERR('ODESSA - EWT(I1) IS R1 .LE. 0.0', + 1 21, 1, 1, I, 0, 1, EWTI, ZERO) + GO TO 700 + 622 CALL XERR('ODESSA - TOUT (=R1) TOO CLOSE TO T(=R2) TO START + 1 INTEGRATION', 22, 1, 0, 0, 0, 2, TOUT, T) + GO TO 700 + 623 CALL XERR('ODESSA - ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU + 1 (= R2)', 23, 1, 1, ITASK, 0, 2, TOUT, TP) + GO TO 700 + 624 CALL XERR('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR + 1 (=R2)', 24, 1, 0, 0, 0, 2, TCRIT, TN) + GO TO 700 + 625 CALL XERR('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT + 1 (=R2)', 25, 1, 0, 0, 0, 2, TCRIT, TOUT) + GO TO 700 + 626 CALL XERR('ODESSA - AT START OF PROBLEM, TOO MUCH ACCURACY', + 1 26, 1, 0, 0, 0, 0, ZERO,ZERO) + CALL XERR('REQUESTED FOR PRECISION OF MACHINE. SEE TOLSF (=R1)', + 1 26, 1, 0, 0, 0, 1, TOLSF, ZERO) + RWORK(14) = TOLSF + GO TO 700 + 627 CALL XERR('ODESSA - TROUBLE FROM ODESSA_INTDY. ITASK = I1, TOUT = + 1 R1', 27, 1, 1, ITASK, 0, 1, TOUT, ZERO) + GO TO 700 +C ERROR STATEMENTS ASSOCIATED WITH SENSITIVITY ANALYSIS. + 628 CALL XERR('ODESSA - NPAR (=I1) .LT. 1', + 1 28, 1, 1, NPAR, 0, 0, ZERO,ZERO) + GO TO 700 + 629 CALL XERR('ODESSA - ISTATE = 3 AND NPAR CHANGED (I1 TO I2)', + 1 29, 1, 2, NP, NPAR, 0, ZERO,ZERO) + GO TO 700 + 630 CALL XERR('ODESSA - MITER (=I1) ILLEGAL', + 1 30, 1, 1, MITER, 0, 0, ZERO,ZERO) + GO TO 700 + 631 CALL XERR('ODESSA - TROUBLE IN ODESSA_SPRIME (IERPJ)', + 1 31, 1, 0, 0, 0, 0, ZERO,ZERO) + GO TO 700 + 632 CALL XERR('ODESSA - TROUBLE IN ODESSA_SPRIME (MITER)', + 1 32, 1, 0, 0, 0, 0, ZERO,ZERO) + GO TO 700 + 633 CALL XERR('ODESSA - FATAL ERROR IN ODESSA_STODE (KFLAG = -3)', + 1 33, 2, 0, 0, 0, 0, ZERO,ZERO) + GO TO 801 +C + 700 IF (ILLIN .EQ. 5) GO TO 710 + ILLIN = ILLIN + 1 + ISTATE = -3 + RETURN + 710 CALL XERR('ODESSA - REPEATED OCCURRENCES OF ILLEGAL INPUT', + 1 302, 1, 0, 0, 0, 0, ZERO,ZERO) +C + 800 CALL XERR('ODESSA - RUN ABORTED.. APPARENT INFINITE LOOP', + 1 303, 2, 0, 0, 0, 0, ZERO,ZERO) + RETURN + 801 CALL XERR('ODESSA - RUN ABORTED', + 1 304, 2, 0, 0, 0, 0, ZERO,ZERO) + RETURN +C-------------------- END OF SUBROUTINE ODESSA ------------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_addx.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,11 @@ + DOUBLE PRECISION FUNCTION ODESSA_ADDX(A,B) + DOUBLE PRECISION A,B +C +C THIS FUNCTION IS NECESSARY TO FORCE OPTIMIZING COMPILERS TO +C EXECUTE AND STORE A SUM, FOR SUCCESSFUL EXECUTION OF THE +C TEST A + B = B. +C + ODESSA_ADDX = A + B + RETURN +C-------------------- END OF FUNCTION SUM ------------------------------ + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_cfode.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,108 @@ + SUBROUTINE ODESSA_CFODE (METH, ELCO, TESCO) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION ELCO(13,12), TESCO(3,12) +C----------------------------------------------------------------------- +C ODESSA_CFODE 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 ODESSA_CFODE 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----------------------------------------------------------------------- + DIMENSION PC(12) + PARAMETER (ONE=1.0D0,ZERO=0.0D0) +C + GO TO (100, 200), METH +C + 100 ELCO(1,1) = ONE + ELCO(2,1) = ONE + TESCO(1,1) = ZERO + TESCO(2,1) = 2.0D0 + TESCO(1,2) = ONE + TESCO(3,12) = ZERO + PC(1) = ONE + RQFAC = ONE + 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/REAL(NQ) + NQM1 = NQ - 1 + FNQM1 = REAL(NQM1) + NQP1 = NQ + 1 +C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ---------------------------------- + PC(NQ) = ZERO + 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.0D0 + TSIGN = ONE + DO 120 I = 2,NQ + TSIGN = -TSIGN + PINT = PINT + TSIGN*PC(I)/REAL(I) + 120 XPIN = XPIN + TSIGN*PC(I)/REAL(I+1) +C STORE COEFFICIENTS IN ELCO AND TESCO. -------------------------------- + ELCO(1,NQ) = PINT*RQ1FAC + ELCO(2,NQ) = ONE + DO 130 I = 2,NQ + 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/REAL(I) + AGAMQ = RQFAC*XPIN + RAGQ = ONE/AGAMQ + TESCO(2,NQ) = RAGQ + IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/REAL(NQP1) + TESCO(3,NQM1) = RAGQ + 140 CONTINUE + RETURN +C + 200 PC(1) = ONE + RQ1FAC = ONE + 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 = REAL(NQ) + NQP1 = NQ + 1 +C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------ + PC(NQP1) = ZERO + 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) = ONE + TESCO(1,NQ) = RQ1FAC + TESCO(2,NQ) = REAL(NQP1)/ELCO(1,NQ) + TESCO(3,NQ) = REAL(NQ+2)/ELCO(1,NQ) + RQ1FAC = RQ1FAC/FNQ + 230 CONTINUE + RETURN +C----------------------- END OF SUBROUTINE ODESSA_CFODE ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_ewset.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,19 @@ + SUBROUTINE ODESSA_EWSET (N, ITOL, RTOL, ATOL, YCUR, EWT) +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----------------------------------------------------------------------- + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION RTOL(*), ATOL(*), YCUR(N), EWT(N) + RTOLI = RTOL(1) + ATOLI = ATOL(1) + DO 10 I = 1,N + IF (ITOL .GE. 3) RTOLI = RTOL(I) + IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) + EWT(I) = RTOLI*ABS(YCUR(I)) + ATOLI + 10 CONTINUE + RETURN +C----------------------- END OF SUBROUTINE ODESSA_EWSET ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_intdy.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,74 @@ + SUBROUTINE ODESSA_INTDY (T, K, YH, NYH, DKY, IFLAG) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION YH(NYH,1), DKY(1) + COMMON /ODE001/ ROWND, ROWNS(173), + 2 RDUM1(38),H, RDUM2(2), HU, RDUM3, TN, UROUND, + 3 IOWND(14), IOWNS(4), + 4 IDUM1(8), L, IDUM2, + 5 IDUM3(5), N, NQ, IDUM4(4) +C----------------------------------------------------------------------- +C ODESSA_INTDY 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----------------------------------------------------------------------- + IFLAG = 0 + IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80 + TP = TN - HU*(1.0D0 + 100.0D0*UROUND) + IF ((T-TP)*(T-TN) .GT. 0.0D0) 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 = REAL(IC) + DO 20 I = 1,NYH + 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 = REAL(IC) + DO 40 I = 1,NYH + 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,NYH + 60 DKY(I) = R*DKY(I) + RETURN +C + 80 CALL XERR('ODESSA_INTDY-- K (=I1) ILLEGAL', + 1 51, 1, 1, K, 0, 0, ZERO,ZERO) + IFLAG = -1 + RETURN + 90 CALL XERR ('ODESSA_INTDY-- T (=R1) ILLEGAL', + 1 52, 1, 0, 0, 0, 1, T, ZERO) + CALL XERR('T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2)', + 1 52, 1, 0, 0, 0, 2, TP, TN) + IFLAG = -2 + RETURN +C------------------ END OF SUBROUTINE ODESSA_INTDY ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_prepd.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,54 @@ + SUBROUTINE ODESSA_PREPDF (NEQ, Y, SRUR, SAVF, FTEM, DFDP, PAR, + 1 F, DF, JPAR) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + EXTERNAL F, DF + DIMENSION NEQ(*), Y(*), SAVF(*), FTEM(*), DFDP(*), PAR(*) + COMMON /ODE001/ ROWND, ROWNS(173), + 1 RDUM1(43), TN, RDUM2, + 2 IOWND1(14), IOWNS(4), + 3 IDUM1(10), MITER, IDUM2(4), N, IDUM3(2), NFE, IDUM4(2) + COMMON /ODE002/ RDUM3(3), + 1 IOWND2(3), IDUM5(2), NDFE, IDUM6, IDF, IDUM7(3) +C----------------------------------------------------------------------- +C ODESSA_PREPDF IS CALLED BY ODESSA_SPRIME AND ODESSA_STESA TO COMPUTE +C THE INHOMOGENEITY VECTORS DF(I)/DP(JPAR). HERE DF/DP IS COMPUTED BY +C THE USER-SUPPLIED ROUTINE DF IF IDF = 1, OR BY FINITE DIFFERENCING IF +C IDF = 0. +C +C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH +C ODESSA_PREPDF USES THE FOLLOWING.. +C Y = REAL ARRAY OF LENGTH NYH CONTAINING DEPENDENT VARIABLES. +C ODESSA_PREPDF USES ONLY THE FIRST N ENTRIES OF Y(*). +C SRUR = SQRT(UROUND) (= WM(1)). +C SAVF = REAL ARRAY OF LENGTH N CONTAINING DERIVATIVES DY/DT. +C FTEM = REAL ARRAY OF LENGTH N USED TO TEMPORARILY STORE DY/DT FOR +C NUMERICAL DIFFERENTIATION. +C DFDP = REAL ARRAY OF LENGTH N USED TO STORE DF(I)/DP(JPAR), I = 1,N. +C PAR = REAL ARRAY OF LENGTH NPAR CONTAINING EQUATION PARAMETERS +C OF INTEREST. +C JPAR = INPUT PARAMETER, 2 .LE. JPAR .LE. NSV, DESIGNATING THE +C APPROPRIATE SOLUTION VECTOR CORRESPONDING TO PAR(JPAR). +C THIS ROUTINE ALSO USES THE COMMON VARIABLES TN, MITER, N, NFE, NDFE, +C AND IDF. +C----------------------------------------------------------------------- + NDFE = NDFE + 1 + IDF1 = IDF + 1 + GO TO (100, 200), IDF1 +C IDF = 0, CALL F TO APPROXIMATE DFDP. --------------------------------- + 100 RPAR = PAR(JPAR) + R = MAX(SRUR*ABS(RPAR),SRUR) + PAR(JPAR) = RPAR + R + FAC = 1.0D0/R + CALL F (NEQ, TN, Y, PAR, FTEM) + DO 110 I = 1,N + 110 DFDP(I) = (FTEM(I) - SAVF(I))*FAC + PAR(JPAR) = RPAR + NFE = NFE + 1 + RETURN +C IDF = 1, CALL USER SUPPLIED DF. -------------------------------------- + 200 DO 210 I = 1,N + 210 DFDP(I) = 0.0D0 + CALL DF (NEQ, TN, Y, PAR, DFDP, JPAR) + RETURN +C -------------------- END OF SUBROUTINE ODESSA_PREPDF ------------------------ + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_prepj.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,166 @@ + SUBROUTINE ODESSA_PREPJ (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, + 1 FTEM, PAR, F, JAC, JOPT) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION NEQ(*), Y(*), YH(NYH,*), WM(*), IWM(*), EWT(*), + 1 SAVF(*), FTEM(*), PAR(*) + EXTERNAL F, JAC + PARAMETER (ZERO=0.0D0,ONE=1.0D0) + COMMON /ODE001/ ROWND, ROWNS(173), + 2 RDUM1(37), EL0, H, RDUM2(4), TN, UROUND, + 3 IOWND(14), IOWNS(4), + 4 IDUM1(3), IERPJ, IDUM2, JCUR, IDUM3(4), + 5 MITER, IDUM4(4), N, IDUM5(2), NFE, NJE, IDUM6 +C----------------------------------------------------------------------- +C ODESSA_PREPJ IS CALLED BY ODESSA_STODE TO COMPUTE AND PROCESS THE MATRIX +C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. +C IF ISOPT = 1, ODESSA_PREPJ IS ALSO CALLED BY ODESSA_SPRIME WITH JOPT = 1. +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 (JOPT = 0) IN PREPARATION FOR LATER +C SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS +C DONE BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5. +C +C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION +C WITH ODESSA_PREPJ USES THE FOLLOWING.. +C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY. +C FTEM = WORK ARRAY OF LENGTH N (ACOR IN ODESSA_STODE). +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 JOPT = INPUT JACOBIAN OPTION, = 1 IF JAC IS DESIRED ONLY. +C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND, +C IERPJ, JCUR, MITER, N, NFE, AND NJE. +C----------------------------------------------------------------------- + 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) = ZERO + CALL JAC (NEQ, TN, Y, PAR, 0, 0, WM(3), N) + IF (JOPT .EQ. 1) RETURN + 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 = ODESSA_VNORM (N, SAVF, EWT) + R0 = 1000.0D0*ABS(H)*UROUND*REAL(N)*FAC + IF (R0 .EQ. ZERO) R0 = ONE + 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, PAR, 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 + IF (JOPT .EQ. 1) RETURN +C ADD IDENTITY MATRIX. ------------------------------------------------- + 240 J = 3 + DO 250 I = 1,N + WM(J) = WM(J) + ONE + 250 J = J + (N + 1) +C DO LU DECOMPOSITION ON P. -------------------------------------------- + CALL DGEFA (WM(3), N, 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.1D0 + DO 310 I = 1,N + 310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) + CALL F (NEQ, TN, Y, PAR, WM(3)) + NFE = NFE + 1 + DO 320 I = 1,N + R0 = H*SAVF(I) - YH(I,2) + DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I)) + WM(I+2) = 1.0D0 + IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320 + IF (ABS(DI) .EQ. ZERO) GO TO 330 + WM(I+2) = 0.1D0*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) = ZERO + CALL JAC (NEQ, TN, Y, PAR, ML, MU, WM(ML3), MEBAND) + IF (JOPT .EQ. 1) RETURN + 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 = ODESSA_VNORM (N, SAVF, EWT) + R0 = 1000.0D0*ABS(H)*UROUND*REAL(N)*FAC + IF (R0 .EQ. ZERO) R0 = ONE + 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, PAR, 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 + IF (JOPT .EQ. 1) RETURN +C ADD IDENTITY MATRIX. ------------------------------------------------- + 570 II = MBAND + 2 + DO 580 I = 1,N + WM(II) = WM(II) + ONE + 580 II = II + MEBAND +C DO LU DECOMPOSITION OF P. -------------------------------------------- + CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER) + IF (IER .NE. 0) IERPJ = 1 + RETURN +C---------------- END OF SUBROUTINE ODESSA_PREPJ ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_rscom.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,29 @@ + SUBROUTINE ODESSA_RSCOM (RSAV, ISAV) +C----------------------------------------------------------------------- +C THIS ROUTINE RESTORES FROM RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS +C ODE001, ODE002 AND EH0001, WHICH ARE USED INTERNALLY IN THE ODESSSA +C PACKAGE. THIS PRESUMES THAT RSAV AND ISAV WERE LOADED BY MEANS +C OF SUBROUTINE ODESSA_SVCOM OR THE EQUIVALENT. +C----------------------------------------------------------------------- + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION RSAV(*), ISAV(*) + COMMON /ODE001/ RODE1(219), IODE1(39) + COMMON /ODE002/ RODE2(3), IODE2(11) + COMMON /EH0001/ IEH(2) + DATA LRODE1/219/, LIODE1/39/, LRODE2/3/, LIODE2/11/ +C + DO 10 I = 1,LRODE1 + 10 RODE1(I) = RSAV(I) + DO 20 I = 1,LRODE2 + J = LRODE1 + I + 20 RODE2(I) = RSAV(J) + DO 30 I = 1,LIODE1 + 30 IODE1(I) = ISAV(I) + DO 40 I = 1,LODE2 + J = LIODE1 + I + 40 IODE2(I) = ISAV(J) + IEH(1) = ISAV(LIODE1+LIODE2+1) + IEH(2) = ISAV(LIODE1+LIODE2+2) + RETURN +C----------------------- END OF SUBROUTINE ODESSA_RSCOM ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_solsy.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,60 @@ + SUBROUTINE ODESSA_SOLSY (WM, IWM, X, TEM) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION WM(*), IWM(*), X(*), TEM(*) + PARAMETER (ZERO=0.0D0,ONE=1.0D0) + COMMON /ODE001/ ROWND, ROWNS(173), + 2 RDUM1(37), EL0, H, RDUM2(6), + 3 IOWND(14), IOWNS(4), + 4 IDUM1(4), IERSL, IDUM2(5), + 5 MITER, IDUM3(4), N, IDUM4(5) +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 DGESL 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 DGBSL. +C COMMUNICATION WITH ODESSA_SOLSY 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----------------------------------------------------------------------- + IERSL = 0 + GO TO (100, 100, 300, 400, 400), MITER + 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0) + 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 = ONE - R*(ONE - ONE/WM(I+2)) + IF (ABS(DI) .EQ. ZERO) GO TO 390 + 320 WM(I+2) = ONE/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 DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0) + RETURN +C----------------------- END OF SUBROUTINE ODESSA_SOLSY ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_sprim.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,125 @@ + SUBROUTINE ODESSA_SPRIME (NEQ, Y, YH, NYH, NROW, NCOL, WM, IWM, + 1 EWT, SAVF, FTEM, DFDP, PAR, F, JAC, DF, PJAC, PDF) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION NEQ(*), Y(*), YH(NROW,NCOL,*), WM(*), IWM(*), + 1 EWT(*), SAVF(*), FTEM(*), DFDP(NROW,*), PAR(*) + EXTERNAL F, JAC, DF, PJAC, PDF + PARAMETER (ONE=1.0D0,ZERO=0.0D0) + COMMON /ODE001/ ROWND, ROWNS(173), + 1 RDUM1(37),EL0, H, RDUM2(6), + 2 IOWND1(14), IOWNS(4), + 3 IDUM1(3), IERPJ, IDUM2(6), + 4 MITER, IDUM3(4), N, IDUM4(5) + COMMON /ODE002/ RDUM3(3), + 1 IOWND2(3), IDUM5, NSV, IDUM6, NSPE, IDUM7, IERSP, JOPT, IDUM8 +C----------------------------------------------------------------------- +C ODESSA_SPRIME IS CALLED BY ODESSA TO INITIALIZE THE YH ARRAY. IT IS ALSO +C CALLED BY ODESSA_STODE TO REEVALUATE FIRST ORDER DERIVATIVES WHEN KFLAG +C .LE. -3. ODESSA_SPRIME COMPUTES THE FIRST DERIVATIVES OF THE SENSITIVITY +C COEFFICIENTS WITH RESPECT TO THE INDEPENDENT VARIABLE T... +C +C ODESSA_SPRIME = D(DY/DP)/DT = JAC*DY/DP + DF/DP +C WHERE JAC = JACOBIAN MATRIX +C DY/DP = SENSITIVITY MATRIX +C DF/DP = INHOMOGENEITY MATRIX +C THIS ROUTINE USES THE COMMON VARIABLES EL0, H, IERPJ, MITER, N, +C NSV, NSPE, IERSP, JOPT +C----------------------------------------------------------------------- +C CALL ODESSA_PREPJ WITH JOPT = 1. +C IF MITER = 2 OR 5, EL0 IS TEMPORARILY SET TO -1.0 AND H IS +C TEMPORARILY SET TO 1.0D0. +C----------------------------------------------------------------------- + NSPE = NSPE + 1 + JOPT = 1 + IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 10 + HTEMP = H + ETEMP = EL0 + H = ONE + EL0 = -ONE + 10 CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, FTEM, + 1 PAR, F, JAC, JOPT) + IF (IERPJ .NE. 0) GO TO 300 + JOPT = 0 + IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 20 + H = HTEMP + EL0 = ETEMP +C----------------------------------------------------------------------- +C CALL ODESSA_PREPDF AND LOAD DFDP(*,JPAR). +C----------------------------------------------------------------------- + 20 DO 30 J = 2,NSV + JPAR = J - 1 + CALL PDF (NEQ, Y, WM, SAVF, FTEM, DFDP(1,JPAR), PAR, + 1 F, DF, JPAR) + 30 CONTINUE +C----------------------------------------------------------------------- +C COMPUTE JAC*DY/DP AND STORE RESULTS IN YH(*,*,2). +C----------------------------------------------------------------------- + GO TO (40,40,310,100,100) MITER +C THE JACOBIAN IS FULL.------------------------------------------------ +C FOR EACH ROW OF THE JACOBIAN.. + 40 DO 70 IROW = 1,N +C AND EACH COLUMN OF THE SENSITIVITY MATRIX.. + DO 60 J = 2,NSV + SUM = ZERO +C TAKE THE VECTOR DOT PRODUCT.. + DO 50 I = 1,N + IPD = IROW + N*(I-1) + 2 + SUM = SUM + WM(IPD)*YH(I,J,1) + 50 CONTINUE + YH(IROW,J,2) = SUM + 60 CONTINUE + 70 CONTINUE + GO TO 200 +C THE JACOBIAN IS BANDED.----------------------------------------------- + 100 ML = IWM(1) + MU = IWM(2) + ICOUNT = 1 + MBAND = ML + MU + 1 + MEBAND = MBAND + ML + NMU = N - MU + ML1 = ML + 1 +C FOR EACH ROW OF THE JACOBIAN.. + DO 160 IROW = 1,N + IF (IROW .GT. ML1) GO TO 110 + IPD = MBAND + IROW + 1 + IYH = 1 + LBAND = MU + IROW + GO TO 120 + 110 ICOUNT = ICOUNT + 1 + IPD = ICOUNT*MEBAND + 2 + IYH = IYH + 1 + LBAND = LBAND - 1 + IF (IROW .LE. NMU) LBAND = MBAND +C AND EACH COLUMN OF THE SENSITIVITY MATRIX.. + 120 DO 150 J = 2,NSV + SUM = ZERO + I1 = IPD + I2 = IYH +C TAKE THE VECTOR DOT PRODUCT. + DO 140 I = 1,LBAND + SUM = SUM + WM(I1)*YH(I2,J,1) + I1 = I1 + MEBAND - 1 + I2 = I2 + 1 + 140 CONTINUE + YH(IROW,J,2) = SUM + 150 CONTINUE + 160 CONTINUE +C----------------------------------------------------------------------- +C ADD THE INHOMOGENEITY TERM, I.E., ADD DFDP(*,JPAR) TO YH(*,JPAR+1,2). +C----------------------------------------------------------------------- + 200 DO 220 J = 2,NSV + JPAR = J - 1 + DO 210 I = 1,N + YH(I,J,2) = YH(I,J,2) + DFDP(I,JPAR) + 210 CONTINUE + 220 CONTINUE + RETURN +C----------------------------------------------------------------------- +C ERROR RETURNS. +C----------------------------------------------------------------------- + 300 IERSP = -1 + RETURN + 310 IERSP = -2 + RETURN +C------------------------END OF SUBROUTINE ODESSA_SPRIME----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_stesa.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,135 @@ + SUBROUTINE ODESSA_STESA (NEQ, Y, NROW, NCOL, YH, WM, IWM, EWT, + 1 SAVF, ACOR, PAR, NRS, F, JAC, DF, PJAC, PDF, SOLVE) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + EXTERNAL F, JAC, DF, PJAC, PDF, SOLVE + DIMENSION NEQ(*), Y(NROW,*), YH(NROW,NCOL,*), WM(*), IWM(*), + 1 EWT(NROW,*), SAVF(*), ACOR(NROW,*), PAR(*), NRS(*) + PARAMETER (ONE=1.0D0,ZERO=0.0D0) + COMMON /ODE001/ ROWND, ROWNS(173), + 1 TESCO(3,12), RDUM1, EL0, H, RDUM2(4), TN, RDUM3, + 2 IOWND1(14), IOWNS(4), + 3 IALTH, LMAX, IDUM1, IERPJ, IERSL, JCUR, IDUM2, KFLAG, L, IDUM3, + 4 MITER, IDUM4(4), N, NQ, IDUM5, NFE, IDUM6(2) + COMMON /ODE002/ DUPS, DSMS, DDNS, + 1 IOWND2(3), IDUM7, NSV, IDUM8(2), IDF, IDUM9, JOPT, KFLAGS +C----------------------------------------------------------------------- +C ODESSA_STESA IS CALLED BY ODESSA_STODE TO PERFORM AN EXPLICIT +C CALCULATION FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS DY(I)/DP(J), +C I = 1,N; J = 1,NPAR. +C +C IN ADDITION TO THE VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION +C WITH ODESSA_STESA USES THE FOLLOWING.. +C Y = AN NROW (=N) BY NCOL (=NSV) REAL ARRAY CONTAINING THE +C CORRECTED DEPENDENT VARIABLES ON OUTPUT.. +C Y(I,1) , I = 1,N = STATE VARIABLES (INPUT); +C Y(I,J) , I = 1,N , J = 2,NSV , +C = SENSITIVITY COEFFICIENTS, DY(I)/DP(J). +C YH = AN N BY NSV BY LMAX REAL ARRAY CONTAINING THE PREDICTED +C DEPENDENT VARIABLES AND THEIR APPROXIMATE SCALED DERIVATIVES. +C SAVF = A REAL ARRAY OF LENGTH N USED TO STORE FIRST DERIVATIVES +C OF DEPENDENT VARIABLES IF MITER = 2 OR 5. +C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING THE EQUATION +C PARAMETERS OF INTEREST. +C NRS = AN INTEGER ARRAY OF LENGTH NPAR + 1 CONTAINING THE NUMBER +C OF REPEATED STEPS (KFLAGS .LT. 0) DUE TO THE SENSITIVITY +C CALCULATIONS.. +C NRS(1) = TOTAL NUMBER OF REPEATED STEPS +C NRS(I) , I = 2,NPAR = NUMBER OF REPEATED STEPS DUE +C TO PARAMETER I. +C NSV = NUMBER OF SOLUTION VECTORS = NPAR + 1. +C KFLAGS = LOCAL ERROR TEST FLAG, = 0 IF TEST PASSES, .LT. 0 IF TEST +C FAILS, AND STEP NEEDS TO BE REPEATED. ERROR TEST IS APPLIED +C TO EACH SOLUTION VECTOR INDEPENDENTLY. +C DUPS, DSMS, DDNS = REAL SCALARS USED FOR COMPUTING RHUP, RHSM, RHDN, +C ON RETURN TO ODESSA_STODE (IALTH .EQ. 1). +C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, IALTH, LMAX, +C IERPJ, IERSL, JCUR, KFLAG, L, MITER, N, NQ, NFE, AND JOPT. +C----------------------------------------------------------------------- + DUPS = ZERO + DSMS = ZERO + DDNS = ZERO + HL0 = H*EL0 + EL0I = ONE/EL0 + TI2 = ONE/TESCO(2,NQ) + TI3 = ONE/TESCO(3,NQ) +C IF MITER = 2 OR 5 (OR IDF = 0), SUPPLY DERIVATIVES AT CORRECTED +C Y(*,1) VALUES FOR NUMERICAL DIFFERENTIATION IN PJAC AND/OR PDF. + IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. IDF .EQ. 0) GO TO 10 + GO TO 15 + 10 CALL F (NEQ, TN, Y, PAR, SAVF) + NFE = NFE + 1 +C IF JCUR = 0, UPDATE THE JACOBIAN MATRIX. +C IF MITER = 5, LOAD CORRECTED Y(*,1) VALUES INTO Y(*,2). + 15 IF (JCUR .EQ. 1) GO TO 30 + IF (MITER .NE. 5) GO TO 25 + DO 20 I = 1,N + 20 Y(I,2) = Y(I,1) + 25 CALL PJAC (NEQ, Y, Y(1,2), N, WM, IWM, EWT, SAVF, ACOR(1,2), + 1 PAR, F, JAC, JOPT) + IF (IERPJ .NE. 0) RETURN +C----------------------------------------------------------------------- +C THIS IS A LOOPING POINT FOR THE SENSITIVITY CALCULATIONS. +C----------------------------------------------------------------------- +C FOR EACH PARAMETER PAR(*), A SENSITIVITY SOLUTION VECTOR IS COMPUTED +C USING THE SAME STEP SIZE (H) AND ORDER (NQ) AS IN ODESSA_STODE. +C A LOCAL ERROR TEST IS APPLIED INDEPENDENTLY TO EACH SOLUTION VECTOR. +C----------------------------------------------------------------------- + 30 DO 100 J = 2,NSV + JPAR = J - 1 +C EVALUATE INHOMOGENEITY TERM, TEMPORARILY LOAD INTO Y(*,JPAR+1). ------ + CALL PDF(NEQ, Y, WM, SAVF, ACOR(1,J), Y(1,J), PAR, + 1 F, DF, JPAR) +C----------------------------------------------------------------------- +C LOAD RHS OF SENSITIVITY SOLUTION (CORRECTOR) EQUATION.. +C +C RHS = DY/DP - EL(1)*H*D(DY/DP)/DT + EL(1)*H*DF/DP +C +C----------------------------------------------------------------------- + DO 40 I = 1,N + 40 Y(I,J) = YH(I,J,1) - EL0*YH(I,J,2) + HL0*Y(I,J) +C----------------------------------------------------------------------- +C SOLVE CORRECTOR EQUATION: THE SOLUTIONS ARE LOCATED IN Y(*,JPAR+1). +C THE EXPLICIT FORMULA IS.. +C +C (I - EL(1)*H*JAC) * DY/DP(CORRECTED) = RHS +C +C----------------------------------------------------------------------- + CALL SOLVE (WM, IWM, Y(1,J), DUM) + IF (IERSL .NE. 0) RETURN +C ESTIMATE LOCAL TRUNCATION ERROR. ------------------------------------- + DO 50 I = 1,N + 50 ACOR(I,J) = (Y(I,J) - YH(I,J,1))*EL0I + ERR = ODESSA_VNORM(N, ACOR(1,J), EWT(1,J))*TI2 + IF (ERR .GT. ONE) GO TO 200 +C----------------------------------------------------------------------- +C LOCAL ERROR TEST PASSED. SET KFLAGS TO 0 TO INDICATE THIS. +C IF IALTH = 1, COMPUTE DSMS, DDNS, AND DUPS (IF L .LT. LMAX). +C----------------------------------------------------------------------- + KFLAGS = 0 + IF (IALTH .GT. 1) GO TO 100 + IF (L .EQ. LMAX) GO TO 70 + DO 60 I= 1,N + 60 Y(I,J) = ACOR(I,J) - YH(I,J,LMAX) + DUPS = MAX(DUPS,ODESSA_VNORM(N,Y(1,J),EWT(1,J))*TI3) + 70 DSMS = MAX(DSMS,ERR) + 100 CONTINUE + RETURN +C----------------------------------------------------------------------- +C THIS SECTION IS REACHED IF THE ERROR TOLERANCE FOR SENSITIVITY +C SOLUTION VECTOR JPAR HAS BEEN VIOLATED. KFLAGS IS MADE NEGATIVE TO +C INDICATE THIS. IF KFLAGS = -1, SET KFLAG EQUAL TO ZERO SO THAT KFLAG +C IS SET TO -1 ON RETURN TO ODESSA_STODE BEFORE REPEATING THE STEP. +C INCREMENT NRS(1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO ALL +C SENSITIVITY SOLUTION VECTORS) BY ONE. +C INCREMENT NRS(JPAR+1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO +C SOLUTION VECTOR JPAR+1) BY ONE. +C LOAD DSMS FOR RH CALCULATION IN ODESSA_STODE. +C----------------------------------------------------------------------- + 200 KFLAGS = KFLAGS - 1 + IF (KFLAGS .EQ. -1) KFLAG = 0 + NRS(1) = NRS(1) + 1 + NRS(J) = NRS(J) + 1 + DSMS = ERR + RETURN +C-------------------- END OF SUBROUTINE ODESSA_STESA ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_stode.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,518 @@ + SUBROUTINE ODESSA_STODE (NEQ, Y, YH, NYH, YH1, WM, IWM, EWT, SAVF, + 1 ACOR, PAR, NRS, F, JAC, DF, PJAC, PDF, SLVS) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + EXTERNAL F, JAC, DF, PJAC, PDF, SLVS + DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), WM(*), IWM(*), EWT(*), + 1 SAVF(*), ACOR(*), PAR(*), NRS(*) + PARAMETER (ONE=1.0D0,ZERO=0.0D0) + COMMON /ODE001/ ROWND, + 1 CONIT, CRATE, EL(13), ELCO(13,12), HOLD, RMAX, + 2 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, + 3 IOWND1(14), IPUP, MEO, NQNYH, NSLP, + 4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, + 5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU + COMMON /ODE002/ DUPS, DSMS, DDNS, + 1 IOWND2(3), ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS +C----------------------------------------------------------------------- +C ODESSA_STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE +C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. +C NOTE.. ODESSA_STODE 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 FOR ISOPT = 1, ODESSA_STODE CALLS ODESSA_STESA FOR SENSITIVITY CALCULATIONS. +C VARIABLES USED FOR COMMUNICATION WITH ODESSA_STESA ARE DESCRIBED IN +C ODESSA_STESA. COMMUNICATION WITH ODESSA_STODE IS DONE WITH THE +C FOLLOWING VARIABLES.. +C +C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND +C NUMBER OF PARAMETERS TO BE CONSIDERED IN THE SENSITIVITY +C ANALYSIS NEQ(2) (FOR ISOPT = 1), AND PASSED AS THE +C NEQ ARGUMENT IN ALL CALLS TO F, JAC, AND DF. +C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN +C ALL CALLS TO F, JAC, AND DF. +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 THE TOTAL NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.. +C NYH = N, ISOPT = 0, +C NYH = N * (NPAR + 1), ISOPT = 1 +C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH. +C EWT = AN ARRAY OF LENGTH NYH 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 NYH, 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 IF ISOPT = 1, PJAC CAN BE CALLED TO CALCULATE JAC BY +C SETTING JOPT = 1. +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, OR MITER. +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, (OR ODESSA_STESA). +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 (= 3, IF ISOPT = 0) +C (= 4, IF ISOPT = 1) +C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0). +C IF ISOPT = 1, PJAC IS CALLED AT LEAST ONCE EVERY STEP. +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 MODEL DIFFERENTIAL EQUATIONS. +C----------------------------------------------------------------------- + KFLAG = 0 + KFLAGS = 0 + TOLD = TN + NCF = 0 + IERPJ = 0 + IERSL = 0 + JCUR = 0 + ICF = 0 + 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 AT 2 +C FOR THE NEXT INCREASE. +C THESE COMPUTATIONS CONSIDER ONLY THE ORIGINAL SOLUTION VECTOR. +C THE SENSITIVITY SOLUTION VECTORS ARE CONSIDERED IN ODESSA_STESA (ISOPT = 1). +C----------------------------------------------------------------------- + LMAX = MAXORD + 1 + NQ = 1 + L = 2 + IALTH = 2 + RMAX = 10000.0D0 + RC = ZERO + EL0 = ONE + CRATE = 0.7D0 + DELP = ZERO + 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, ODESSA_CFODE 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 ODESSA_CFODE (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.5D0/REAL(NQ+2) + DDN = ODESSA_VNORM (N, SAVF, EWT)/TESCO(1,L) + EXDN = ONE/REAL(L) + RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0) + RH = MIN(RHDN,ONE) + IREDO = 3 + IF (H .EQ. HOLD) GO TO 170 + RH = MIN(RH,ABS(H/HOLD)) + H = HOLD + GO TO 175 +C----------------------------------------------------------------------- +C ODESSA_CFODE 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 ODESSA_CFODE (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.5D0/REAL(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(ONE,ABS(H)*HMXI*RH) + R = ONE + DO 180 J = 2,L + R = R*RH + DO 180 I = 1,NYH + 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 FOR ISOPT = 0, +C AND AT LEAST ONCE EVERY STEP FOR ISOPT = 1. +C----------------------------------------------------------------------- + 200 IF (ABS(RC-ONE) .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 + DO 210 I = I1,NQNYH + 210 YH1(I) = YH1(I) + YH1(I+NYH) + 215 CONTINUE +C----------------------------------------------------------------------- +C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. (= 3, FOR ISOPT = 0; +C = 4, FOR ISOPT = 1). A CONVERGENCE TEST IS MADE ON THE R.M.S. NORM +C OF EACH CORRECTION, WEIGHTED BY THE ERROR WEIGHT VECTOR EWT. THE SUM +C OF THE CORRECTIONS IS ACCUMULATED IN THE VECTOR ACOR(I), I = 1,N. +C (ACOR(I), I = N+1,NYH IS LOADED IN SUBROUTINE ODESSA_STESA (ISOPT = 1).) +C 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, PAR, 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----------------------------------------------------------------------- + IPUP = 0 + RC = ONE + NSLP = NST + CRATE = 0.7D0 + CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, ACOR, PAR, + 1 F, JAC, JOPT) + IF (IERPJ .NE. 0) GO TO 430 + 250 DO 260 I = 1,N + 260 ACOR(I) = ZERO + 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 (IF ISOPT = 1, FUNCTIONAL ITERATION IS NOT ALLOWED.) +C----------------------------------------------------------------------- + DO 290 I = 1,N + SAVF(I) = H*SAVF(I) - YH(I,2) + 290 Y(I) = SAVF(I) - ACOR(I) + DEL = ODESSA_VNORM (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 = ODESSA_VNORM (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.2D0*CRATE,DEL/DELP) + DCON = DEL*MIN(ONE,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) + IF (DCON .LE. ONE) GO TO 450 + M = M + 1 + IF (M .EQ. MAXCOR) GO TO 410 + IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410 + DELP = DEL + CALL F (NEQ, TN, Y, PAR, SAVF) + NFE = NFE + 1 + GO TO 270 +C----------------------------------------------------------------------- +C THE CORRECTOR ITERATION FAILED TO CONVERGE IN MAXCOR TRIES. +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.0D0 + TN = TOLD + I1 = NQNYH + 1 + DO 445 JB = 1,NQ + I1 = I1 - NYH + 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.00001D0) GO TO 670 + IF (NCF .EQ. MXNCF) GO TO 670 + RH = 0.25D0 + IPUP = MITER + IREDO = 1 + GO TO 170 +C----------------------------------------------------------------------- +C THE CORRECTOR HAS CONVERGED. +C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500 +C IF IT FAILS. OTHERWISE, ODESSA_STESA IS CALLED (ISOPT = 1) TO PERFORM +C SENSITIVITY CALCULATIONS AT CURRENT STEP SIZE AND ORDER. +C----------------------------------------------------------------------- + 450 CONTINUE + IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) + IF (M .GT. 0) DSM = ODESSA_VNORM (N, ACOR, EWT)/TESCO(2,NQ) + IF (DSM .GT. ONE) GO TO 500 +C + IF (ISOPT .EQ. 0) GO TO 460 +C----------------------------------------------------------------------- +C CALL ODESSA_STESA TO PERFORM EXPLICIT SENSITIVITY ANALYSIS. +C IF THE LOCAL ERROR TEST FAILS (WITHIN ODESSA_STESA) FOR ANY SOLUTION +C VECTOR, KFLAGS IS SET .LT. 0 AND CONTROL PASSES TO STATEMENT 500 UPON +C RETURN. IN EITHER CASE, JCUR IS SET TO ZERO TO SIGNAL THAT THE +C JACOBIAN MAY NEED UPDATING LATER. +C----------------------------------------------------------------------- + CALL ODESSA_STESA (NEQ, Y, N, NSV, YH, WM, IWM, EWT, SAVF, ACOR, + 1 PAR, NRS, F, JAC, DF, PJAC, PDF, SLVS) + IF (IERPJ .NE. 0 .OR. IERSL .NE. 0) GO TO 680 + IF (KFLAGS .LT. 0) 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----------------------------------------------------------------------- + 460 JCUR = 0 + KFLAG = 0 + IREDO = 0 + NST = NST + 1 + HU = H + NQU = NQ + DO 470 J = 1,L + DO 470 I = 1,NYH + 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,NYH + 490 YH(I,LMAX) = ACOR(I) + GO TO 700 +C----------------------------------------------------------------------- +C THE ERROR TEST FAILED IN EITHER ODESSA_STODE OR ODESSA_STESA. +C 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 + JCUR = 0 + TN = TOLD + I1 = NQNYH + 1 + DO 515 JB = 1,NQ + I1 = I1 - NYH + DO 510 I = I1,NQNYH + 510 YH1(I) = YH1(I) - YH1(I+NYH) + 515 CONTINUE + RMAX = 2.0D0 + IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660 + IF (KFLAG .LE. -3) GO TO 640 + IREDO = 2 + RHUP = ZERO + 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 FOR ISOPT = 1, DUPS AND DSMS ARE LOADED WITH THE LARGEST RMS-NORMS +C OBTAINED BY CONSIDERING SEPARATELY THE SENSITIVITY SOLUTION VECTORS. +C----------------------------------------------------------------------- + 520 RHUP = ZERO + IF (L .EQ. LMAX) GO TO 540 + DO 530 I = 1,N + 530 SAVF(I) = ACOR(I) - YH(I,LMAX) + DUP = ODESSA_VNORM (N, SAVF, EWT)/TESCO(3,NQ) + DUP = MAX(DUP,DUPS) + EXUP = ONE/REAL(L+1) + RHUP = ONE/(1.4D0*DUP**EXUP + 0.0000014D0) + 540 EXSM = ONE/REAL(L) + DSM = MAX(DSM,DSMS) + RHSM = ONE/(1.2D0*DSM**EXSM + 0.0000012D0) + RHDN = ZERO + IF (NQ .EQ. 1) GO TO 560 + JPOINT = 1 + DO 550 J = 1,NSV + DDN = ODESSA_VNORM (N, YH(JPOINT,L), EWT(JPOINT))/TESCO(1,NQ) + DDNS = MAX(DDNS,DDN) + JPOINT = JPOINT + N + 550 CONTINUE + DDN = DDNS + DDNS = ZERO + EXDN = ONE/REAL(NQ) + RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0) + 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. ONE) RH = ONE + GO TO 620 + 590 NEWQ = L + RH = RHUP + IF (RH .LT. 1.1D0) GO TO 610 + R = EL(L)/REAL(L) + DO 600 I = 1,NYH + 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.1D0)) GO TO 610 + IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) +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.1D0 + RH = MAX(HMIN/ABS(H),RH) + H = H*RH + DO 645 I = 1,NYH + 645 Y(I) = YH(I,1) + CALL F (NEQ, TN, Y, PAR, SAVF) + NFE = NFE + 1 + IF (ISOPT .EQ. 0) GO TO 649 + CALL ODESSA_SPRIME (NEQ, Y, YH, NYH, N, NSV, WM, IWM, EWT, SAVF, + 1 ACOR, ACOR(N+1), PAR, F, JAC, DF, PJAC, PDF) + IF (IERSP .LT. 0) GO TO 680 + DO 646 I = N+1,NYH + 646 YH(I,2) = H*YH(I,2) + 649 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.0D0 + 700 R = ONE/TESCO(2,NQU) + DO 710 I = 1,NYH + 710 ACOR(I) = ACOR(I)*R + 720 HOLD = H + JSTART = 1 + RETURN +C----------------------- END OF SUBROUTINE ODESSA_STODE ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_svcom.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,30 @@ + SUBROUTINE ODESSA_SVCOM (RSAV, ISAV) +C----------------------------------------------------------------------- +C THIS ROUTINE STORES IN RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS +C ODE001, ODE002 AND EH0001, WHICH ARE USED INTERNALLY IN THE ODESSA +C PACKAGE. +C RSAV = REAL ARRAY OF LENGTH 222 OR MORE. +C ISAV = INTEGER ARRAY OF LENGTH 52 OR MORE. +C----------------------------------------------------------------------- + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION RSAV(*), ISAV(*) + COMMON /ODE001/ RODE1(219), IODE1(39) + COMMON /ODE002/ RODE2(3), IODE2(11) + COMMON /EH0001/ IEH(2) + DATA LRODE1/219/, LIODE1/39/, LRODE2/3/, LIODE2/11/ +C + DO 10 I = 1,LRODE1 + 10 RSAV(I) = RODE1(I) + DO 20 I = 1,LRODE2 + J = LRODE1 + I + 20 RSAV(J) = RODE2(I) + DO 30 I = 1,LIODE1 + 30 ISAV(I) = IODE1(I) + DO 40 I = 1,LIODE2 + J = LIODE1 + I + 40 ISAV(J) = IODE2(I) + ISAV(LIODE1+LIODE2+1) = IEH(1) + ISAV(LIODE1+LIODE2+2) = IEH(2) + RETURN +C----------------------- END OF SUBROUTINE ODESSA_SVCOM ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_vnorm.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,72 @@ + DOUBLE PRECISION FUNCTION ODESSA_VNORM (N, V, W) +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 ODESSA_VNORM = SQRT( (1/N) * SUM( V(I)*W(I) )**2 ) +C PROTECTION FOR UNDERFLOW/OVERFLOW IS ACCOMPLISHED USING TWO +C CONSTANTS WHICH ARE HOPEFULLY APPLICABLE FOR ALL MACHINES. +C THESE ARE: +C CUTLO = maximum of SQRT(U/EPS) over all known machines +C CUTHI = minimum of SQRT(Z) over all known machines +C WHERE +C EPS = smallest number s.t. EPS + 1 .GT. 1 +C U = smallest positive number (underflow limit) +C Z = largest number (overflow limit) +C +C DETAILS OF THE ALGORITHM AND OF VALUES OF CUTLO AND CUTHI ARE +C FOUND IN THE BLAS ROUTINE SNRM2 (SEE ALSO ALGORITHM 539, TRANS. +C MATH. SOFTWARE, VOL. 5 NO. 3, 1979, 308-323. +C FOR SINGLE PRECISION, THE FOLLOWING VALUES SHOULD BE UNIVERSAL: +C DATA CUTLO,CUTHI /4.441E-16,1.304E19/ +C FOR DOUBLE PRECISION, USE +C DATA CUTLO,CUTHI /8.232D-11,1.304D19/ +C +C----------------------------------------------------------------------- + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + INTEGER NEXT,I,J,N + DIMENSION V(N),W(N) + DATA CUTLO,CUTHI /8.232D-11,1.304D19/ + DATA ZERO,ONE/0.0D0,1.0D0/ +C BLAS ALGORITHM + NEXT = 1 + SUM = ZERO + I = 1 +20 SX = V(I)*W(I) + GO TO (30,40,70,80),NEXT +30 IF (ABS(SX).GT.CUTLO) GO TO 110 + NEXT = 2 + XMAX = ZERO +40 IF (SX.EQ.ZERO) GO TO 130 + IF (ABS(SX).GT.CUTLO) GO TO 110 + NEXT = 3 + GO TO 60 +50 I=J + NEXT = 4 + SUM = (SUM/SX)/SX +60 XMAX = ABS(SX) + GO TO 90 +70 IF(ABS(SX).GT.CUTLO) GO TO 100 +80 IF(ABS(SX).LE.XMAX) GO TO 90 + SUM = ONE + SUM * (XMAX/SX)**2 + XMAX = ABS(SX) + GO TO 130 +90 SUM = SUM + (SX/XMAX)**2 + GO TO 130 +100 SUM = (SUM*XMAX)*XMAX +110 HITEST = CUTHI/REAL(N) + DO 120 J = I,N + SX = V(J)*W(J) + IF(ABS(SX).GE.HITEST) GO TO 50 + SUM = SUM + SX**2 +120 CONTINUE + ODESSA_VNORM = SQRT(SUM) + GO TO 140 +130 CONTINUE + I = I + 1 + IF (I.LE.N) GO TO 20 + ODESSA_VNORM = XMAX * SQRT(SUM) +140 CONTINUE + RETURN +C----------------------- END OF FUNCTION ODESSA_VNORM ------------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/odessa_xsetf.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,11 @@ + SUBROUTINE ODESSA_XSETF (MFLAG) +C +C THIS ROUTINE RESETS THE PRINT CONTROL FLAG MFLAG. +C + INTEGER MFLAG, MESFLG, LUNIT + COMMON /EH0001/ MESFLG, LUNIT +C + IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) MESFLG = MFLAG + RETURN +C----------------------- END OF SUBROUTINE ODESSA_XSETF ----------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/xerr.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,63 @@ + SUBROUTINE XERR (MSG, NERR, IERT, NI, I1, I2, NR, R1, R2) + INTEGER NERR, IERT, NI, I1, I2, NR, + 1 LUN, LUNIT, MESFLG + DOUBLE PRECISION R1, R2 + CHARACTER*(*) MSG +C------------------------------------------------------------------- +C +C ALL ARGUMENTS ARE INPUT ARGUMENTS. +C +C MSG = THE MESSAGE (CHARACTER VARIABLE) +C NERR = THE ERROR NUMBER (NOT USED). +C IERT = THE ERROR TYPE.. +C 1 MEANS RECOVERABLE (CONTROL RETURNS TO CALLER). +C 2 MEANS FATAL (RUN IS ABORTED--SEE NOTE BELOW). +C NI = NUMBER OF INTEGERS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE. +C I1,I2 = INTEGERS TO BE PRINTED, DEPENDING ON NI. +C NR = NUMBER OF REALS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE. +C R1,R2 = REALS TO BE PRINTED, DEPENDING ON NR. +C +C NOTES: +C 1. THE DIMENSION OF MSG IS ASSUMED TO BE AT MOST 60. +C (MULTI-LINE MESSAGES ARE GENERATED BY REPEATED CALLS.) +C 2. IF IERT = 2, CONTROL PASSES TO THE STATEMENT STOP +C TO ABORT THE RUN. THIS STATEMENT MAY BE MACHINE-DEPENDENT. +C 3. R1 AND R2 ARE ASSUMED TO BE IN DOUBLE PRECISION AND ARE PRINTED +C IN D21.13 FORMAT. +C 4. THE COMMON BLOCK /EH0001/ BELOW IS DATA-LOADED (A MACHINE- +C DEPENDENT FEATURE) WITH DEFAULT VALUES. +C THIS BLOCK IS NEEDED FOR PROPER RETENTION OF PARAMETERS USED BY +C THIS ROUTINE WHICH THE USER CAN RESET BY CALLING ODESSA_XSETF OR XSETUN. +C THE VARIABLES IN THIS BLOCK ARE AS FOLLOWS.. +C MESFLG = PRINT CONTROL FLAG.. +C 1 MEANS PRINT ALL MESSAGES (THE DEFAULT). +C 0 MEANS NO PRINTING. +C LUNIT = LOGICAL UNIT NUMBER FOR MESSAGES. +C THE DEFAULT IS 6 (MACHINE-DEPENDENT). +C 5. TO CHANGE THE DEFAULT OUTPUT UNIT, CHANGE THE DATA STATEMENT +C IN THE BLOCK DATA SUBPROGRAM BELOW. +C +C FOR A DIFFERENT RUN-ABORT COMMAND, CHANGE THE STATEMENT FOLLOWING +C STATEMENT 100 AT THE END. +C----------------------------------------------------------------------- + COMMON /EH0001/ MESFLG, LUNIT + IF (MESFLG .EQ. 0) GO TO 100 +C GET LOGICAL UNIT NUMBER. --------------------------------------------- + LUN = LUNIT +C WRITE THE MESSAGE. --------------------------------------------------- + WRITE (LUN, 10) MSG + 10 FORMAT(1X,A) +C----------------------------------------------------------------------- + IF (NI .EQ. 1) WRITE (LUN, 20) I1 + 20 FORMAT(6X,'IN ABOVE MESSAGE, I1 = ',I10) + IF (NI .EQ. 2) WRITE (LUN, 30) I1,I2 + 30 FORMAT(6X,'IN ABOVE MESSAGE, I1 = ',I10,3X,'I2 = ',I10) + IF (NR .EQ. 1) WRITE (LUN, 40) R1 + 40 FORMAT(6X,'IN ABOVE MESSAGE, R1 = ',D21.13) + IF (NR .EQ. 2) WRITE (LUN, 50) R1,R2 + 50 FORMAT(6X,'IN ABOVE, R1 = ',D21.13,3X,'R2 = ',D21.13) +C ABORT THE RUN IF IERT = 2. ------------------------------------------- + 100 IF (IERT .NE. 2) RETURN + STOP +C----------------------- END OF SUBROUTINE XERR ---------------------- + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/odessa/xsetun.f Thu Jul 11 02:37:44 2002 +0000 @@ -0,0 +1,11 @@ + SUBROUTINE XSETUN (LUN) +C +C THIS ROUTINE RESETS THE LOGICAL UNIT NUMBER FOR MESSAGES. +C + INTEGER LUN, MESFLG, LUNIT + COMMON /EH0001/ MESFLG, LUNIT +C + IF (LUN .GT. 0) LUNIT = LUN + RETURN +C----------------------- END OF SUBROUTINE XSETUN ---------------------- + END