Mercurial > hg > octave-max
changeset 5173:1278a2bc1527
[project @ 2005-03-02 01:33:37 by jwe]
| author | jwe |
|---|---|
| date | Wed, 02 Mar 2005 01:38:31 +0000 |
| parents | 9d9bbda4f00c |
| children | 4515ad040dc7 |
| files | ChangeLog configure.in doc/interpreter/diffeq.txi libcruft/ChangeLog libcruft/Makefile.in libcruft/odessa/Makefile.in libcruft/odessa/dodessa.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 liboctave/ChangeLog liboctave/Makefile.in liboctave/ODESSA-opts.in liboctave/ODESSA.cc liboctave/ODESSA.h src/ChangeLog src/Makefile.in |
| diffstat | 27 files changed, 27 insertions(+), 4336 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +2005-03-01 John W. Eaton <jwe@octave.org> + + * configure.in (AC_CONFIG_FILES): Remove libcruft/odessa/Makefile + from the list. + 2005-03-01 Todd Neal <tolchz@gmail.com> * examples/make_int.cc: DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA now
--- a/configure.in +++ b/configure.in @@ -29,7 +29,7 @@ EXTERN_CXXFLAGS="$CXXFLAGS" AC_INIT -AC_REVISION($Revision: 1.461 $) +AC_REVISION($Revision: 1.462 $) AC_PREREQ(2.57) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -1439,7 +1439,7 @@ libcruft/dasrt/Makefile libcruft/dassl/Makefile \ libcruft/fftpack/Makefile libcruft/lapack/Makefile \ libcruft/minpack/Makefile libcruft/misc/Makefile \ - libcruft/odepack/Makefile libcruft/odessa/Makefile \ + libcruft/odepack/Makefile \ libcruft/ordered-qz/Makefile libcruft/quadpack/Makefile \ libcruft/ranlib/Makefile libcruft/slatec-fn/Makefile \ libcruft/slatec-err/Makefile libcruft/villad/Makefile \
--- a/doc/interpreter/diffeq.txi +++ b/doc/interpreter/diffeq.txi @@ -87,14 +87,6 @@ @end group @end example -Octave also includes @sc{Odessa}, a modification of @sc{Lsode} that -allows for the computation of parametric sensitivity information -simultaneously with the solution of the set of ODEs. - -@DOCSTRING(odessa) - -@DOCSTRING(odessa_options) - See Alan C. Hindmarsh, @cite{ODEPACK, A Systematized Collection of ODE Solvers}, in Scientific Computing, R. S. Stepleman, editor, (1983) for more information about the inner workings of @code{lsode}.
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,5 +1,8 @@ 2005-03-01 John W. Eaton <jwe@octave.org> + * Makefile.in (CRUFT_DIRS): Remove it from the list. + * odessa: Delete directory. + * misc/machar.c (rmachar): Declare local REAL variables volatile. 2005-02-25 John W. Eaton <jwe@octave.org>
--- a/libcruft/Makefile.in +++ b/libcruft/Makefile.in @@ -30,7 +30,7 @@ CRUFT_DIRS = amos @BLAS_DIR@ blas-xtra daspk dasrt dassl \ @FFT_DIR@ @LAPACK_DIR@ lapack-xtra minpack \ - misc odepack odessa ordered-qz quadpack ranlib \ + misc odepack ordered-qz quadpack ranlib \ slatec-err slatec-fn villad SUBDIRS = $(CRUFT_DIRS)
deleted file mode 100644 --- a/libcruft/odessa/Makefile.in +++ /dev/null @@ -1,19 +0,0 @@ -# -# 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
deleted file mode 100644 --- a/libcruft/odessa/dodessa.f +++ /dev/null @@ -1,1910 +0,0 @@ -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 -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 -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 DGETRF AND DGETRS ARE ROUTINES FROM LAPACK FOR SOLVING FULL -C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS. -C DGBTRF AND DGBTRS ARE ROUTINES FROM LAPACK FOR SOLVING BANDED -C LINEAR SYSTEMS. -C 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 XERRWD, 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 DODESSA (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) = DSQRT(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) = DSQRT(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 = DABS(TOUT - T) - W0 = DMAX1(DABS(T),DABS(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 = DMAX1(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 = DABS(Y(I)) - IF (AYI .NE. ZERO) TOL = DMAX1(TOL,ATOLI/AYI) - 150 CONTINUE - 160 TOL = DMAX1(TOL,100.0D0*UROUND) - TOL = DMIN1(TOL,0.001D0) - SUM = ODESSA_VNORM (N, RWORK(LF0), RWORK(LEWT)) - SUM = ONE/(TOL*W0*W0) + TOL*SUM**2 - H0 = ONE/DSQRT(SUM) - H0 = MIN(H0,TDIST) - H0 = DSIGN(H0,TOUT-T) -C ADJUST H0 IF NECESSARY TO MEET HMAX BOUND. --------------------------- - 180 RH = DABS(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 = DABS(TN) + DABS(H) - IHIT = DABS(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 XERRWD ('ODESSA - WARNING..INTERNAL T (=R1) AND H (=R2) ARE', - 1 50, 101, 1, 0, 0, 0, 0, ZERO, ZERO) - CALL XERRWD - 1 ('SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP', - 1 52, 101, 1, 0, 0, 0, 0, ZERO, ZERO) - CALL XERRWD ('(H = STEP SIZE). SOLVER WILL CONTINUE ANYWAY', - 1 44, 101, 1, 0, 0, 0, 2, TN, H) - IF (NHNIL .LT. MXHNIL) GO TO 290 - CALL XERRWD ('ODESSA - ABOVE WARNING HAS BEEN ISSUED I1 TIMES.', - 1 48, 102, 1, 0, 0, 0, 0, ZERO, ZERO) - CALL XERRWD ('IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM', - 1 44, 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 = DABS(TN) + DABS(H) - IHIT = DABS(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 = DABS(TN) + DABS(H) - IHIT = DABS(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 XERRWD ('ODESSA -- REPEATED CALLS WITH ISTATE = 1 AND - 1 TOUT = T (=R1)', 59, 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 XERRWD ('ODESSA - AT CURRENT T (=R1), MXSTEP (=I1) STEPS', - 1 47, 201, 1, 0, 0, 0, 0, ZERO,ZERO) - CALL XERRWD ('TAKEN ON THIS CALL BEFORE REACHING TOUT', - 1 39, 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 XERRWD ('ODESSA - AT T (=R1), EWT(I1) HAS BECOME R2 .LE. 0.', - 1 50, 202, 1, 1, I, 0, 2, TN, EWTI) - ISTATE = -6 - GO TO 580 -C TOO MUCH ACCURACY REQUESTED FOR MACHINE PRECISION. ------------------- - 520 CALL XERRWD ('ODESSA - AT T (=R1), TOO MUCH ACCURACY REQUESTED', - 1 48, 203, 1, 0, 0, 0, 0, ZERO,ZERO) - CALL XERRWD ('FOR PRECISION OF MACHINE.. SEE TOLSF (=R2)', - 1 43, 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 XERRWD ('ODESSA - AT T(=R1) AND STEP SIZE H(=R2), THE ERROR', - 1 50, 204, 1, 0, 0, 0, 0, ZERO, ZERO) - CALL XERRWD ('TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN', - 1 44, 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 XERRWD ('ODESSA - AT T (=R1) AND STEP SIZE H (=R2), THE', - 1 46, 205, 1, 0, 0, 0, 0, ZERO,ZERO) - CALL XERRWD ('CORRECTOR CONVERGENCE FAILED REPEATEDLY', - 1 39, 205, 1, 0, 0, 0, 0, ZERO,ZERO) - CALL XERRWD ('OR WITH ABS(H) = HMIN', - 1 21, 0, 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 = DABS(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 XERRWD ('ODESSA - ISTATE (=I1) ILLEGAL', - 1 29, 1, 1, 1, ISTATE, 0, 0, ZERO,ZERO) - GO TO 700 - 602 CALL XERRWD ('ODESSA - ITASK (=I1) ILLEGAL', - 1 28, 2, 1, 1, ITASK, 0, 0, ZERO,ZERO) - GO TO 700 - 603 CALL XERRWD ('ODESSA - ISTATE .GT. 1 BUT ODESSA NOT INITIALIZED', - 1 49, 3, 1, 0, 0, 0, 0, ZERO,ZERO) - GO TO 700 - 604 CALL XERRWD ('ODESSA - NEQ (=I1) .LT. 1', - 1 25, 4, 1, 1, NEQ(1), 0, 0, ZERO,ZERO) - GO TO 700 - 605 CALL XERRWD ('ODESSA - ISTATE = 3 AND NEQ CHANGED. (I1 TO I2)', - 1 48, 5, 1, 2, N, NEQ(1), 0, ZERO,ZERO) - GO TO 700 - 606 CALL XERRWD ('ODESSA - ITOL (=I1) ILLEGAL', - 1 27, 6, 1, 1, ITOL, 0, 0, ZERO,ZERO) - GO TO 700 - 607 CALL XERRWD ('ODESSA - IOPT (=I1) ILLEGAL', - 1 27, 7, 1, 1, IOPT, 0, 0, ZERO,ZERO) - GO TO 700 - 608 CALL XERRWD('ODESSA - MF (=I1) ILLEGAL', - 1 25, 8, 1, 1, MF, 0, 0, ZERO,ZERO) - GO TO 700 - 609 CALL XERRWD('ODESSA - ML (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)', - 1 50, 9, 1, 2, ML, NEQ(1), 0, ZERO,ZERO) - GO TO 700 - 610 CALL XERRWD('ODESSA - MU (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)', - 1 50, 10, 1, 2, MU, NEQ(1), 0, ZERO,ZERO) - GO TO 700 - 611 CALL XERRWD('ODESSA - MAXORD (=I1) .LT. 0', - 1 28, 11, 1, 1, MAXORD, 0, 0, ZERO,ZERO) - GO TO 700 - 612 CALL XERRWD('ODESSA - MXSTEP (=I1) .LT. 0', - 1 28, 12, 1, 1, MXSTEP, 0, 0, ZERO,ZERO) - GO TO 700 - 613 CALL XERRWD('ODESSA - MXHNIL (=I1) .LT. 0', - 1 28, 13, 1, 1, MXHNIL, 0, 0, ZERO,ZERO) - GO TO 700 - 614 CALL XERRWD('ODESSA - TOUT (=R1) BEHIND T (=R2)', - 1 34, 14, 1, 0, 0, 0, 2, TOUT, T) - CALL XERRWD('INTEGRATION DIRECTION IS GIVEN BY H0 (=R1)', - 1 42, 14, 1, 0, 0, 0, 1, H0, ZERO) - GO TO 700 - 615 CALL XERRWD('ODESSA - HMAX (=R1) .LT. 0.0', - 1 28, 15, 1, 0, 0, 0, 1, HMAX, ZERO) - GO TO 700 - 616 CALL XERRWD('ODESSA - HMIN (=R1) .LT. 0.0', - 1 28, 16, 1, 0, 0, 0, 1, HMIN, ZERO) - GO TO 700 - 617 CALL XERRWD('ODESSA - RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS - 1 LRW (=I2)', 60, 17, 1, 2, LENRW, LRW, 0, ZERO,ZERO) - GO TO 700 - 618 CALL XERRWD('ODESSA - IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS - 1 LIW (=I2)', 60, 18, 1, 2, LENIW, LIW, 0, ZERO,ZERO) - GO TO 700 - 619 CALL XERRWD('ODESSA - RTOL(I1) IS R1 .LT. 0.0', - 1 32, 19, 1, 1, I, 0, 1, RTOLI, ZREO) - GO TO 700 - 620 CALL XERRWD('ODESSA - ATOL(I1) IS R1 .LT. 0.0', - 1 32, 20, 1, 1, I, 0, 1, ATOLI, ZERO) - GO TO 700 -* - 621 EWTI = RWORK(LEWT+I-1) - CALL XERRWD('ODESSA - EWT(I1) IS R1 .LE. 0.0', - 1 31, 21, 1, 1, I, 0, 1, EWTI, ZERO) - GO TO 700 - 622 CALL XERRWD('ODESSA - TOUT (=R1) TOO CLOSE TO T(=R2) TO START - 1 INTEGRATION', 60, 22, 1, 0, 0, 0, 2, TOUT, T) - GO TO 700 - 623 CALL XERRWD('ODESSA - ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU - 1 (= R2)', 58, 23, 1, 1, ITASK, 0, 2, TOUT, TP) - GO TO 700 - 624 CALL XERRWD('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR - 1 (=R2)', 57, 24, 1, 0, 0, 0, 2, TCRIT, TN) - GO TO 700 - 625 CALL XERRWD('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT - 1 (=R2)', 57, 25, 1, 0, 0, 0, 2, TCRIT, TOUT) - GO TO 700 - 626 CALL XERRWD('ODESSA - AT START OF PROBLEM, TOO MUCH ACCURACY', - 1 47, 26, 1, 0, 0, 0, 0, ZERO,ZERO) - CALL XERRWD('REQUESTED FOR PRECISION OF MACHINE. SEE TOLSF (=R1)', - 1 51, 26, 1, 0, 0, 0, 1, TOLSF, ZERO) - RWORK(14) = TOLSF - GO TO 700 - 627 CALL XERRWD - 1 ('ODESSA - TROUBLE FROM ODESSA_INTDY. ITASK = I1, TOUT = R1', - 1 57, 27, 1, 1, ITASK, 0, 1, TOUT, ZERO) - GO TO 700 -C ERROR STATEMENTS ASSOCIATED WITH SENSITIVITY ANALYSIS. - 628 CALL XERRWD('ODESSA - NPAR (=I1) .LT. 1', - 1 26, 28, 1, 1, NPAR, 0, 0, ZERO,ZERO) - GO TO 700 - 629 CALL XERRWD('ODESSA - ISTATE = 3 AND NPAR CHANGED (I1 TO I2)', - 1 47, 29, 1, 2, NP, NPAR, 0, ZERO,ZERO) - GO TO 700 - 630 CALL XERRWD('ODESSA - MITER (=I1) ILLEGAL', - 1 28, 30, 1, 1, MITER, 0, 0, ZERO,ZERO) - GO TO 700 - 631 CALL XERRWD('ODESSA - TROUBLE IN ODESSA_SPRIME (IERPJ)', - 1 41, 31, 1, 0, 0, 0, 0, ZERO,ZERO) - GO TO 700 - 632 CALL XERRWD('ODESSA - TROUBLE IN ODESSA_SPRIME (MITER)', - 1 41, 32, 1, 0, 0, 0, 0, ZERO,ZERO) - GO TO 700 - 633 CALL XERRWD('ODESSA - FATAL ERROR IN ODESSA_STODE (KFLAG = -3)', - 1 49, 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 XERRWD('ODESSA - REPEATED OCCURRENCES OF ILLEGAL INPUT', - 1 46, 302, 1, 0, 0, 0, 0, ZERO,ZERO) -C - 800 CALL XERRWD('ODESSA - RUN ABORTED.. APPARENT INFINITE LOOP', - 1 45, 303, 2, 0, 0, 0, 0, ZERO,ZERO) - RETURN - 801 CALL XERRWD('ODESSA - RUN ABORTED', - 1 20, 304, 2, 0, 0, 0, 0, ZERO,ZERO) - RETURN -C-------------------- END OF SUBROUTINE ODESSA ------------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_addx.f +++ /dev/null @@ -1,11 +0,0 @@ - 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
deleted file mode 100644 --- a/libcruft/odessa/odessa_cfode.f +++ /dev/null @@ -1,108 +0,0 @@ - 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/DBLE(NQ) - NQM1 = NQ - 1 - FNQM1 = DBLE(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)/DBLE(I) - 120 XPIN = XPIN + TSIGN*PC(I)/DBLE(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)/DBLE(I) - AGAMQ = RQFAC*XPIN - RAGQ = ONE/AGAMQ - TESCO(2,NQ) = RAGQ - IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/DBLE(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 = DBLE(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) = DBLE(NQP1)/ELCO(1,NQ) - TESCO(3,NQ) = DBLE(NQ+2)/ELCO(1,NQ) - RQ1FAC = RQ1FAC/FNQ - 230 CONTINUE - RETURN -C----------------------- END OF SUBROUTINE ODESSA_CFODE ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_ewset.f +++ /dev/null @@ -1,19 +0,0 @@ - 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*DABS(YCUR(I)) + ATOLI - 10 CONTINUE - RETURN -C----------------------- END OF SUBROUTINE ODESSA_EWSET ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_intdy.f +++ /dev/null @@ -1,74 +0,0 @@ - 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 = DBLE(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 = DBLE(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 XERRWD('ODESSA_INTDY-- K (=I1) ILLEGAL', - 1 31, 51, 1, 1, K, 0, 0, ZERO,ZERO) - IFLAG = -1 - RETURN - 90 CALL XERRWD ('ODESSA_INTDY-- T (=R1) ILLEGAL', - 1 31, 52, 1, 0, 0, 0, 1, T, ZERO) - CALL XERRWD('T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2)', - 1 48, 52, 1, 0, 0, 0, 2, TP, TN) - IFLAG = -2 - RETURN -C------------------ END OF SUBROUTINE ODESSA_INTDY ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_prepd.f +++ /dev/null @@ -1,63 +0,0 @@ - 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 = DMAX1(SRUR*DABS(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 - - - - - - - - -
deleted file mode 100644 --- a/libcruft/odessa/odessa_prepj.f +++ /dev/null @@ -1,166 +0,0 @@ - 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 DGETRF IF MITER = 1 OR 2, AND BY DGBTRF 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*DABS(H)*UROUND*DBLE(N)*FAC - IF (R0 .EQ. ZERO) R0 = ONE - SRUR = WM(1) - J1 = 2 - DO 230 J = 1,N - YJ = Y(J) - R = DMAX1(SRUR*DABS(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 DGETRF ( N, N, WM(3), N, IWM(21), IER) - IF (IER .NE. 0) IERPJ = 1 - RETURN -C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. --------- - 300 WM(2) = HL0 - R = EL0*0.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 (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320 - IF (DABS(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 = MIN0(MBAND,N) - MEBAND = MBAND + ML - MEB1 = MEBAND - 1 - SRUR = WM(1) - FAC = ODESSA_VNORM (N, SAVF, EWT) - R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC - IF (R0 .EQ. ZERO) R0 = ONE - DO 560 J = 1,MBA - DO 530 I = J,N,MBAND - YI = Y(I) - R = DMAX1(SRUR*DABS(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 = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ)) - FAC = -HL0/R - I1 = MAX0(JJ-MU,1) - I2 = MIN0(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 DGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER) - IF (IER .NE. 0) IERPJ = 1 - RETURN -C---------------- END OF SUBROUTINE ODESSA_PREPJ ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_rscom.f +++ /dev/null @@ -1,26 +0,0 @@ - SUBROUTINE ODESSA_RSCOM (RSAV, ISAV) -C----------------------------------------------------------------------- -C THIS ROUTINE RESTORES FROM RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS -C ODE001 AND ODE002, 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) - 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,LIODE2 - J = LIODE1 + I - 40 IODE2(I) = ISAV(J) - RETURN -C----------------------- END OF SUBROUTINE ODESSA_RSCOM ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_solsy.f +++ /dev/null @@ -1,61 +0,0 @@ - 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 DGETRS 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 DGBTRS. -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 DGETRS ( 'N', N, 1, WM(3), N, IWM(21), X, N, INLPCK) - RETURN -C - 300 PHL0 = WM(2) - HL0 = H*EL0 - WM(2) = HL0 - IF (HL0 .EQ. PHL0) GO TO 330 - R = HL0/PHL0 - DO 320 I = 1,N - DI = ONE - R*(ONE - ONE/WM(I+2)) - IF (DABS(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 DGBTRS ( 'N', N, ML, MU, 1, WM(3), MEBAND, IWM(21), X, N, - * INLPCK) - RETURN -C----------------------- END OF SUBROUTINE ODESSA_SOLSY ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_sprim.f +++ /dev/null @@ -1,125 +0,0 @@ - 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
deleted file mode 100644 --- a/libcruft/odessa/odessa_stesa.f +++ /dev/null @@ -1,135 +0,0 @@ - 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 = DMAX1(DUPS,ODESSA_VNORM(N,Y(1,J),EWT(1,J))*TI3) - 70 DSMS = DMAX1(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
deleted file mode 100644 --- a/libcruft/odessa/odessa_stode.f +++ /dev/null @@ -1,518 +0,0 @@ - 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/DBLE(NQ+2) - DDN = ODESSA_VNORM (N, SAVF, EWT)/TESCO(1,L) - EXDN = ONE/DBLE(L) - RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0) - RH = DMIN1(RHDN,ONE) - IREDO = 3 - IF (H .EQ. HOLD) GO TO 170 - RH = DMIN1(RH,DABS(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/DBLE(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 = DMAX1(RH,HMIN/DABS(H)) - 175 RH = DMIN1(RH,RMAX) - RH = RH/DMAX1(ONE,DABS(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 (DABS(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 = DMAX1(0.2D0*CRATE,DEL/DELP) - DCON = DEL*DMIN1(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 (DABS(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 (DABS(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 = DMAX1(DUP,DUPS) - EXUP = ONE/DBLE(L+1) - RHUP = ONE/(1.4D0*DUP**EXUP + 0.0000014D0) - 540 EXSM = ONE/DBLE(L) - DSM = DMAX1(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 = DMAX1(DDNS,DDN) - JPOINT = JPOINT + N - 550 CONTINUE - DDN = DDNS - DDNS = ZERO - EXDN = ONE/DBLE(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)/DBLE(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 = DMIN1(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 = DMAX1(HMIN/DABS(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
deleted file mode 100644 --- a/libcruft/odessa/odessa_svcom.f +++ /dev/null @@ -1,27 +0,0 @@ - SUBROUTINE ODESSA_SVCOM (RSAV, ISAV) -C----------------------------------------------------------------------- -C THIS ROUTINE STORES IN RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS -C ODE001 AND ODE002, 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) - 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) - RETURN -C----------------------- END OF SUBROUTINE ODESSA_SVCOM ----------------------- - END
deleted file mode 100644 --- a/libcruft/odessa/odessa_vnorm.f +++ /dev/null @@ -1,72 +0,0 @@ - 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 (DABS(SX).GT.CUTLO) GO TO 110 - NEXT = 2 - XMAX = ZERO -40 IF (SX.EQ.ZERO) GO TO 130 - IF (DABS(SX).GT.CUTLO) GO TO 110 - NEXT = 3 - GO TO 60 -50 I=J - NEXT = 4 - SUM = (SUM/SX)/SX -60 XMAX = DABS(SX) - GO TO 90 -70 IF(DABS(SX).GT.CUTLO) GO TO 100 -80 IF(DABS(SX).LE.XMAX) GO TO 90 - SUM = ONE + SUM * (XMAX/SX)**2 - XMAX = DABS(SX) - GO TO 130 -90 SUM = SUM + (SX/XMAX)**2 - GO TO 130 -100 SUM = (SUM*XMAX)*XMAX -110 HITEST = CUTHI/DBLE(N) - DO 120 J = I,N - SX = V(J)*W(J) - IF(DABS(SX).GE.HITEST) GO TO 50 - SUM = SUM + SX**2 -120 CONTINUE - ODESSA_VNORM = DSQRT(SUM) - GO TO 140 -130 CONTINUE - I = I + 1 - IF (I.LE.N) GO TO 20 - ODESSA_VNORM = XMAX * DSQRT(SUM) -140 CONTINUE - RETURN -C----------------------- END OF FUNCTION ODESSA_VNORM ------------------------- - END
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2005-03-01 John W. Eaton <jwe@octave.org> + + * ODESSA.h, ODESSA.cc, ODESSA-opts.in: Delete. + * Makefile.in: Remove them from the lists. + 2005-02-28 John W. Eaton <jwe@octave.org> * Makefile.in (LINK_DEPS): Remove -lglob from the list.
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -51,7 +51,7 @@ SPARSE_MX_OP_INC := $(shell $(AWK) -f $(srcdir)/sparse-mk-ops.awk prefix=smx list_h_files=1 $(srcdir)/sparse-mx-ops) OPTS_INC_DATA := DASPK-opts.in DASRT-opts.in DASSL-opts.in \ - LSODE-opts.in NLEqn-opts.in ODESSA-opts.in Quad-opts.in + LSODE-opts.in NLEqn-opts.in Quad-opts.in OPTS_INC := $(OPTS_INC_DATA:.in=.h) @@ -59,7 +59,7 @@ DAERTFunc.h DASPK.h DASRT.h DASSL.h FEGrid.h \ LinConst.h LP.h LPsolve.h LSODE.h NLConst.h NLEqn.h \ NLFunc.h NLP.h ODE.h ODEFunc.h ODES.h ODESFunc.h \ - ODESSA.h Objective.h QP.h Quad.h Range.h base-dae.h \ + Objective.h QP.h Quad.h Range.h base-dae.h \ base-de.h base-min.h byte-swap.h cmd-edit.h cmd-hist.h \ data-conv.h dir-ops.h file-ops.h file-stat.h getopt.h \ glob-match.h idx-vector.h kpse-xfns.h \ @@ -109,7 +109,7 @@ LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc \ DASPK.cc DASRT.cc DASSL.cc FEGrid.cc LinConst.cc \ - LPsolve.cc LSODE.cc NLEqn.cc ODES.cc ODESSA.cc \ + LPsolve.cc LSODE.cc NLEqn.cc ODES.cc \ Quad.cc Range.cc data-conv.cc dir-ops.cc \ file-ops.cc file-stat.cc glob-match.cc idx-vector.cc \ lo-ieee.cc lo-mappers.cc lo-specfun.cc lo-sysdep.cc \
deleted file mode 100644 --- a/liboctave/ODESSA-opts.in +++ /dev/null @@ -1,132 +0,0 @@ -CLASS = "ODESSA" - -INCLUDE = "ODES.h" - -OPTION - NAME = "absolute tolerance" - DOC_ITEM -Absolute tolerance. May be either vector or scalar. If a vector, it -must match the dimension of the state vector. - END_DOC_ITEM - TYPE = "Array<double>" - SET_ARG_TYPE = "const $TYPE&" - INIT_BODY - $OPTVAR.resize (1); - $OPTVAR(0) = ::sqrt (DBL_EPSILON); - END_INIT_BODY - SET_CODE - void set_$OPT (double val) - { - $OPTVAR.resize (1); - $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); - reset = true; - } - - void set_$OPT (const $TYPE& val) - { $OPTVAR = val; reset = true; } - END_SET_CODE -END_OPTION - -OPTION - NAME = "relative tolerance" - DOC_ITEM -Relative tolerance parameter. Unlike the absolute tolerance, this -parameter may only be a scalar. - -The local error test applied at each integration step is - -@example - abs (local error in x(i)) - <= rtol * abs (y(i)) + atol(i) -@end example - END_DOC_ITEM - TYPE = "double" - INIT_VALUE = "::sqrt (DBL_EPSILON)" - SET_EXPR = "(val > 0.0) ? val : ::sqrt (DBL_EPSILON)" -END_OPTION - -OPTION - NAME = "integration method" - DOC_ITEM -A string specifing the method of integration to use to solve the ODE -system. Valid values are - -@table @asis -@item \"adams\" -@itemx \"non-stiff\" -No Jacobian used (even if it is available). -@item \"bdf\" -@item \"stiff\" -Use stiff backward differentiation formula (BDF) method. If a -function to compute the Jacobian is not supplied, @code{lsode} will -compute a finite difference approximation of the Jacobian matrix. -@end table - END_DOC_ITEM - TYPE = "std::string" - SET_ARG_TYPE = "const $TYPE&" - INIT_VALUE = {"stiff"} - SET_BODY - if (val == "stiff" || val == "bdf") - $OPTVAR = "stiff"; - else if (val == "non-stiff" || val == "adams") - $OPTVAR = "non-stiff"; - else - (*current_liboctave_error_handler) - ("lsode_options: method must be \"stiff\", \"bdf\", \"non-stiff\", or \"adams\""); - END_SET_BODY -END_OPTION - -OPTION - NAME = "initial step size" - DOC_ITEM -The step size to be attempted on the first step (default is determined -automatically). - END_DOC_ITEM - TYPE = "double" - INIT_VALUE = "-1.0" - SET_EXPR = "(val >= 0.0) ? val : -1.0" -END_OPTION - -OPTION - NAME = "maximum order" - DOC_ITEM -Restrict the maximum order of the solution method. If using the Adams -method, this option must be between 1 and 12. Otherwise, it must be -between 1 and 5, inclusive. - END_DOC_ITEM - TYPE = "int" - INIT_VALUE = "-1" - SET_EXPR = "val" -END_OPTION - -OPTION - NAME = "maximum step size" - DOC_ITEM -Setting the maximum stepsize will avoid passing over very large -regions (default is not specified). - END_DOC_ITEM - TYPE = "double" - INIT_VALUE = "-1.0" - SET_EXPR = "(val >= 0.0) ? val : -1.0" -END_OPTION - -OPTION - NAME = "minimum step size" - DOC_ITEM -The minimum absolute step size allowed (default is 0). - END_DOC_ITEM - TYPE = "double" - INIT_VALUE = "0.0" - SET_EXPR = "(val >= 0.0) ? val : 0.0" -END_OPTION - -OPTION - NAME = "step limit" - DOC_ITEM -Maximum number of steps allowed (default is 100000). - END_DOC_ITEM - TYPE = "int" - INIT_VALUE = "100000" - SET_EXPR = "val" -END_OPTION -
deleted file mode 100644 --- a/liboctave/ODESSA.cc +++ /dev/null @@ -1,720 +0,0 @@ -/* - -Copyright (C) 2002 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cfloat> -#include <cmath> - -// For instantiating the Array<Matrix> object. -#include "Array.h" -#include "Array.cc" - -#include "ODESSA.h" -#include "f77-fcn.h" -#include "lo-error.h" -#include "lo-sstream.h" -#include "quit.h" - -typedef int (*odessa_fcn_ptr) (int*, const double&, double*, - double*, double*); - -typedef int (*odessa_jac_ptr) (int*, const double&, double*, - double*, const int&, const int&, - double*, const int&); - -typedef int (*odessa_dfdp_ptr) (int*, const double&, double*, - double*, double*, const int&); - - -extern "C" -{ - F77_RET_T - F77_FUNC (dodessa, DODESSA) (odessa_fcn_ptr, odessa_dfdp_ptr, int*, - double*, double*, double&, double&, - int&, double&, const double*, int&, - int&, int*, double*, int&, int*, int&, - odessa_jac_ptr, int&); -} - -INSTANTIATE_ARRAY (Matrix); - -static ODESFunc::ODES_fsub user_fsub; -static ODESFunc::ODES_bsub user_bsub; -static ODESFunc::ODES_jsub user_jsub; - - -static int -odessa_f (int* neq, const double& t, double *state, - double *par, double *fval) -{ - BEGIN_INTERRUPT_WITH_EXCEPTIONS; - - int n = neq[0]; - int n_par = neq[1]; - - // Load the state and parameter arrays as Octave objects - - ColumnVector tmp_state (n); - for (int i = 0; i < n; i++) - tmp_state(i) = state[i]; - - ColumnVector tmp_param (n_par); - for (int i = 0; i < n_par; i++) - tmp_param(i) = par[i]; - - ColumnVector tmp_fval = user_fsub (tmp_state, t, tmp_param); - - if (tmp_fval.length () == 0) - octave_jump_to_enclosing_context (); - else - { - for (int i = 0; i < n; i++) - fval[i] = tmp_fval(i); - } - - END_INTERRUPT_WITH_EXCEPTIONS; - - return 0; -} - -static int -odessa_j (int* neq, const double& t, double *state, - double *par, const int& /* ml */, const int& /* mu */, - double *pd, const int& nrowpd) -{ - BEGIN_INTERRUPT_WITH_EXCEPTIONS; - - int n = neq[0]; - int n_par = neq[1]; - - // Load the state and parameter arrays as Octave objects - ColumnVector tmp_state (n); - for (int i = 0; i < n; i++) - tmp_state(i) = state[i]; - - ColumnVector tmp_param (n_par); - for (int i = 0; i < n_par; i++) - tmp_param(i) = par[i]; - - Matrix tmp_fval = user_jsub (tmp_state, t, tmp_param); - - if (tmp_fval.length () == 0) - octave_jump_to_enclosing_context (); - else - { - for (int j = 0; j < n; j++) - for (int i = 0; i < nrowpd; i++) - pd[nrowpd*j+i] = tmp_fval(i,j); - } - - END_INTERRUPT_WITH_EXCEPTIONS; - - return 0; -} - -static int -odessa_b (int* neq, const double& t, double *state, - double *par, double *dfdp, const int& jpar) - -{ - BEGIN_INTERRUPT_WITH_EXCEPTIONS; - - int n = neq[0]; - int n_par = neq[1]; - - // Load the state and parameter arrays as Octave objects - ColumnVector tmp_state (n); - for (int i = 0; i < n; i++) - tmp_state(i) = state[i]; - - ColumnVector tmp_param (n_par); - for (int i = 0; i < n_par; i++) - tmp_param(i) = par[i]; - - ColumnVector tmp_fval = user_bsub (tmp_state, t, tmp_param, jpar); - - if (tmp_fval.length () == 0) - octave_jump_to_enclosing_context (); - else - { - for (int i = 0; i < n; i++) - dfdp[i] = tmp_fval(i); - } - - END_INTERRUPT_WITH_EXCEPTIONS; - - return 0; -} - -ODESSA::ODESSA (void) : ODES (), ODESSA_options () -{ - initialized = false; - - neq.resize(2); - n = size (); - - iopt.resize(4); - - itask = 1; - iopt(0) = 0; - isopt = 0; - iopt(1) = isopt; - npar = 0; - neq(0) = n; - neq(1) = npar; - - sanity_checked = false; -} - -ODESSA::ODESSA (const ColumnVector& s, double tm, ODESFunc& f) - : ODES (s, tm, f), ODESSA_options () -{ - initialized = false; - - neq.resize(2); - n = size (); - - iopt.resize(4); - itask = 1; - iopt(0) = 0; - isopt = 0; - iopt(1) = isopt; - - sanity_checked = false; - - npar = 0; - neq(0) = n; - neq(1) = npar; - - y.resize (n, 1, 0.0); -} - -ODESSA::ODESSA (const ColumnVector& s, const ColumnVector& xtheta, - const Matrix& sensitivity_guess, double tm, ODESFunc& f) - : ODES (s, xtheta, tm, f) -{ - initialized = false; - - neq.resize(2); - n = s.length(); - npar = xtheta.length(); - - neq(0) = n; - neq(1) = npar; - - sx0 = sensitivity_guess; - par.resize (npar); - - for (int i = 0; i < npar; i++) - { - par(i) = xtheta(i); - } - - sanity_checked = false; - - npar = xtheta.length (); - - iopt.resize(4); - itask = 1; - iopt(0) = 0; - isopt = 1; - iopt(1) = isopt; - - y.resize (n, npar+1, 0.0); -} - -void -ODESSA::integrate (double tout) -{ - ODESSA_result retval; - - if (! initialized) - { - - for (int i = 0; i < n; i++) - y(i,0) = x(i); - - if (npar > 0) - { - for (int j = 0; j < npar; j++) - for (int i = 0; i < n; i++) - y(i,j+1) = sx0(i,j); - } - - integration_error = false; - - user_fsub = ODESFunc::fsub_function (); - user_bsub = ODESFunc::bsub_function (); - user_jsub = ODESFunc::jsub_function (); - - int idf; - - if (user_bsub) - idf = 1; - else - idf = 0; - - iopt(2) = idf; - - if (restart) - { - restart = false; - istate = 1; - } - - int max_maxord = 0; - - if (integration_method () == "stiff") - { - if (user_jsub) - method_flag = 21; - else - method_flag = 22; - - max_maxord = 5; - - if (isopt) - { - liw = 21 + n + npar; - lrw = 22 + 8*(npar+1)*n + n*n + n; - } - else - { - liw = 20 + n; - lrw = 22 + 9*n + n*n; - } - } - else - { - max_maxord = 12; - - if (isopt) - { - if (user_jsub) - method_flag = 11; - else - method_flag = 12; - liw = 21 + n + npar; - lrw = 22 + 15*(npar+1)*n + n*n + n; - } - else - { - method_flag = 10; - liw = 20 + n; - lrw = 22 + 16 * n; - } - } - - if (iwork.length () != liw) - { - iwork.resize (liw); - - for (int i = 4; i < 10; i++) - iwork.elem (i) = 0; - } - - if (rwork.length () != lrw) - { - rwork.resize (lrw); - - for (int i = 4; i < 10; i++) - rwork.elem (i) = 0.0; - } - - maxord = maximum_order (); - - if (maxord >= 0) - { - if (maxord > 0 && maxord <= max_maxord) - { - iwork(4) = maxord; - iopt(0) = 1; - } - else - { - (*current_liboctave_error_handler) - ("odessa: invalid value for maximum order"); - integration_error = true; - return; - } - } - - initialized = true; - } - - integration_error = false; - - // NOTE: this won't work if LSODE passes copies of the state vector. - // In that case we have to create a temporary vector object - // and copy. - - - if (! sanity_checked) - { - ColumnVector fval = user_fsub (x, t, theta); - - if (fval.length () != x.length ()) - { - (*current_liboctave_error_handler) - ("odessa: inconsistent sizes for state and residual vectors"); - - integration_error = true; - return; - } - - sanity_checked = true; - } - - if (stop_time_set) - { - itask = 4; - rwork.elem (0) = stop_time; - iopt(0) = 1; - } - else - { - itask = 1; - } - - double rel_tol = relative_tolerance (); - - Array<double> abs_tol = absolute_tolerance (); //note; this should - // be a matrix, not a vector - - int abs_tol_len = abs_tol.length (); - - int itol; - - if (abs_tol_len == 1) - itol = 1; - else if (abs_tol_len == n) - itol = 2; - else - { - (*current_liboctave_error_handler) - ("odessa: inconsistent sizes for state and absolute tolerance vectors"); - - integration_error = 1; - return; - } - - if (initial_step_size () >= 0.0) - { - rwork.elem (4) = initial_step_size (); - iopt(0) = 1; - } - - if (maximum_step_size () >= 0.0) - { - rwork.elem (5) = maximum_step_size (); - iopt(0) = 1; - } - - if (minimum_step_size () >= 0.0) - { - rwork.elem (6) = minimum_step_size (); - iopt(0) = 1; - } - - if (step_limit () > 0) - { - iwork.elem (5) = step_limit (); - iopt(0) = 1; - } - - py = y.fortran_vec (); - piwork = iwork.fortran_vec (); - prwork = rwork.fortran_vec (); - ppar = par.fortran_vec (); - piopt = iopt.fortran_vec (); - pneq = neq.fortran_vec (); - - const double *pabs_tol = abs_tol.fortran_vec (); - - F77_XFCN (dodessa, DODESSA, (odessa_f, odessa_b, pneq, py, ppar, t, - tout, itol, rel_tol, pabs_tol, itask, - istate, piopt, prwork, lrw, piwork, liw, - odessa_j, method_flag)); - - if (f77_exception_encountered) - { - integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in odessa"); - } - else - { - switch (istate) - { - case 1: // prior to initial integration step. - case 2: // odessa was successful. - t = tout; - break; - - case -1: // excess work done on this call (perhaps wrong mf). - case -2: // excess accuracy requested (tolerances too small). - case -3: // illegal input detected (see printed message). - case -4: // repeated error test failures (check all inputs). - case -5: // repeated convergence failures (perhaps bad jacobian - // supplied or wrong choice of mf or tolerances). - case -6: // error weight became zero during problem. (solution - // component i vanished, and atol or atol(i) = 0.) - case -13: // Return requested in user-supplied function. - integration_error = 1; - break; - - default: - integration_error = 1; - (*current_liboctave_error_handler) - ("unrecognized value of istate (= %d) returned from odessa", - istate); - break; - } - } -} - -std::string -ODESSA::error_message (void) const -{ - std::string retval; - - OSSTREAM buf; - buf << t << OSSTREAM_ENDS; - std::string t_curr = OSSTREAM_STR (buf); - OSSTREAM_FREEZE (buf); - - switch (istate) - { - case 1: - retval = "prior to initial integration step"; - break; - - case 2: - retval = "successful exit"; - break; - - case 3: - retval = "prior to continuation call with modified parameters"; - break; - - case -1: - retval = std::string ("excess work on this call (t = ") - + t_curr + "; perhaps wrong integration method)"; - break; - - case -2: - retval = "excess accuracy requested (tolerances too small)"; - break; - - case -3: - retval = "invalid input detected (see printed message)"; - break; - - case -4: - retval = std::string ("repeated error test failures (t = ") - + t_curr + "check all inputs)"; - break; - - case -5: - retval = std::string ("repeated convergence failures (t = ") - + t_curr - + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)"; - break; - - case -6: - retval = std::string ("error weight became zero during problem. (t = ") - + t_curr - + "; solution component i vanished, and atol or atol(i) == 0)"; - break; - - case -13: - retval = "return requested in user-supplied function (t = " - + t_curr + ")"; - break; - - default: - retval = "unknown error state"; - break; - } - - return retval; -} - - -ODESSA_result -ODESSA::integrate (const ColumnVector& tout) -{ - ODESSA_result retval; - - Matrix x_out; - - Array<Matrix> x_s_out; - - int n_out = tout.capacity (); - - if (n_out > 0 && n > 0) - { - x_out.resize (n_out, n); - - x_s_out.resize (npar, Matrix (n_out, n, 0.0)); - - for (int j = 0; j < n_out; j++) - { - integrate (tout(j)); - - if (integration_error) - { - retval = ODESSA_result (x_out, x_s_out); - return retval; - } - - for (int i = 0; i < n; i++) - { - x_out(j,i) = y(i,0); - - for (int k = 0; k < npar; k++) - { - x_s_out(k)(j,i) = y(i,k+1); - } - } - } - } - - retval = ODESSA_result (x_out, x_s_out); - - return retval; -} - -ODESSA_result -ODESSA::integrate (const ColumnVector& tout, const ColumnVector& tcrit) -{ - ODESSA_result retval; - - Matrix x_out; - - Array<Matrix> x_s_out; - - int n_out = tout.capacity (); - - if (n_out > 0 && n > 0) - { - x_out.resize (n_out, n); - - x_s_out.resize (npar, Matrix (n_out, n, 0.0)); - - int n_crit = tcrit.capacity (); - - if (n_crit > 0) - { - int i_crit = 0; - int i_out = 0; - double next_crit = tcrit(0); - double next_out; - while (i_out < n_out) - { - bool do_restart = false; - - next_out = tout(i_out); - if (i_crit < n_crit) - next_crit = tcrit(i_crit); - - int save_output; - double t_out; - - if (next_crit == next_out) - { - set_stop_time (next_crit); - t_out = next_out; - save_output = true; - i_out++; - i_crit++; - do_restart = true; - } - else if (next_crit < next_out) - { - if (i_crit < n_crit) - { - set_stop_time (next_crit); - t_out = next_crit; - save_output = false; - i_crit++; - do_restart = true; - } - else - { - clear_stop_time (); - t_out = next_out; - save_output = true; - i_out++; - } - } - else - { - set_stop_time (next_crit); - t_out = next_out; - save_output = true; - i_out++; - } - integrate (t_out); - if (integration_error) - { - retval = ODESSA_result (x_out, x_s_out); - return retval; - } - if (save_output) - { - for (int i = 0; i < n; i++) - { - x_out(i_out-1,i) = y(i,0); - - for (int k = 0; k < npar; k++) - { - x_s_out(k)(i_out-1,i) = y(i,k+1); - } - } - } - - if (do_restart) - force_restart (); - } - - retval = ODESSA_result (x_out, x_s_out); - } - else - { - retval = integrate (tout); - - if (integration_error) - return retval; - } - } - - return retval; -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/liboctave/ODESSA.h +++ /dev/null @@ -1,134 +0,0 @@ -/* - -Copyright (C) 2002 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -#if !defined (octave_ODESSA_h) -#define octave_ODESSA_h 1 - -#include <cfloat> -#include <cmath> - -#include "ODESSA-opts.h" - -class -ODESSA_result -{ -public: - - ODESSA_result (void) { } - - ODESSA_result (const Matrix& xx, - const Array<Matrix>& xx_s) - - : x (xx), x_s (xx_s) { } - - ODESSA_result (const ODESSA_result& r) - : x (r.x), x_s (r.x_s) { } - - ODESSA_result& operator = (const ODESSA_result& r) - { - if (this != &r) - { - x = r.x; - x_s = r.x_s; - } - return *this; - } - - ~ODESSA_result (void) { } - - Matrix state (void) const { return x; } - Array<Matrix> state_sensitivity (void) const { return x_s; } - -private: - - Matrix x; - Array<Matrix> x_s; -}; - -class -ODESSA : public ODES, public ODESSA_options -{ -public: - - ODESSA (void); - - ODESSA (const ColumnVector& x, double time, ODESFunc& f); - - ODESSA (const ColumnVector& x, const ColumnVector& theta, - const Matrix& sensitivity_guess, double time, ODESFunc& f); - - ~ODESSA (void) { } - - ODESSA_result integrate (const ColumnVector& tout); - - ODESSA_result integrate (const ColumnVector& tout, - const ColumnVector& tcrit); - - std::string error_message (void) const; - -private: - - bool initialized; - - bool sanity_checked; - - int liw; - int lrw; - int method_flag; - int maxord; - Array<int> iwork; - Array<double> rwork; - int itask; - Array<int> iopt; - int isopt; - - Array<int> neq; - - int n; - int npar; - - // XXX FIXME XXX -- ??? - Array<double> par; - - Matrix sx0; - - Matrix y; - - double *py; - double *ppar; - int *piwork; - int *piopt; - int *pneq; - double *prwork; - - void init_work_size (int); - - void integrate (double t); -}; - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,9 @@ +2005-03-01 John W. Eaton <jwe@octave.org> + + * DLD-FUNCTIONS/odessa.cc: Delete. + * Makefile.in (DLD_XSRC): Remove it from the list. + (OPT_HANDLERS): Remove ODESSA-opts.cc from the list. + 2005-02-28 John W. Eaton <jwe@octave.org> * toplev.cc (octave_config_info): Remove LIBGLOB and GLOB_INCFLAGS
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -40,14 +40,14 @@ endif OPT_HANDLERS := DASPK-opts.cc DASRT-opts.cc DASSL-opts.cc \ - LSODE-opts.cc NLEqn-opts.cc ODESSA-opts.cc Quad-opts.cc + LSODE-opts.cc NLEqn-opts.cc Quad-opts.cc DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colamd.cc \ colloc.cc daspk.cc dasrt.cc dassl.cc det.cc dispatch.cc \ eig.cc expm.cc fft.cc fft2.cc fftn.cc fftw_wisdom.cc \ filter.cc find.cc fsolve.cc gammainc.cc gcd.cc getgrent.cc \ getpwent.cc getrusage.cc givens.cc hess.cc inv.cc kron.cc \ - lpsolve.cc lsode.cc lu.cc minmax.cc odessa.cc pinv.cc qr.cc \ + lpsolve.cc lsode.cc lu.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc spdet.cc \ splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l