Mercurial > octave-nkf
annotate liboctave/cruft/dassl/ddassl.f @ 19627:446c46af4b42 stable
strip trailing whitespace from most source files
* Makefile.am, NEWS, build-aux/common.mk, configure.ac,
doc/Makefile.am, doc/doxyhtml/Makefile.am,
doc/interpreter/Makefile.am, doc/interpreter/arith.txi,
doc/interpreter/audio.txi, doc/interpreter/basics.txi,
doc/interpreter/bugs.txi, doc/interpreter/container.txi,
doc/interpreter/cp-idx.txi, doc/interpreter/data.txi,
doc/interpreter/debug.txi, doc/interpreter/diagperm.txi,
doc/interpreter/diffeq.txi, doc/interpreter/doccheck/README,
doc/interpreter/doccheck/spellcheck, doc/interpreter/emacs.txi,
doc/interpreter/errors.txi, doc/interpreter/eval.txi,
doc/interpreter/expr.txi, doc/interpreter/external.txi,
doc/interpreter/fn-idx.txi, doc/interpreter/func.txi,
doc/interpreter/geometry.txi, doc/interpreter/geometryimages.m,
doc/interpreter/gpl.txi, doc/interpreter/grammar.txi,
doc/interpreter/gui.txi, doc/interpreter/image.txi,
doc/interpreter/install.txi, doc/interpreter/interp.txi,
doc/interpreter/interpimages.m, doc/interpreter/intro.txi,
doc/interpreter/io.txi, doc/interpreter/java.txi,
doc/interpreter/linalg.txi, doc/interpreter/macros.texi,
doc/interpreter/matrix.txi, doc/interpreter/munge-texi.pl,
doc/interpreter/nonlin.txi, doc/interpreter/numbers.txi,
doc/interpreter/obsolete.txi, doc/interpreter/octave-config.1,
doc/interpreter/octave.texi, doc/interpreter/oop.txi,
doc/interpreter/op-idx.txi, doc/interpreter/optim.txi,
doc/interpreter/package.txi, doc/interpreter/plot.txi,
doc/interpreter/poly.txi, doc/interpreter/preface.txi,
doc/interpreter/quad.txi, doc/interpreter/set.txi,
doc/interpreter/signal.txi, doc/interpreter/sparse.txi,
doc/interpreter/sparseimages.m, doc/interpreter/splineimages.m,
doc/interpreter/stats.txi, doc/interpreter/stmt.txi,
doc/interpreter/strings.txi, doc/interpreter/system.txi,
doc/interpreter/testfun.txi, doc/interpreter/tips.txi,
doc/interpreter/var.txi, doc/interpreter/vectorize.txi,
doc/liboctave/Makefile.am, doc/liboctave/array.texi,
doc/liboctave/bugs.texi, doc/liboctave/cp-idx.texi,
doc/liboctave/dae.texi, doc/liboctave/diffeq.texi,
doc/liboctave/error.texi, doc/liboctave/factor.texi,
doc/liboctave/fn-idx.texi, doc/liboctave/gpl.texi,
doc/liboctave/install.texi, doc/liboctave/intro.texi,
doc/liboctave/liboctave.texi, doc/liboctave/matvec.texi,
doc/liboctave/nleqn.texi, doc/liboctave/nlfunc.texi,
doc/liboctave/ode.texi, doc/liboctave/optim.texi,
doc/liboctave/preface.texi, doc/liboctave/quad.texi,
doc/liboctave/range.texi, doc/refcard/Makefile.am,
doc/refcard/refcard.tex, etc/HACKING, etc/NEWS.1, etc/NEWS.2,
etc/NEWS.3, etc/OLD-ChangeLogs/ChangeLog,
etc/OLD-ChangeLogs/doc-ChangeLog,
etc/OLD-ChangeLogs/scripts-ChangeLog,
etc/OLD-ChangeLogs/src-ChangeLog, etc/OLD-ChangeLogs/test-ChangeLog,
etc/PROJECTS, etc/README.Cygwin, etc/README.MacOS, etc/README.MinGW,
etc/README.gnuplot, etc/gdbinit, etc/icons/Makefile.am,
examples/@polynomial/end.m, examples/@polynomial/subsasgn.m,
examples/Makefile.am, examples/standalonebuiltin.cc,
libgui/Makefile.am, libgui/qterminal/libqterminal/README,
libgui/qterminal/libqterminal/unix/BlockArray.cpp,
libgui/qterminal/libqterminal/unix/BlockArray.h,
libgui/qterminal/libqterminal/unix/Character.h,
libgui/qterminal/libqterminal/unix/CharacterColor.h,
libgui/qterminal/libqterminal/unix/Emulation.cpp,
libgui/qterminal/libqterminal/unix/Emulation.h,
libgui/qterminal/libqterminal/unix/Filter.cpp,
libgui/qterminal/libqterminal/unix/Filter.h,
libgui/qterminal/libqterminal/unix/History.cpp,
libgui/qterminal/libqterminal/unix/History.h,
libgui/qterminal/libqterminal/unix/KeyboardTranslator.cpp,
libgui/qterminal/libqterminal/unix/KeyboardTranslator.h,
libgui/qterminal/libqterminal/unix/LineFont.h,
libgui/qterminal/libqterminal/unix/QUnixTerminalImpl.cpp,
libgui/qterminal/libqterminal/unix/QUnixTerminalImpl.h,
libgui/qterminal/libqterminal/unix/Screen.cpp,
libgui/qterminal/libqterminal/unix/Screen.h,
libgui/qterminal/libqterminal/unix/ScreenWindow.cpp,
libgui/qterminal/libqterminal/unix/ScreenWindow.h,
libgui/qterminal/libqterminal/unix/TerminalCharacterDecoder.cpp,
libgui/qterminal/libqterminal/unix/TerminalCharacterDecoder.h,
libgui/qterminal/libqterminal/unix/Vt102Emulation.h,
libgui/qterminal/libqterminal/win32/QWinTerminalImpl.cpp,
libgui/qterminal/qterminal/main.cpp,
libgui/src/m-editor/file-editor-tab.cc, libgui/src/octave-gui.cc,
libgui/src/octave-qt-link.cc, libinterp/corefcn/data.cc,
libinterp/corefcn/defun-int.h, libinterp/corefcn/det.cc,
libinterp/corefcn/gl2ps-renderer.cc, libinterp/corefcn/graphics.cc,
libinterp/corefcn/graphics.in.h, libinterp/corefcn/ls-mat5.cc,
libinterp/corefcn/lu.cc, libinterp/corefcn/oct-tex-parser.yy,
libinterp/corefcn/oct-tex-symbols.in, libinterp/corefcn/quadcc.cc,
libinterp/corefcn/zfstream.cc, libinterp/dldfcn/__eigs__.cc,
libinterp/dldfcn/__voronoi__.cc, libinterp/gendoc.pl,
libinterp/genprops.awk, libinterp/mk-errno-list, libinterp/mk-pkg-add,
libinterp/mkbuiltins, libinterp/mkdefs, libinterp/mkdocs,
libinterp/mkops, libinterp/octave-value/ov-java.cc,
libinterp/parse-tree/lex.ll, libinterp/parse-tree/oct-parse.in.yy,
libinterp/parse-tree/octave.gperf, liboctave/Makefile.am,
liboctave/array/Array.cc, liboctave/array/module.mk,
liboctave/cruft/daspk/datv.f, liboctave/cruft/daspk/dcnst0.f,
liboctave/cruft/daspk/dcnstr.f, liboctave/cruft/daspk/ddasic.f,
liboctave/cruft/daspk/ddasid.f, liboctave/cruft/daspk/ddasik.f,
liboctave/cruft/daspk/ddaspk.f, liboctave/cruft/daspk/ddstp.f,
liboctave/cruft/daspk/ddwnrm.f, liboctave/cruft/daspk/dfnrmd.f,
liboctave/cruft/daspk/dfnrmk.f, liboctave/cruft/daspk/dhels.f,
liboctave/cruft/daspk/dheqr.f, liboctave/cruft/daspk/dinvwt.f,
liboctave/cruft/daspk/dlinsd.f, liboctave/cruft/daspk/dlinsk.f,
liboctave/cruft/daspk/dmatd.f, liboctave/cruft/daspk/dnedd.f,
liboctave/cruft/daspk/dnedk.f, liboctave/cruft/daspk/dnsd.f,
liboctave/cruft/daspk/dnsid.f, liboctave/cruft/daspk/dnsik.f,
liboctave/cruft/daspk/dnsk.f, liboctave/cruft/daspk/dorth.f,
liboctave/cruft/daspk/dslvd.f, liboctave/cruft/daspk/dslvk.f,
liboctave/cruft/daspk/dspigm.f, liboctave/cruft/daspk/dyypnw.f,
liboctave/cruft/dasrt/ddasrt.f, liboctave/cruft/dasrt/drchek.f,
liboctave/cruft/dassl/ddaslv.f, liboctave/cruft/dassl/ddassl.f,
liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/module.mk,
liboctave/cruft/odepack/cfode.f, liboctave/cruft/odepack/dlsode.f,
liboctave/cruft/odepack/ewset.f, liboctave/cruft/odepack/intdy.f,
liboctave/cruft/odepack/prepj.f, liboctave/cruft/odepack/sintdy.f,
liboctave/cruft/odepack/slsode.f, liboctave/cruft/odepack/solsy.f,
liboctave/cruft/odepack/ssolsy.f, liboctave/cruft/odepack/stode.f,
liboctave/cruft/odepack/vnorm.f, liboctave/cruft/ranlib/Basegen.doc,
liboctave/cruft/ranlib/README, liboctave/cruft/ranlib/genbet.f,
liboctave/cruft/ranlib/genexp.f, liboctave/cruft/ranlib/gennch.f,
liboctave/cruft/ranlib/gennf.f, liboctave/cruft/ranlib/gennor.f,
liboctave/cruft/ranlib/getsd.f, liboctave/cruft/ranlib/initgn.f,
liboctave/cruft/ranlib/phrtsd.f, liboctave/cruft/ranlib/randlib.fdoc,
liboctave/cruft/ranlib/setsd.f, liboctave/cruft/ranlib/tstgmn.for,
liboctave/cruft/ranlib/tstmid.for, liboctave/cruft/slatec-fn/atanh.f,
liboctave/cruft/slatec-fn/datanh.f,
liboctave/cruft/slatec-fn/xgmainc.f,
liboctave/cruft/slatec-fn/xsgmainc.f, liboctave/numeric/module.mk,
liboctave/operators/mk-ops.awk, liboctave/operators/mx-ops,
liboctave/operators/sparse-mk-ops.awk,
liboctave/operators/sparse-mx-ops, liboctave/operators/vx-ops,
liboctave/util/module.mk, run-octave.in, scripts/@ftp/ftp.m,
scripts/audio/wavread.m, scripts/deprecated/java_convert_matrix.m,
scripts/deprecated/java_debug.m, scripts/deprecated/java_invoke.m,
scripts/deprecated/java_new.m,
scripts/deprecated/java_unsigned_conversion.m,
scripts/deprecated/javafields.m, scripts/deprecated/javamethods.m,
scripts/deprecated/shell_cmd.m, scripts/general/accumarray.m,
scripts/general/display.m, scripts/general/fieldnames.m,
scripts/general/interp1.m, scripts/general/interp2.m,
scripts/general/interp3.m, scripts/general/isa.m,
scripts/general/methods.m, scripts/general/sortrows.m,
scripts/geometry/convhull.m, scripts/geometry/delaunay.m,
scripts/geometry/delaunay3.m, scripts/geometry/delaunayn.m,
scripts/geometry/griddata.m, scripts/geometry/griddatan.m,
scripts/geometry/voronoi.m, scripts/geometry/voronoin.m,
scripts/gui/guihandles.m, scripts/gui/inputdlg.m,
scripts/gui/listdlg.m, scripts/gui/msgbox.m, scripts/gui/questdlg.m,
scripts/gui/uigetfile.m, scripts/gui/waitbar.m, scripts/gui/warndlg.m,
scripts/help/doc.m, scripts/help/help.m, scripts/help/type.m,
scripts/image/bone.m, scripts/image/cmpermute.m,
scripts/image/cmunique.m, scripts/image/colorcube.m,
scripts/image/colormap.m, scripts/image/contrast.m,
scripts/image/gray2ind.m, scripts/image/image.m,
scripts/image/imshow.m, scripts/image/ind2gray.m, scripts/image/jet.m,
scripts/image/rgb2ntsc.m, scripts/image/spinmap.m,
scripts/io/importdata.m, scripts/io/strread.m, scripts/io/textread.m,
scripts/io/textscan.m, scripts/java/java_get.m,
scripts/java/java_set.m, scripts/java/javaaddpath.m,
scripts/java/javaclasspath.m, scripts/java/javamem.m,
scripts/linear-algebra/linsolve.m, scripts/linear-algebra/qzhess.m,
scripts/miscellaneous/debug.m, scripts/miscellaneous/desktop.m,
scripts/miscellaneous/dir.m, scripts/miscellaneous/dos.m,
scripts/miscellaneous/edit.m, scripts/miscellaneous/fact.m,
scripts/miscellaneous/getappdata.m, scripts/miscellaneous/inputname.m,
scripts/miscellaneous/license.m, scripts/miscellaneous/ls_command.m,
scripts/miscellaneous/run.m, scripts/miscellaneous/setfield.m,
scripts/miscellaneous/unix.m, scripts/miscellaneous/ver.m,
scripts/mk-pkg-add, scripts/mkdoc.pl,
scripts/optimization/fminsearch.m, scripts/optimization/optimset.m,
scripts/optimization/sqp.m, scripts/pkg/pkg.m,
scripts/pkg/private/create_pkgadddel.m,
scripts/pkg/private/fix_depends.m, scripts/pkg/private/install.m,
scripts/plot/appearance/axis.m, scripts/plot/appearance/box.m,
scripts/plot/appearance/clabel.m, scripts/plot/appearance/daspect.m,
scripts/plot/appearance/datetick.m, scripts/plot/appearance/grid.m,
scripts/plot/appearance/legend.m, scripts/plot/appearance/orient.m,
scripts/plot/appearance/shading.m, scripts/plot/appearance/text.m,
scripts/plot/appearance/title.m, scripts/plot/appearance/xlabel.m,
scripts/plot/appearance/ylabel.m, scripts/plot/appearance/zlabel.m,
scripts/plot/draw/area.m, scripts/plot/draw/bar.m,
scripts/plot/draw/barh.m, scripts/plot/draw/colorbar.m,
scripts/plot/draw/contour.m, scripts/plot/draw/contour3.m,
scripts/plot/draw/contourf.m, scripts/plot/draw/ellipsoid.m,
scripts/plot/draw/errorbar.m, scripts/plot/draw/ezcontour.m,
scripts/plot/draw/ezcontourf.m, scripts/plot/draw/ezmesh.m,
scripts/plot/draw/ezpolar.m, scripts/plot/draw/fill.m,
scripts/plot/draw/fplot.m, scripts/plot/draw/hist.m,
scripts/plot/draw/meshc.m, scripts/plot/draw/meshz.m,
scripts/plot/draw/pareto.m, scripts/plot/draw/patch.m,
scripts/plot/draw/peaks.m, scripts/plot/draw/pie.m,
scripts/plot/draw/pie3.m, scripts/plot/draw/plot.m,
scripts/plot/draw/plotyy.m, scripts/plot/draw/private/__bar__.m,
scripts/plot/draw/private/__contour__.m,
scripts/plot/draw/private/__errplot__.m,
scripts/plot/draw/private/__ezplot__.m,
scripts/plot/draw/private/__patch__.m,
scripts/plot/draw/private/__stem__.m, scripts/plot/draw/rectangle.m,
scripts/plot/draw/ribbon.m, scripts/plot/draw/rose.m,
scripts/plot/draw/scatter.m, scripts/plot/draw/scatter3.m,
scripts/plot/draw/semilogx.m, scripts/plot/draw/shrinkfaces.m,
scripts/plot/draw/sombrero.m, scripts/plot/draw/sphere.m,
scripts/plot/draw/stairs.m, scripts/plot/draw/stem.m,
scripts/plot/draw/stemleaf.m, scripts/plot/draw/surf.m,
scripts/plot/draw/surface.m, scripts/plot/draw/surfc.m,
scripts/plot/draw/surfl.m, scripts/plot/draw/surfnorm.m,
scripts/plot/draw/tetramesh.m, scripts/plot/draw/trimesh.m,
scripts/plot/draw/triplot.m, scripts/plot/draw/trisurf.m,
scripts/plot/util/__gnuplot_drawnow__.m,
scripts/plot/util/__plt_get_axis_arg__.m, scripts/plot/util/axes.m,
scripts/plot/util/clf.m, scripts/plot/util/copyobj.m,
scripts/plot/util/figure.m, scripts/plot/util/gcbo.m,
scripts/plot/util/graphics_toolkit.m, scripts/plot/util/hggroup.m,
scripts/plot/util/meshgrid.m, scripts/plot/util/newplot.m,
scripts/plot/util/print.m,
scripts/plot/util/private/__add_default_menu__.m,
scripts/plot/util/private/__fltk_print__.m,
scripts/plot/util/private/__gnuplot_print__.m,
scripts/plot/util/private/__print_parse_opts__.m,
scripts/plot/util/refreshdata.m, scripts/plot/util/subplot.m,
scripts/polynomial/conv.m, scripts/polynomial/poly.m,
scripts/polynomial/polyeig.m, scripts/polynomial/polyfit.m,
scripts/polynomial/polyval.m,
scripts/polynomial/private/__splinefit__.m,
scripts/polynomial/spline.m, scripts/prefs/prefdir.m,
scripts/prefs/preferences.m, scripts/prefs/private/prefsfile.m,
scripts/prefs/rmpref.m, scripts/signal/freqz.m,
scripts/signal/module.mk, scripts/sparse/eigs.m, scripts/sparse/pcg.m,
scripts/sparse/private/__sprand_impl__.m, scripts/sparse/sprand.m,
scripts/sparse/sprandn.m, scripts/sparse/spy.m, scripts/sparse/svds.m,
scripts/specfun/expint.m, scripts/specfun/factor.m,
scripts/special-matrix/gallery.m, scripts/special-matrix/hankel.m,
scripts/special-matrix/toeplitz.m, scripts/startup/inputrc,
scripts/statistics/base/kurtosis.m, scripts/statistics/base/moment.m,
scripts/statistics/base/qqplot.m, scripts/statistics/base/var.m,
scripts/statistics/distributions/betarnd.m,
scripts/statistics/distributions/binoinv.m,
scripts/statistics/distributions/binopdf.m,
scripts/statistics/distributions/binornd.m,
scripts/statistics/distributions/cauchy_rnd.m,
scripts/statistics/distributions/chi2rnd.m,
scripts/statistics/distributions/discrete_pdf.m,
scripts/statistics/distributions/discrete_rnd.m,
scripts/statistics/distributions/empirical_rnd.m,
scripts/statistics/distributions/exprnd.m,
scripts/statistics/distributions/frnd.m,
scripts/statistics/distributions/gamrnd.m,
scripts/statistics/distributions/geornd.m,
scripts/statistics/distributions/hygernd.m,
scripts/statistics/distributions/kolmogorov_smirnov_cdf.m,
scripts/statistics/distributions/laplace_cdf.m,
scripts/statistics/distributions/laplace_pdf.m,
scripts/statistics/distributions/logistic_cdf.m,
scripts/statistics/distributions/logistic_pdf.m,
scripts/statistics/distributions/lognrnd.m,
scripts/statistics/distributions/nbincdf.m,
scripts/statistics/distributions/nbininv.m,
scripts/statistics/distributions/nbinpdf.m,
scripts/statistics/distributions/nbinrnd.m,
scripts/statistics/distributions/normrnd.m,
scripts/statistics/distributions/poissinv.m,
scripts/statistics/distributions/poissrnd.m,
scripts/statistics/distributions/tinv.m,
scripts/statistics/distributions/trnd.m,
scripts/statistics/distributions/unidcdf.m,
scripts/statistics/distributions/unidpdf.m,
scripts/statistics/distributions/unidrnd.m,
scripts/statistics/distributions/unifrnd.m,
scripts/statistics/distributions/wblrnd.m,
scripts/statistics/models/module.mk,
scripts/statistics/tests/kruskal_wallis_test.m,
scripts/strings/base2dec.m, scripts/strings/deblank.m,
scripts/strings/dec2base.m, scripts/strings/dec2bin.m,
scripts/strings/dec2hex.m, scripts/strings/mat2str.m,
scripts/strings/ostrsplit.m, scripts/strings/regexptranslate.m,
scripts/strings/str2num.m, scripts/strings/strcat.m,
scripts/strings/strjoin.m, scripts/strings/strsplit.m,
scripts/strings/strtok.m, scripts/strings/strtrim.m,
scripts/strings/strtrunc.m, scripts/strings/substr.m,
scripts/testfun/__run_test_suite__.m, scripts/testfun/speed.m,
scripts/testfun/test.m, scripts/time/asctime.m,
scripts/time/datenum.m, scripts/time/datevec.m,
scripts/time/weekday.m, src/Makefile.am, test/Makefile.am,
test/build-bc-overload-tests.sh, test/build-sparse-tests.sh,
test/jit.tst, test/line-continue.tst: Strip trailing whitespace.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 20 Jan 2015 08:26:57 -0500 |
| parents | 648dabbb4c6b |
| children |
| rev | line source |
|---|---|
| 2329 | 1 SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, |
| 2 + IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) | |
| 3 C***BEGIN PROLOGUE DDASSL | |
| 4 C***PURPOSE This code solves a system of differential/algebraic | |
| 5 C equations of the form G(T,Y,YPRIME) = 0. | |
| 6 C***LIBRARY SLATEC (DASSL) | |
| 7 C***CATEGORY I1A2 | |
| 8 C***TYPE DOUBLE PRECISION (SDASSL-S, DDASSL-D) | |
| 9 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS, | |
| 10 C IMPLICIT DIFFERENTIAL SYSTEMS | |
| 11 C***AUTHOR PETZOLD, LINDA R., (LLNL) | |
| 12 C COMPUTING AND MATHEMATICS RESEARCH DIVISION | |
| 13 C LAWRENCE LIVERMORE NATIONAL LABORATORY | |
| 14 C L - 316, P.O. BOX 808, | |
| 15 C LIVERMORE, CA. 94550 | |
| 16 C***DESCRIPTION | |
| 17 C | |
| 18 C *Usage: | |
| 19 C | |
| 20 C EXTERNAL RES, JAC | |
| 21 C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR | |
| 22 C DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL, | |
| 23 C * RWORK(LRW), RPAR | |
| 24 C | |
| 25 C CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, | |
| 26 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) | |
| 27 C | |
| 28 C | |
| 29 C *Arguments: | |
| 30 C (In the following, all real arrays should be type DOUBLE PRECISION.) | |
| 31 C | |
| 32 C RES:EXT This is a subroutine which you provide to define the | |
| 33 C differential/algebraic system. | |
| 34 C | |
| 35 C NEQ:IN This is the number of equations to be solved. | |
| 36 C | |
| 37 C T:INOUT This is the current value of the independent variable. | |
| 38 C | |
| 39 C Y(*):INOUT This array contains the solution components at T. | |
| 40 C | |
| 41 C YPRIME(*):INOUT This array contains the derivatives of the solution | |
| 42 C components at T. | |
| 43 C | |
| 44 C TOUT:IN This is a point at which a solution is desired. | |
| 45 C | |
| 46 C INFO(N):IN The basic task of the code is to solve the system from T | |
| 47 C to TOUT and return an answer at TOUT. INFO is an integer | |
| 48 C array which is used to communicate exactly how you want | |
| 49 C this task to be carried out. (See below for details.) | |
| 50 C N must be greater than or equal to 15. | |
| 51 C | |
| 52 C RTOL,ATOL:INOUT These quantities represent relative and absolute | |
| 53 C error tolerances which you provide to indicate how | |
| 54 C accurately you wish the solution to be computed. You | |
| 55 C may choose them to be both scalars or else both vectors. | |
| 56 C Caution: In Fortran 77, a scalar is not the same as an | |
| 57 C array of length 1. Some compilers may object | |
| 58 C to using scalars for RTOL,ATOL. | |
| 59 C | |
| 60 C IDID:OUT This scalar quantity is an indicator reporting what the | |
| 61 C code did. You must monitor this integer variable to | |
| 62 C decide what action to take next. | |
| 63 C | |
| 64 C RWORK:WORK A real work array of length LRW which provides the | |
| 65 C code with needed storage space. | |
| 66 C | |
| 67 C LRW:IN The length of RWORK. (See below for required length.) | |
| 68 C | |
| 69 C IWORK:WORK An integer work array of length LIW which probides the | |
| 70 C code with needed storage space. | |
| 71 C | |
| 72 C LIW:IN The length of IWORK. (See below for required length.) | |
| 73 C | |
| 74 C RPAR,IPAR:IN These are real and integer parameter arrays which | |
| 75 C you can use for communication between your calling | |
| 76 C program and the RES subroutine (and the JAC subroutine) | |
| 77 C | |
| 78 C JAC:EXT This is the name of a subroutine which you may choose | |
| 79 C to provide for defining a matrix of partial derivatives | |
| 80 C described below. | |
| 81 C | |
| 82 C Quantities which may be altered by DDASSL are: | |
| 83 C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL, | |
| 84 C IDID, RWORK(*) AND IWORK(*) | |
| 85 C | |
| 86 C *Description | |
| 87 C | |
| 88 C Subroutine DDASSL uses the backward differentiation formulas of | |
| 89 C orders one through five to solve a system of the above form for Y and | |
| 90 C YPRIME. Values for Y and YPRIME at the initial time must be given as | |
| 91 C input. These values must be consistent, (that is, if T,Y,YPRIME are | |
| 92 C the given initial values, they must satisfy G(T,Y,YPRIME) = 0.). The | |
| 93 C subroutine solves the system from T to TOUT. It is easy to continue | |
| 94 C the solution to get results at additional TOUT. This is the interval | |
| 95 C mode of operation. Intermediate results can also be obtained easily | |
| 96 C by using the intermediate-output capability. | |
| 97 C | |
| 98 C The following detailed description is divided into subsections: | |
| 99 C 1. Input required for the first call to DDASSL. | |
| 100 C 2. Output after any return from DDASSL. | |
| 101 C 3. What to do to continue the integration. | |
| 102 C 4. Error messages. | |
| 103 C | |
| 104 C | |
| 105 C -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------ | |
| 106 C | |
| 107 C The first call of the code is defined to be the start of each new | |
| 108 C problem. Read through the descriptions of all the following items, | |
| 109 C provide sufficient storage space for designated arrays, set | |
| 110 C appropriate variables for the initialization of the problem, and | |
| 111 C give information about how you want the problem to be solved. | |
| 112 C | |
| 113 C | |
| 114 C RES -- Provide a subroutine of the form | |
| 115 C SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) | |
| 116 C to define the system of differential/algebraic | |
| 117 C equations which is to be solved. For the given values | |
| 118 C of T,Y and YPRIME, the subroutine should | |
| 119 C return the residual of the defferential/algebraic | |
| 120 C system | |
| 121 C DELTA = G(T,Y,YPRIME) | |
| 122 C (DELTA(*) is a vector of length NEQ which is | |
| 123 C output for RES.) | |
| 124 C | |
| 125 C Subroutine RES must not alter T,Y or YPRIME. | |
| 126 C You must declare the name RES in an external | |
| 127 C statement in your program that calls DDASSL. | |
| 128 C You must dimension Y,YPRIME and DELTA in RES. | |
| 129 C | |
| 130 C IRES is an integer flag which is always equal to | |
| 131 C zero on input. Subroutine RES should alter IRES | |
| 132 C only if it encounters an illegal value of Y or | |
| 133 C a stop condition. Set IRES = -1 if an input value | |
| 134 C is illegal, and DDASSL will try to solve the problem | |
| 135 C without getting IRES = -1. If IRES = -2, DDASSL | |
| 136 C will return control to the calling program | |
| 137 C with IDID = -11. | |
| 138 C | |
| 139 C RPAR and IPAR are real and integer parameter arrays which | |
| 140 C you can use for communication between your calling program | |
| 141 C and subroutine RES. They are not altered by DDASSL. If you | |
| 142 C do not need RPAR or IPAR, ignore these parameters by treat- | |
| 143 C ing them as dummy arguments. If you do choose to use them, | |
| 144 C dimension them in your calling program and in RES as arrays | |
| 145 C of appropriate length. | |
| 146 C | |
| 147 C NEQ -- Set it to the number of differential equations. | |
| 148 C (NEQ .GE. 1) | |
| 149 C | |
| 150 C T -- Set it to the initial point of the integration. | |
| 151 C T must be defined as a variable. | |
| 152 C | |
| 153 C Y(*) -- Set this vector to the initial values of the NEQ solution | |
| 154 C components at the initial point. You must dimension Y of | |
| 155 C length at least NEQ in your calling program. | |
| 156 C | |
| 157 C YPRIME(*) -- Set this vector to the initial values of the NEQ | |
| 158 C first derivatives of the solution components at the initial | |
| 159 C point. You must dimension YPRIME at least NEQ in your | |
| 160 C calling program. If you do not know initial values of some | |
| 161 C of the solution components, see the explanation of INFO(11). | |
| 162 C | |
| 163 C TOUT -- Set it to the first point at which a solution | |
| 164 C is desired. You can not take TOUT = T. | |
| 165 C integration either forward in T (TOUT .GT. T) or | |
| 166 C backward in T (TOUT .LT. T) is permitted. | |
| 167 C | |
| 168 C The code advances the solution from T to TOUT using | |
| 169 C step sizes which are automatically selected so as to | |
| 170 C achieve the desired accuracy. If you wish, the code will | |
| 171 C return with the solution and its derivative at | |
| 172 C intermediate steps (intermediate-output mode) so that | |
| 173 C you can monitor them, but you still must provide TOUT in | |
| 174 C accord with the basic aim of the code. | |
| 175 C | |
| 176 C The first step taken by the code is a critical one | |
| 177 C because it must reflect how fast the solution changes near | |
| 178 C the initial point. The code automatically selects an | |
| 179 C initial step size which is practically always suitable for | |
| 180 C the problem. By using the fact that the code will not step | |
| 181 C past TOUT in the first step, you could, if necessary, | |
| 182 C restrict the length of the initial step size. | |
| 183 C | |
| 184 C For some problems it may not be permissible to integrate | |
| 185 C past a point TSTOP because a discontinuity occurs there | |
| 186 C or the solution or its derivative is not defined beyond | |
| 187 C TSTOP. When you have declared a TSTOP point (SEE INFO(4) | |
| 188 C and RWORK(1)), you have told the code not to integrate | |
| 189 C past TSTOP. In this case any TOUT beyond TSTOP is invalid | |
| 190 C input. | |
| 191 C | |
| 192 C INFO(*) -- Use the INFO array to give the code more details about | |
| 193 C how you want your problem solved. This array should be | |
| 194 C dimensioned of length 15, though DDASSL uses only the first | |
| 195 C eleven entries. You must respond to all of the following | |
| 196 C items, which are arranged as questions. The simplest use | |
| 197 C of the code corresponds to answering all questions as yes, | |
| 198 C i.e. setting all entries of INFO to 0. | |
| 199 C | |
| 200 C INFO(1) - This parameter enables the code to initialize | |
| 201 C itself. You must set it to indicate the start of every | |
| 202 C new problem. | |
| 203 C | |
| 204 C **** Is this the first call for this problem ... | |
| 205 C Yes - Set INFO(1) = 0 | |
| 206 C No - Not applicable here. | |
| 207 C See below for continuation calls. **** | |
| 208 C | |
| 209 C INFO(2) - How much accuracy you want of your solution | |
| 210 C is specified by the error tolerances RTOL and ATOL. | |
| 211 C The simplest use is to take them both to be scalars. | |
| 212 C To obtain more flexibility, they can both be vectors. | |
| 213 C The code must be told your choice. | |
| 214 C | |
| 215 C **** Are both error tolerances RTOL, ATOL scalars ... | |
| 216 C Yes - Set INFO(2) = 0 | |
| 217 C and input scalars for both RTOL and ATOL | |
| 218 C No - Set INFO(2) = 1 | |
| 219 C and input arrays for both RTOL and ATOL **** | |
| 220 C | |
| 221 C INFO(3) - The code integrates from T in the direction | |
| 222 C of TOUT by steps. If you wish, it will return the | |
| 223 C computed solution and derivative at the next | |
| 224 C intermediate step (the intermediate-output mode) or | |
| 225 C TOUT, whichever comes first. This is a good way to | |
| 226 C proceed if you want to see the behavior of the solution. | |
| 227 C If you must have solutions at a great many specific | |
| 228 C TOUT points, this code will compute them efficiently. | |
| 229 C | |
| 230 C **** Do you want the solution only at | |
| 231 C TOUT (and not at the next intermediate step) ... | |
| 232 C Yes - Set INFO(3) = 0 | |
| 233 C No - Set INFO(3) = 1 **** | |
| 234 C | |
| 235 C INFO(4) - To handle solutions at a great many specific | |
| 236 C values TOUT efficiently, this code may integrate past | |
| 237 C TOUT and interpolate to obtain the result at TOUT. | |
| 238 C Sometimes it is not possible to integrate beyond some | |
| 239 C point TSTOP because the equation changes there or it is | |
| 240 C not defined past TSTOP. Then you must tell the code | |
| 241 C not to go past. | |
| 242 C | |
| 243 C **** Can the integration be carried out without any | |
| 244 C restrictions on the independent variable T ... | |
| 245 C Yes - Set INFO(4)=0 | |
| 246 C No - Set INFO(4)=1 | |
| 247 C and define the stopping point TSTOP by | |
| 248 C setting RWORK(1)=TSTOP **** | |
| 249 C | |
| 250 C INFO(5) - To solve differential/algebraic problems it is | |
| 251 C necessary to use a matrix of partial derivatives of the | |
| 252 C system of differential equations. If you do not | |
| 253 C provide a subroutine to evaluate it analytically (see | |
| 254 C description of the item JAC in the call list), it will | |
| 255 C be approximated by numerical differencing in this code. | |
| 256 C although it is less trouble for you to have the code | |
| 257 C compute partial derivatives by numerical differencing, | |
| 258 C the solution will be more reliable if you provide the | |
| 259 C derivatives via JAC. Sometimes numerical differencing | |
| 260 C is cheaper than evaluating derivatives in JAC and | |
| 261 C sometimes it is not - this depends on your problem. | |
| 262 C | |
| 263 C **** Do you want the code to evaluate the partial | |
| 264 C derivatives automatically by numerical differences ... | |
| 265 C Yes - Set INFO(5)=0 | |
| 266 C No - Set INFO(5)=1 | |
| 267 C and provide subroutine JAC for evaluating the | |
| 268 C matrix of partial derivatives **** | |
| 269 C | |
| 270 C INFO(6) - DDASSL will perform much better if the matrix of | |
| 271 C partial derivatives, DG/DY + CJ*DG/DYPRIME, | |
| 272 C (here CJ is a scalar determined by DDASSL) | |
| 273 C is banded and the code is told this. In this | |
| 274 C case, the storage needed will be greatly reduced, | |
| 275 C numerical differencing will be performed much cheaper, | |
| 276 C and a number of important algorithms will execute much | |
| 277 C faster. The differential equation is said to have | |
| 278 C half-bandwidths ML (lower) and MU (upper) if equation i | |
| 279 C involves only unknowns Y(J) with | |
| 280 C I-ML .LE. J .LE. I+MU | |
| 281 C for all I=1,2,...,NEQ. Thus, ML and MU are the widths | |
| 282 C of the lower and upper parts of the band, respectively, | |
| 283 C with the main diagonal being excluded. If you do not | |
| 284 C indicate that the equation has a banded matrix of partial | |
| 285 C derivatives, the code works with a full matrix of NEQ**2 | |
| 286 C elements (stored in the conventional way). Computations | |
| 287 C with banded matrices cost less time and storage than with | |
| 288 C full matrices if 2*ML+MU .LT. NEQ. If you tell the | |
| 289 C code that the matrix of partial derivatives has a banded | |
| 290 C structure and you want to provide subroutine JAC to | |
| 291 C compute the partial derivatives, then you must be careful | |
| 292 C to store the elements of the matrix in the special form | |
| 293 C indicated in the description of JAC. | |
| 294 C | |
| 295 C **** Do you want to solve the problem using a full | |
| 296 C (dense) matrix (and not a special banded | |
| 297 C structure) ... | |
| 298 C Yes - Set INFO(6)=0 | |
| 299 C No - Set INFO(6)=1 | |
| 300 C and provide the lower (ML) and upper (MU) | |
| 301 C bandwidths by setting | |
| 302 C IWORK(1)=ML | |
| 303 C IWORK(2)=MU **** | |
| 304 C | |
| 305 C | |
| 306 C INFO(7) -- You can specify a maximum (absolute value of) | |
| 307 C stepsize, so that the code | |
| 308 C will avoid passing over very | |
| 309 C large regions. | |
| 310 C | |
| 311 C **** Do you want the code to decide | |
| 312 C on its own maximum stepsize? | |
| 313 C Yes - Set INFO(7)=0 | |
| 314 C No - Set INFO(7)=1 | |
| 315 C and define HMAX by setting | |
| 316 C RWORK(2)=HMAX **** | |
| 317 C | |
| 318 C INFO(8) -- Differential/algebraic problems | |
| 319 C may occaisionally suffer from | |
| 320 C severe scaling difficulties on the | |
| 321 C first step. If you know a great deal | |
| 322 C about the scaling of your problem, you can | |
| 323 C help to alleviate this problem by | |
| 324 C specifying an initial stepsize HO. | |
| 325 C | |
| 326 C **** Do you want the code to define | |
| 327 C its own initial stepsize? | |
| 328 C Yes - Set INFO(8)=0 | |
| 329 C No - Set INFO(8)=1 | |
| 330 C and define HO by setting | |
| 331 C RWORK(3)=HO **** | |
| 332 C | |
| 333 C INFO(9) -- If storage is a severe problem, | |
| 334 C you can save some locations by | |
| 335 C restricting the maximum order MAXORD. | |
| 336 C the default value is 5. for each | |
| 337 C order decrease below 5, the code | |
| 338 C requires NEQ fewer locations, however | |
| 339 C it is likely to be slower. In any | |
| 340 C case, you must have 1 .LE. MAXORD .LE. 5 | |
| 341 C **** Do you want the maximum order to | |
| 342 C default to 5? | |
| 343 C Yes - Set INFO(9)=0 | |
| 344 C No - Set INFO(9)=1 | |
| 345 C and define MAXORD by setting | |
| 346 C IWORK(3)=MAXORD **** | |
| 347 C | |
| 348 C INFO(10) --If you know that the solutions to your equations | |
| 349 C will always be nonnegative, it may help to set this | |
| 350 C parameter. However, it is probably best to | |
| 351 C try the code without using this option first, | |
| 352 C and only to use this option if that doesn't | |
| 353 C work very well. | |
| 354 C **** Do you want the code to solve the problem without | |
| 355 C invoking any special nonnegativity constraints? | |
| 356 C Yes - Set INFO(10)=0 | |
| 357 C No - Set INFO(10)=1 | |
| 358 C | |
| 359 C INFO(11) --DDASSL normally requires the initial T, | |
| 360 C Y, and YPRIME to be consistent. That is, | |
| 361 C you must have G(T,Y,YPRIME) = 0 at the initial | |
| 362 C time. If you do not know the initial | |
| 363 C derivative precisely, you can let DDASSL try | |
| 364 C to compute it. | |
| 365 C **** Are the initialHE INITIAL T, Y, YPRIME consistent? | |
| 366 C Yes - Set INFO(11) = 0 | |
| 367 C No - Set INFO(11) = 1, | |
| 368 C and set YPRIME to an initial approximation | |
| 369 C to YPRIME. (If you have no idea what | |
| 370 C YPRIME should be, set it to zero. Note | |
| 371 C that the initial Y should be such | |
| 372 C that there must exist a YPRIME so that | |
| 373 C G(T,Y,YPRIME) = 0.) | |
| 374 C | |
| 375 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL | |
| 376 C error tolerances to tell the code how accurately you | |
| 377 C want the solution to be computed. They must be defined | |
| 378 C as variables because the code may change them. You | |
| 379 C have two choices -- | |
| 380 C Both RTOL and ATOL are scalars. (INFO(2)=0) | |
| 381 C Both RTOL and ATOL are vectors. (INFO(2)=1) | |
| 382 C in either case all components must be non-negative. | |
| 383 C | |
| 384 C The tolerances are used by the code in a local error | |
| 385 C test at each step which requires roughly that | |
| 386 C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL | |
| 387 C for each vector component. | |
| 388 C (More specifically, a root-mean-square norm is used to | |
| 389 C measure the size of vectors, and the error test uses the | |
| 390 C magnitude of the solution at the beginning of the step.) | |
| 391 C | |
| 392 C The true (global) error is the difference between the | |
| 393 C true solution of the initial value problem and the | |
| 394 C computed approximation. Practically all present day | |
| 395 C codes, including this one, control the local error at | |
| 396 C each step and do not even attempt to control the global | |
| 397 C error directly. | |
| 398 C Usually, but not always, the true accuracy of the | |
| 399 C computed Y is comparable to the error tolerances. This | |
| 400 C code will usually, but not always, deliver a more | |
| 401 C accurate solution if you reduce the tolerances and | |
| 402 C integrate again. By comparing two such solutions you | |
| 403 C can get a fairly reliable idea of the true error in the | |
| 404 C solution at the bigger tolerances. | |
| 405 C | |
| 406 C Setting ATOL=0. results in a pure relative error test on | |
| 407 C that component. Setting RTOL=0. results in a pure | |
| 408 C absolute error test on that component. A mixed test | |
| 409 C with non-zero RTOL and ATOL corresponds roughly to a | |
| 410 C relative error test when the solution component is much | |
| 411 C bigger than ATOL and to an absolute error test when the | |
| 412 C solution component is smaller than the threshhold ATOL. | |
| 413 C | |
| 414 C The code will not attempt to compute a solution at an | |
| 415 C accuracy unreasonable for the machine being used. It will | |
| 416 C advise you if you ask for too much accuracy and inform | |
| 417 C you as to the maximum accuracy it believes possible. | |
| 418 C | |
| 419 C RWORK(*) -- Dimension this real work array of length LRW in your | |
| 420 C calling program. | |
| 421 C | |
| 422 C LRW -- Set it to the declared length of the RWORK array. | |
| 423 C You must have | |
| 424 C LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2 | |
| 425 C for the full (dense) JACOBIAN case (when INFO(6)=0), or | |
| 426 C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ | |
| 427 C for the banded user-defined JACOBIAN case | |
| 428 C (when INFO(5)=1 and INFO(6)=1), or | |
| 429 C LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ | |
| 430 C +2*(NEQ/(ML+MU+1)+1) | |
| 431 C for the banded finite-difference-generated JACOBIAN case | |
| 432 C (when INFO(5)=0 and INFO(6)=1) | |
| 433 C | |
| 434 C IWORK(*) -- Dimension this integer work array of length LIW in | |
| 435 C your calling program. | |
| 436 C | |
| 437 C LIW -- Set it to the declared length of the IWORK array. | |
| 4429 | 438 C You must have LIW .GE. 21+NEQ |
| 2329 | 439 C |
| 440 C RPAR, IPAR -- These are parameter arrays, of real and integer | |
| 441 C type, respectively. You can use them for communication | |
| 442 C between your program that calls DDASSL and the | |
| 443 C RES subroutine (and the JAC subroutine). They are not | |
| 444 C altered by DDASSL. If you do not need RPAR or IPAR, | |
| 445 C ignore these parameters by treating them as dummy | |
| 446 C arguments. If you do choose to use them, dimension | |
| 447 C them in your calling program and in RES (and in JAC) | |
| 448 C as arrays of appropriate length. | |
| 449 C | |
| 450 C JAC -- If you have set INFO(5)=0, you can ignore this parameter | |
| 451 C by treating it as a dummy argument. Otherwise, you must | |
| 452 C provide a subroutine of the form | |
| 453 C SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR) | |
| 454 C to define the matrix of partial derivatives | |
| 455 C PD=DG/DY+CJ*DG/DYPRIME | |
| 456 C CJ is a scalar which is input to JAC. | |
| 457 C For the given values of T,Y,YPRIME, the | |
| 458 C subroutine must evaluate the non-zero partial | |
| 459 C derivatives for each equation and each solution | |
| 460 C component, and store these values in the | |
| 461 C matrix PD. The elements of PD are set to zero | |
| 462 C before each call to JAC so only non-zero elements | |
| 463 C need to be defined. | |
| 464 C | |
| 465 C Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ. | |
| 466 C You must declare the name JAC in an EXTERNAL statement in | |
| 467 C your program that calls DDASSL. You must dimension Y, | |
| 468 C YPRIME and PD in JAC. | |
| 469 C | |
| 470 C The way you must store the elements into the PD matrix | |
| 471 C depends on the structure of the matrix which you | |
| 472 C indicated by INFO(6). | |
| 473 C *** INFO(6)=0 -- Full (dense) matrix *** | |
| 474 C Give PD a first dimension of NEQ. | |
| 475 C When you evaluate the (non-zero) partial derivative | |
| 476 C of equation I with respect to variable J, you must | |
| 477 C store it in PD according to | |
| 478 C PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)" | |
| 479 C *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU | |
| 480 C upper diagonal bands (refer to INFO(6) description | |
| 481 C of ML and MU) *** | |
| 482 C Give PD a first dimension of 2*ML+MU+1. | |
| 483 C when you evaluate the (non-zero) partial derivative | |
| 484 C of equation I with respect to variable J, you must | |
| 485 C store it in PD according to | |
| 486 C IROW = I - J + ML + MU + 1 | |
| 487 C PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)" | |
| 488 C | |
| 489 C RPAR and IPAR are real and integer parameter arrays | |
| 490 C which you can use for communication between your calling | |
| 491 C program and your JACOBIAN subroutine JAC. They are not | |
| 492 C altered by DDASSL. If you do not need RPAR or IPAR, | |
| 493 C ignore these parameters by treating them as dummy | |
| 494 C arguments. If you do choose to use them, dimension | |
| 495 C them in your calling program and in JAC as arrays of | |
| 496 C appropriate length. | |
| 497 C | |
| 498 C | |
| 499 C OPTIONALLY REPLACEABLE NORM ROUTINE: | |
| 500 C | |
| 501 C DDASSL uses a weighted norm DDANRM to measure the size | |
| 502 C of vectors such as the estimated error in each step. | |
| 503 C A FUNCTION subprogram | |
| 504 C DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR) | |
| 505 C DIMENSION V(NEQ),WT(NEQ) | |
| 506 C is used to define this norm. Here, V is the vector | |
| 507 C whose norm is to be computed, and WT is a vector of | |
| 508 C weights. A DDANRM routine has been included with DDASSL | |
| 509 C which computes the weighted root-mean-square norm | |
| 510 C given by | |
| 511 C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2) | |
| 512 C this norm is suitable for most problems. In some | |
| 513 C special cases, it may be more convenient and/or | |
| 514 C efficient to define your own norm by writing a function | |
| 515 C subprogram to be called instead of DDANRM. This should, | |
| 516 C however, be attempted only after careful thought and | |
| 517 C consideration. | |
| 518 C | |
| 519 C | |
| 520 C -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL --------------------- | |
| 521 C | |
| 522 C The principal aim of the code is to return a computed solution at | |
| 523 C TOUT, although it is also possible to obtain intermediate results | |
| 524 C along the way. To find out whether the code achieved its goal | |
| 525 C or if the integration process was interrupted before the task was | |
| 526 C completed, you must check the IDID parameter. | |
| 527 C | |
| 528 C | |
| 529 C T -- The solution was successfully advanced to the | |
| 530 C output value of T. | |
| 531 C | |
| 532 C Y(*) -- Contains the computed solution approximation at T. | |
| 533 C | |
| 534 C YPRIME(*) -- Contains the computed derivative | |
| 535 C approximation at T. | |
| 536 C | |
| 537 C IDID -- Reports what the code did. | |
| 538 C | |
| 539 C *** Task completed *** | |
| 540 C Reported by positive values of IDID | |
| 541 C | |
| 542 C IDID = 1 -- A step was successfully taken in the | |
| 543 C intermediate-output mode. The code has not | |
| 544 C yet reached TOUT. | |
| 545 C | |
| 546 C IDID = 2 -- The integration to TSTOP was successfully | |
| 547 C completed (T=TSTOP) by stepping exactly to TSTOP. | |
| 548 C | |
| 549 C IDID = 3 -- The integration to TOUT was successfully | |
| 550 C completed (T=TOUT) by stepping past TOUT. | |
| 551 C Y(*) is obtained by interpolation. | |
| 552 C YPRIME(*) is obtained by interpolation. | |
| 553 C | |
| 554 C *** Task interrupted *** | |
| 555 C Reported by negative values of IDID | |
| 556 C | |
| 557 C IDID = -1 -- A large amount of work has been expended. | |
| 558 C (About 500 steps) | |
| 559 C | |
| 560 C IDID = -2 -- The error tolerances are too stringent. | |
| 561 C | |
| 562 C IDID = -3 -- The local error test cannot be satisfied | |
| 563 C because you specified a zero component in ATOL | |
| 564 C and the corresponding computed solution | |
| 565 C component is zero. Thus, a pure relative error | |
| 566 C test is impossible for this component. | |
| 567 C | |
| 568 C IDID = -6 -- DDASSL had repeated error test | |
| 569 C failures on the last attempted step. | |
| 570 C | |
| 571 C IDID = -7 -- The corrector could not converge. | |
| 572 C | |
| 573 C IDID = -8 -- The matrix of partial derivatives | |
| 574 C is singular. | |
| 575 C | |
| 576 C IDID = -9 -- The corrector could not converge. | |
| 577 C there were repeated error test failures | |
| 578 C in this step. | |
| 579 C | |
| 580 C IDID =-10 -- The corrector could not converge | |
| 581 C because IRES was equal to minus one. | |
| 582 C | |
| 583 C IDID =-11 -- IRES equal to -2 was encountered | |
| 584 C and control is being returned to the | |
| 585 C calling program. | |
| 586 C | |
| 587 C IDID =-12 -- DDASSL failed to compute the initial | |
| 588 C YPRIME. | |
| 589 C | |
| 590 C | |
| 591 C | |
| 592 C IDID = -13,..,-32 -- Not applicable for this code | |
| 593 C | |
| 594 C *** Task terminated *** | |
| 595 C Reported by the value of IDID=-33 | |
| 596 C | |
| 597 C IDID = -33 -- The code has encountered trouble from which | |
| 598 C it cannot recover. A message is printed | |
| 599 C explaining the trouble and control is returned | |
| 600 C to the calling program. For example, this occurs | |
| 601 C when invalid input is detected. | |
| 602 C | |
| 603 C RTOL, ATOL -- These quantities remain unchanged except when | |
| 604 C IDID = -2. In this case, the error tolerances have been | |
| 605 C increased by the code to values which are estimated to | |
| 606 C be appropriate for continuing the integration. However, | |
| 607 C the reported solution at T was obtained using the input | |
| 608 C values of RTOL and ATOL. | |
| 609 C | |
| 610 C RWORK, IWORK -- Contain information which is usually of no | |
| 611 C interest to the user but necessary for subsequent calls. | |
| 612 C However, you may find use for | |
| 613 C | |
| 614 C RWORK(3)--Which contains the step size H to be | |
| 615 C attempted on the next step. | |
| 616 C | |
| 617 C RWORK(4)--Which contains the current value of the | |
| 618 C independent variable, i.e., the farthest point | |
| 619 C integration has reached. This will be different | |
| 620 C from T only when interpolation has been | |
| 621 C performed (IDID=3). | |
| 622 C | |
| 623 C RWORK(7)--Which contains the stepsize used | |
| 624 C on the last successful step. | |
| 625 C | |
| 626 C IWORK(7)--Which contains the order of the method to | |
| 627 C be attempted on the next step. | |
| 628 C | |
| 629 C IWORK(8)--Which contains the order of the method used | |
| 630 C on the last step. | |
| 631 C | |
| 632 C IWORK(11)--Which contains the number of steps taken so | |
| 633 C far. | |
| 634 C | |
| 635 C IWORK(12)--Which contains the number of calls to RES | |
| 636 C so far. | |
| 637 C | |
| 638 C IWORK(13)--Which contains the number of evaluations of | |
| 639 C the matrix of partial derivatives needed so | |
| 640 C far. | |
| 641 C | |
| 642 C IWORK(14)--Which contains the total number | |
| 643 C of error test failures so far. | |
| 644 C | |
| 645 C IWORK(15)--Which contains the total number | |
| 646 C of convergence test failures so far. | |
| 647 C (includes singular iteration matrix | |
| 648 C failures.) | |
| 649 C | |
| 650 C | |
| 651 C -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------ | |
| 652 C (CALLS AFTER THE FIRST) | |
| 653 C | |
| 654 C This code is organized so that subsequent calls to continue the | |
| 655 C integration involve little (if any) additional effort on your | |
| 656 C part. You must monitor the IDID parameter in order to determine | |
| 657 C what to do next. | |
| 658 C | |
| 659 C Recalling that the principal task of the code is to integrate | |
| 660 C from T to TOUT (the interval mode), usually all you will need | |
| 661 C to do is specify a new TOUT upon reaching the current TOUT. | |
| 662 C | |
| 663 C Do not alter any quantity not specifically permitted below, | |
| 664 C in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*) | |
| 665 C or the differential equation in subroutine RES. Any such | |
| 666 C alteration constitutes a new problem and must be treated as such, | |
| 667 C i.e., you must start afresh. | |
| 668 C | |
| 669 C You cannot change from vector to scalar error control or vice | |
| 670 C versa (INFO(2)), but you can change the size of the entries of | |
| 671 C RTOL, ATOL. Increasing a tolerance makes the equation easier | |
| 672 C to integrate. Decreasing a tolerance will make the equation | |
| 673 C harder to integrate and should generally be avoided. | |
| 674 C | |
| 675 C You can switch from the intermediate-output mode to the | |
| 676 C interval mode (INFO(3)) or vice versa at any time. | |
| 677 C | |
| 678 C If it has been necessary to prevent the integration from going | |
| 679 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the | |
| 680 C code will not integrate to any TOUT beyond the currently | |
| 681 C specified TSTOP. Once TSTOP has been reached you must change | |
| 682 C the value of TSTOP or set INFO(4)=0. You may change INFO(4) | |
| 683 C or TSTOP at any time but you must supply the value of TSTOP in | |
| 684 C RWORK(1) whenever you set INFO(4)=1. | |
| 685 C | |
| 686 C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2) | |
| 687 C unless you are going to restart the code. | |
| 688 C | |
| 689 C *** Following a completed task *** | |
| 690 C If | |
| 691 C IDID = 1, call the code again to continue the integration | |
| 692 C another step in the direction of TOUT. | |
| 693 C | |
| 694 C IDID = 2 or 3, define a new TOUT and call the code again. | |
| 695 C TOUT must be different from T. You cannot change | |
| 696 C the direction of integration without restarting. | |
| 697 C | |
| 698 C *** Following an interrupted task *** | |
| 699 C To show the code that you realize the task was | |
| 700 C interrupted and that you want to continue, you | |
| 701 C must take appropriate action and set INFO(1) = 1 | |
| 702 C If | |
| 703 C IDID = -1, The code has taken about 500 steps. | |
| 704 C If you want to continue, set INFO(1) = 1 and | |
| 705 C call the code again. An additional 500 steps | |
| 706 C will be allowed. | |
| 707 C | |
| 708 C IDID = -2, The error tolerances RTOL, ATOL have been | |
| 709 C increased to values the code estimates appropriate | |
| 710 C for continuing. You may want to change them | |
| 711 C yourself. If you are sure you want to continue | |
| 712 C with relaxed error tolerances, set INFO(1)=1 and | |
| 713 C call the code again. | |
| 714 C | |
| 715 C IDID = -3, A solution component is zero and you set the | |
| 716 C corresponding component of ATOL to zero. If you | |
| 717 C are sure you want to continue, you must first | |
| 718 C alter the error criterion to use positive values | |
| 719 C for those components of ATOL corresponding to zero | |
| 720 C solution components, then set INFO(1)=1 and call | |
| 721 C the code again. | |
| 722 C | |
| 723 C IDID = -4,-5 --- Cannot occur with this code. | |
| 724 C | |
| 725 C IDID = -6, Repeated error test failures occurred on the | |
| 726 C last attempted step in DDASSL. A singularity in the | |
| 727 C solution may be present. If you are absolutely | |
| 728 C certain you want to continue, you should restart | |
| 729 C the integration. (Provide initial values of Y and | |
| 730 C YPRIME which are consistent) | |
| 731 C | |
| 732 C IDID = -7, Repeated convergence test failures occurred | |
| 733 C on the last attempted step in DDASSL. An inaccurate | |
| 734 C or ill-conditioned JACOBIAN may be the problem. If | |
| 735 C you are absolutely certain you want to continue, you | |
| 736 C should restart the integration. | |
| 737 C | |
| 738 C IDID = -8, The matrix of partial derivatives is singular. | |
| 739 C Some of your equations may be redundant. | |
| 740 C DDASSL cannot solve the problem as stated. | |
| 741 C It is possible that the redundant equations | |
| 742 C could be removed, and then DDASSL could | |
| 743 C solve the problem. It is also possible | |
| 744 C that a solution to your problem either | |
| 745 C does not exist or is not unique. | |
| 746 C | |
| 747 C IDID = -9, DDASSL had multiple convergence test | |
| 748 C failures, preceeded by multiple error | |
| 749 C test failures, on the last attempted step. | |
| 750 C It is possible that your problem | |
| 751 C is ill-posed, and cannot be solved | |
| 752 C using this code. Or, there may be a | |
| 753 C discontinuity or a singularity in the | |
| 754 C solution. If you are absolutely certain | |
| 755 C you want to continue, you should restart | |
| 756 C the integration. | |
| 757 C | |
| 758 C IDID =-10, DDASSL had multiple convergence test failures | |
| 759 C because IRES was equal to minus one. | |
| 760 C If you are absolutely certain you want | |
| 761 C to continue, you should restart the | |
| 762 C integration. | |
| 763 C | |
| 764 C IDID =-11, IRES=-2 was encountered, and control is being | |
| 765 C returned to the calling program. | |
| 766 C | |
| 767 C IDID =-12, DDASSL failed to compute the initial YPRIME. | |
| 768 C This could happen because the initial | |
| 769 C approximation to YPRIME was not very good, or | |
| 770 C if a YPRIME consistent with the initial Y | |
| 771 C does not exist. The problem could also be caused | |
| 772 C by an inaccurate or singular iteration matrix. | |
| 773 C | |
| 774 C IDID = -13,..,-32 --- Cannot occur with this code. | |
| 775 C | |
| 776 C | |
| 777 C *** Following a terminated task *** | |
| 778 C | |
| 779 C If IDID= -33, you cannot continue the solution of this problem. | |
| 780 C An attempt to do so will result in your | |
| 781 C run being terminated. | |
| 782 C | |
| 783 C | |
| 784 C -------- ERROR MESSAGES --------------------------------------------- | |
| 785 C | |
| 786 C The SLATEC error print routine XERMSG is called in the event of | |
| 787 C unsuccessful completion of a task. Most of these are treated as | |
| 788 C "recoverable errors", which means that (unless the user has directed | |
| 789 C otherwise) control will be returned to the calling program for | |
| 790 C possible action after the message has been printed. | |
| 791 C | |
| 792 C In the event of a negative value of IDID other than -33, an appro- | |
| 793 C priate message is printed and the "error number" printed by XERMSG | |
| 794 C is the value of IDID. There are quite a number of illegal input | |
| 795 C errors that can lead to a returned value IDID=-33. The conditions | |
| 796 C and their printed "error numbers" are as follows: | |
| 797 C | |
| 798 C Error number Condition | |
| 799 C | |
| 800 C 1 Some element of INFO vector is not zero or one. | |
| 801 C 2 NEQ .le. 0 | |
| 802 C 3 MAXORD not in range. | |
| 803 C 4 LRW is less than the required length for RWORK. | |
| 804 C 5 LIW is less than the required length for IWORK. | |
| 805 C 6 Some element of RTOL is .lt. 0 | |
| 806 C 7 Some element of ATOL is .lt. 0 | |
| 807 C 8 All elements of RTOL and ATOL are zero. | |
| 808 C 9 INFO(4)=1 and TSTOP is behind TOUT. | |
| 809 C 10 HMAX .lt. 0.0 | |
| 810 C 11 TOUT is behind T. | |
| 811 C 12 INFO(8)=1 and H0=0.0 | |
| 812 C 13 Some element of WT is .le. 0.0 | |
| 813 C 14 TOUT is too close to T to start integration. | |
| 814 C 15 INFO(4)=1 and TSTOP is behind T. | |
| 815 C 16 --( Not used in this version )-- | |
| 816 C 17 ML illegal. Either .lt. 0 or .gt. NEQ | |
| 817 C 18 MU illegal. Either .lt. 0 or .gt. NEQ | |
| 818 C 19 TOUT = T. | |
| 819 C | |
| 820 C If DDASSL is called again without any action taken to remove the | |
| 821 C cause of an unsuccessful return, XERMSG will be called with a fatal | |
| 822 C error flag, which will cause unconditional termination of the | |
| 823 C program. There are two such fatal errors: | |
| 824 C | |
| 825 C Error number -998: The last step was terminated with a negative | |
| 826 C value of IDID other than -33, and no appropriate action was | |
| 827 C taken. | |
| 828 C | |
| 829 C Error number -999: The previous call was terminated because of | |
| 830 C illegal input (IDID=-33) and there is illegal input in the | |
| 831 C present call, as well. (Suspect infinite loop.) | |
| 832 C | |
| 833 C --------------------------------------------------------------------- | |
| 834 C | |
| 835 C***REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC | |
| 836 C SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637, | |
| 837 C SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982. | |
| 838 C***ROUTINES CALLED D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, | |
| 839 C XERMSG | |
| 840 C***REVISION HISTORY (YYMMDD) | |
| 841 C 830315 DATE WRITTEN | |
| 842 C 880387 Code changes made. All common statements have been | |
| 843 C replaced by a DATA statement, which defines pointers into | |
| 844 C RWORK, and PARAMETER statements which define pointers | |
| 845 C into IWORK. As well the documentation has gone through | |
| 846 C grammatical changes. | |
| 847 C 881005 The prologue has been changed to mixed case. | |
| 848 C The subordinate routines had revision dates changed to | |
| 849 C this date, although the documentation for these routines | |
| 850 C is all upper case. No code changes. | |
| 851 C 890511 Code changes made. The DATA statement in the declaration | |
| 852 C section of DDASSL was replaced with a PARAMETER | |
| 853 C statement. Also the statement S = 100.D0 was removed | |
| 854 C from the top of the Newton iteration in DDASTP. | |
| 855 C The subordinate routines had revision dates changed to | |
| 856 C this date. | |
| 857 C 890517 The revision date syntax was replaced with the revision | |
| 858 C history syntax. Also the "DECK" comment was added to | |
| 859 C the top of all subroutines. These changes are consistent | |
| 860 C with new SLATEC guidelines. | |
| 861 C The subordinate routines had revision dates changed to | |
| 862 C this date. No code changes. | |
| 863 C 891013 Code changes made. | |
| 864 C Removed all occurrances of FLOAT or DBLE. All operations | |
| 865 C are now performed with "mixed-mode" arithmetic. | |
| 866 C Also, specific function names were replaced with generic | |
| 867 C function names to be consistent with new SLATEC guidelines. | |
| 868 C In particular: | |
| 869 C Replaced DSQRT with SQRT everywhere. | |
| 870 C Replaced DABS with ABS everywhere. | |
| 871 C Replaced DMIN1 with MIN everywhere. | |
| 872 C Replaced MIN0 with MIN everywhere. | |
| 873 C Replaced DMAX1 with MAX everywhere. | |
| 874 C Replaced MAX0 with MAX everywhere. | |
| 875 C Replaced DSIGN with SIGN everywhere. | |
| 876 C Also replaced REVISION DATE with REVISION HISTORY in all | |
| 877 C subordinate routines. | |
| 878 C 901004 Miscellaneous changes to prologue to complete conversion | |
| 879 C to SLATEC 4.0 format. No code changes. (F.N.Fritsch) | |
| 880 C 901009 Corrected GAMS classification code and converted subsidiary | |
| 881 C routines to 4.0 format. No code changes. (F.N.Fritsch) | |
| 882 C 901010 Converted XERRWV calls to XERMSG calls. (R.Clemens,AFWL) | |
| 883 C 901019 Code changes made. | |
| 884 C Merged SLATEC 4.0 changes with previous changes made | |
| 885 C by C. Ulrich. Below is a history of the changes made by | |
| 886 C C. Ulrich. (Changes in subsidiary routines are implied | |
| 887 C by this history) | |
| 888 C 891228 Bug was found and repaired inside the DDASSL | |
| 889 C and DDAINI routines. DDAINI was incorrectly | |
| 890 C returning the initial T with Y and YPRIME | |
| 891 C computed at T+H. The routine now returns T+H | |
| 892 C rather than the initial T. | |
| 893 C Cosmetic changes made to DDASTP. | |
| 894 C 900904 Three modifications were made to fix a bug (inside | |
| 895 C DDASSL) re interpolation for continuation calls and | |
| 896 C cases where TN is very close to TSTOP: | |
| 897 C | |
| 898 C 1) In testing for whether H is too large, just | |
| 899 C compare H to (TSTOP - TN), rather than | |
| 900 C (TSTOP - TN) * (1-4*UROUND), and set H to | |
| 901 C TSTOP - TN. This will force DDASTP to step | |
| 902 C exactly to TSTOP under certain situations | |
| 903 C (i.e. when H returned from DDASTP would otherwise | |
| 904 C take TN beyond TSTOP). | |
| 905 C | |
| 906 C 2) Inside the DDASTP loop, interpolate exactly to | |
| 907 C TSTOP if TN is very close to TSTOP (rather than | |
| 908 C interpolating to within roundoff of TSTOP). | |
| 909 C | |
| 910 C 3) Modified IDID description for IDID = 2 to say that | |
| 911 C the solution is returned by stepping exactly to | |
| 912 C TSTOP, rather than TOUT. (In some cases the | |
| 913 C solution is actually obtained by extrapolating | |
| 914 C over a distance near unit roundoff to TSTOP, | |
| 915 C but this small distance is deemed acceptable in | |
| 916 C these circumstances.) | |
| 917 C 901026 Added explicit declarations for all variables and minor | |
| 918 C cosmetic changes to prologue, removed unreferenced labels, | |
| 919 C and improved XERMSG calls. (FNF) | |
| 920 C 901030 Added ERROR MESSAGES section and reworked other sections to | |
| 921 C be of more uniform format. (FNF) | |
| 922 C 910624 Fixed minor bug related to HMAX (five lines ending in | |
| 923 C statement 526 in DDASSL). (LRP) | |
| 924 C | |
| 925 C***END PROLOGUE DDASSL | |
| 926 C | |
| 927 C**End | |
| 928 C | |
| 929 C Declare arguments. | |
| 930 C | |
| 931 INTEGER NEQ, INFO(15), IDID, LRW, IWORK(*), LIW, IPAR(*) | |
| 932 DOUBLE PRECISION | |
| 933 * T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), RWORK(*), | |
| 934 * RPAR(*) | |
| 935 EXTERNAL RES, JAC | |
| 936 C | |
| 937 C Declare externals. | |
| 938 C | |
| 939 EXTERNAL D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, XERMSG | |
| 940 DOUBLE PRECISION D1MACH, DDANRM | |
| 941 C | |
| 942 C Declare local variables. | |
| 943 C | |
| 944 INTEGER I, ITEMP, LALPHA, LBETA, LCJ, LCJOLD, LCTF, LDELTA, | |
| 4429 | 945 * LENIW, LENPD, LENRW, LE, LETF, LGAMMA, LH, LHMAX, LHOLD, |
| 946 * LMXSTP, LIPVT, | |
| 2329 | 947 * LJCALC, LK, LKOLD, LIWM, LML, LMTYPE, LMU, LMXORD, LNJE, LNPD, |
| 948 * LNRE, LNS, LNST, LNSTL, LPD, LPHASE, LPHI, LPSI, LROUND, LS, | |
| 949 * LSIGMA, LTN, LTSTOP, LWM, LWT, MBAND, MSAVE, MXORD, NPD, NTEMP, | |
| 950 * NZFLG | |
| 951 DOUBLE PRECISION | |
| 952 * ATOLI, H, HMAX, HMIN, HO, R, RH, RTOLI, TDIST, TN, TNEXT, | |
| 953 * TSTOP, UROUND, YPNORM | |
| 954 LOGICAL DONE | |
| 955 C Auxiliary variables for conversion of values to be included in | |
| 956 C error messages. | |
| 957 CHARACTER*8 XERN1, XERN2 | |
| 958 CHARACTER*16 XERN3, XERN4 | |
| 959 C | |
| 960 C SET POINTERS INTO IWORK | |
| 961 PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11, | |
| 4429 | 962 * LNRE=12, LNJE=13, LETF=14, LCTF=15, LNPD=16, LMXSTP=21, |
| 963 * LIPVT=22, LJCALC=5, LPHASE=6, LK=7, LKOLD=8, | |
| 2329 | 964 * LNS=9, LNSTL=10, LIWM=1) |
| 965 C | |
| 966 C SET RELATIVE OFFSET INTO RWORK | |
| 967 PARAMETER (NPD=1) | |
| 968 C | |
| 969 C SET POINTERS INTO RWORK | |
| 970 PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, | |
| 971 * LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9, | |
| 972 * LALPHA=11, LBETA=17, LGAMMA=23, | |
| 973 * LPSI=29, LSIGMA=35, LDELTA=41) | |
| 974 C | |
| 975 C***FIRST EXECUTABLE STATEMENT DDASSL | |
| 976 IF(INFO(1).NE.0)GO TO 100 | |
| 977 C | |
| 978 C----------------------------------------------------------------------- | |
| 979 C THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY. | |
| 980 C IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS. | |
| 981 C----------------------------------------------------------------------- | |
| 982 C | |
| 983 C FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO | |
| 984 C ARE EITHER ZERO OR ONE. | |
| 985 DO 10 I=2,11 | |
| 986 IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701 | |
| 987 10 CONTINUE | |
| 988 C | |
| 989 IF(NEQ.LE.0)GO TO 702 | |
| 990 C | |
| 991 C CHECK AND COMPUTE MAXIMUM ORDER | |
| 992 MXORD=5 | |
| 993 IF(INFO(9).EQ.0)GO TO 20 | |
| 994 MXORD=IWORK(LMXORD) | |
| 995 IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703 | |
| 996 20 IWORK(LMXORD)=MXORD | |
| 997 C | |
| 998 C COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU. | |
| 999 IF(INFO(6).NE.0)GO TO 40 | |
| 1000 LENPD=NEQ**2 | |
| 1001 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD | |
| 1002 IF(INFO(5).NE.0)GO TO 30 | |
| 1003 IWORK(LMTYPE)=2 | |
| 1004 GO TO 60 | |
| 1005 30 IWORK(LMTYPE)=1 | |
| 1006 GO TO 60 | |
| 1007 40 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 | |
| 1008 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 | |
| 1009 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ | |
| 1010 IF(INFO(5).NE.0)GO TO 50 | |
| 1011 IWORK(LMTYPE)=5 | |
| 1012 MBAND=IWORK(LML)+IWORK(LMU)+1 | |
| 1013 MSAVE=(NEQ/MBAND)+1 | |
| 1014 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE | |
| 1015 GO TO 60 | |
| 1016 50 IWORK(LMTYPE)=4 | |
| 1017 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD | |
| 1018 C | |
| 1019 C CHECK LENGTHS OF RWORK AND IWORK | |
| 4429 | 1020 60 LENIW=21+NEQ |
| 2329 | 1021 IWORK(LNPD)=LENPD |
| 1022 IF(LRW.LT.LENRW)GO TO 704 | |
| 1023 IF(LIW.LT.LENIW)GO TO 705 | |
| 1024 C | |
| 1025 C CHECK TO SEE THAT TOUT IS DIFFERENT FROM T | |
| 1026 IF(TOUT .EQ. T)GO TO 719 | |
| 1027 C | |
| 1028 C CHECK HMAX | |
| 1029 IF(INFO(7).EQ.0)GO TO 70 | |
| 1030 HMAX=RWORK(LHMAX) | |
| 1031 IF(HMAX.LE.0.0D0)GO TO 710 | |
| 1032 70 CONTINUE | |
| 1033 C | |
| 4429 | 1034 C CHECK AND COMPUTE MAXIMUM STEPS |
| 1035 MXSTP=500 | |
| 1036 IF(INFO(12).EQ.0)GO TO 80 | |
| 1037 MXSTP=IWORK(LMXSTP) | |
| 1038 IF(MXSTP.LT.0)GO TO 716 | |
| 1039 80 IWORK(LMXSTP)=MXSTP | |
| 1040 C | |
| 2329 | 1041 C INITIALIZE COUNTERS |
| 1042 IWORK(LNST)=0 | |
| 1043 IWORK(LNRE)=0 | |
| 1044 IWORK(LNJE)=0 | |
| 1045 C | |
| 1046 IWORK(LNSTL)=0 | |
| 1047 IDID=1 | |
| 1048 GO TO 200 | |
| 1049 C | |
| 1050 C----------------------------------------------------------------------- | |
| 1051 C THIS BLOCK IS FOR CONTINUATION CALLS | |
| 1052 C ONLY. HERE WE CHECK INFO(1),AND IF THE | |
| 1053 C LAST STEP WAS INTERRUPTED WE CHECK WHETHER | |
| 1054 C APPROPRIATE ACTION WAS TAKEN. | |
| 1055 C----------------------------------------------------------------------- | |
| 1056 C | |
| 1057 100 CONTINUE | |
| 1058 IF(INFO(1).EQ.1)GO TO 110 | |
| 1059 IF(INFO(1).NE.-1)GO TO 701 | |
| 1060 C | |
| 1061 C IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED | |
| 1062 C BY AN ERROR CONDITION FROM DDASTP,AND | |
| 1063 C APPROPRIATE ACTION WAS NOT TAKEN. THIS | |
| 1064 C IS A FATAL ERROR. | |
| 1065 WRITE (XERN1, '(I8)') IDID | |
| 1066 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1067 * 'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' // | |
| 1068 * XERN1 // ' AND NO APPROPRIATE ACTION WAS TAKEN. ' // | |
| 1069 * 'RUN TERMINATED', -998, 2) | |
| 1070 RETURN | |
| 1071 110 CONTINUE | |
| 1072 IWORK(LNSTL)=IWORK(LNST) | |
| 1073 C | |
| 1074 C----------------------------------------------------------------------- | |
| 1075 C THIS BLOCK IS EXECUTED ON ALL CALLS. | |
| 1076 C THE ERROR TOLERANCE PARAMETERS ARE | |
| 1077 C CHECKED, AND THE WORK ARRAY POINTERS | |
| 1078 C ARE SET. | |
| 1079 C----------------------------------------------------------------------- | |
| 1080 C | |
| 1081 200 CONTINUE | |
| 1082 C CHECK RTOL,ATOL | |
| 1083 NZFLG=0 | |
| 1084 RTOLI=RTOL(1) | |
| 1085 ATOLI=ATOL(1) | |
| 1086 DO 210 I=1,NEQ | |
| 1087 IF(INFO(2).EQ.1)RTOLI=RTOL(I) | |
| 1088 IF(INFO(2).EQ.1)ATOLI=ATOL(I) | |
| 1089 IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1 | |
| 1090 IF(RTOLI.LT.0.0D0)GO TO 706 | |
| 1091 IF(ATOLI.LT.0.0D0)GO TO 707 | |
| 1092 210 CONTINUE | |
| 1093 IF(NZFLG.EQ.0)GO TO 708 | |
| 1094 C | |
| 1095 C SET UP RWORK STORAGE.IWORK STORAGE IS FIXED | |
| 1096 C IN DATA STATEMENT. | |
| 1097 LE=LDELTA+NEQ | |
| 1098 LWT=LE+NEQ | |
| 1099 LPHI=LWT+NEQ | |
| 1100 LPD=LPHI+(IWORK(LMXORD)+1)*NEQ | |
| 1101 LWM=LPD | |
| 1102 NTEMP=NPD+IWORK(LNPD) | |
| 1103 IF(INFO(1).EQ.1)GO TO 400 | |
| 1104 C | |
| 1105 C----------------------------------------------------------------------- | |
| 1106 C THIS BLOCK IS EXECUTED ON THE INITIAL CALL | |
| 1107 C ONLY. SET THE INITIAL STEP SIZE, AND | |
| 1108 C THE ERROR WEIGHT VECTOR, AND PHI. | |
| 1109 C COMPUTE INITIAL YPRIME, IF NECESSARY. | |
| 1110 C----------------------------------------------------------------------- | |
| 1111 C | |
| 1112 TN=T | |
| 1113 IDID=1 | |
| 1114 C | |
| 1115 C SET ERROR WEIGHT VECTOR WT | |
| 1116 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) | |
| 1117 DO 305 I = 1,NEQ | |
| 1118 IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713 | |
| 1119 305 CONTINUE | |
| 1120 C | |
| 1121 C COMPUTE UNIT ROUNDOFF AND HMIN | |
| 1122 UROUND = D1MACH(4) | |
| 1123 RWORK(LROUND) = UROUND | |
| 1124 HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT)) | |
| 1125 C | |
| 1126 C CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH | |
| 1127 TDIST = ABS(TOUT - T) | |
| 1128 IF(TDIST .LT. HMIN) GO TO 714 | |
| 1129 C | |
| 1130 C CHECK HO, IF THIS WAS INPUT | |
| 1131 IF (INFO(8) .EQ. 0) GO TO 310 | |
| 1132 HO = RWORK(LH) | |
| 1133 IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711 | |
| 1134 IF (HO .EQ. 0.0D0) GO TO 712 | |
| 1135 GO TO 320 | |
| 1136 310 CONTINUE | |
| 1137 C | |
| 1138 C COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER | |
| 1139 C DDASTP OR DDAINI, DEPENDING ON INFO(11) | |
| 1140 HO = 0.001D0*TDIST | |
| 1141 YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR) | |
| 1142 IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM | |
| 1143 HO = SIGN(HO,TOUT-T) | |
| 1144 C ADJUST HO IF NECESSARY TO MEET HMAX BOUND | |
| 1145 320 IF (INFO(7) .EQ. 0) GO TO 330 | |
| 1146 RH = ABS(HO)/RWORK(LHMAX) | |
| 1147 IF (RH .GT. 1.0D0) HO = HO/RH | |
| 1148 C COMPUTE TSTOP, IF APPLICABLE | |
| 1149 330 IF (INFO(4) .EQ. 0) GO TO 340 | |
| 1150 TSTOP = RWORK(LTSTOP) | |
| 1151 IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715 | |
| 1152 IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T | |
| 1153 IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709 | |
| 1154 C | |
| 1155 C COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE | |
| 1156 340 IF (INFO(11) .EQ. 0) GO TO 350 | |
| 1157 CALL DDAINI(TN,Y,YPRIME,NEQ, | |
| 1158 * RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR, | |
| 1159 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), | |
| 1160 * RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND), | |
| 1161 * INFO(10),NTEMP) | |
| 1162 IF (IDID .LT. 0) GO TO 390 | |
| 1163 C | |
| 1164 C LOAD H WITH HO. STORE H IN RWORK(LH) | |
| 1165 350 H = HO | |
| 1166 RWORK(LH) = H | |
| 1167 C | |
| 1168 C LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2) | |
| 1169 ITEMP = LPHI + NEQ | |
| 1170 DO 370 I = 1,NEQ | |
| 1171 RWORK(LPHI + I - 1) = Y(I) | |
| 1172 370 RWORK(ITEMP + I - 1) = H*YPRIME(I) | |
| 1173 C | |
| 1174 390 GO TO 500 | |
| 1175 C | |
| 1176 C------------------------------------------------------- | |
| 1177 C THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS | |
| 1178 C PURPOSE IS TO CHECK STOP CONDITIONS BEFORE | |
| 1179 C TAKING A STEP. | |
| 1180 C ADJUST H IF NECESSARY TO MEET HMAX BOUND | |
| 1181 C------------------------------------------------------- | |
| 1182 C | |
| 1183 400 CONTINUE | |
| 1184 UROUND=RWORK(LROUND) | |
| 1185 DONE = .FALSE. | |
| 1186 TN=RWORK(LTN) | |
| 1187 H=RWORK(LH) | |
| 1188 IF(INFO(7) .EQ. 0) GO TO 410 | |
| 1189 RH = ABS(H)/RWORK(LHMAX) | |
| 1190 IF(RH .GT. 1.0D0) H = H/RH | |
| 1191 410 CONTINUE | |
| 1192 IF(T .EQ. TOUT) GO TO 719 | |
| 1193 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 | |
| 1194 IF(INFO(4) .EQ. 1) GO TO 430 | |
| 1195 IF(INFO(3) .EQ. 1) GO TO 420 | |
| 1196 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 | |
| 1197 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1198 * RWORK(LPHI),RWORK(LPSI)) | |
| 1199 T=TOUT | |
| 1200 IDID = 3 | |
| 1201 DONE = .TRUE. | |
| 1202 GO TO 490 | |
| 1203 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 | |
| 1204 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 | |
| 1205 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1206 * RWORK(LPHI),RWORK(LPSI)) | |
| 1207 T = TN | |
| 1208 IDID = 1 | |
| 1209 DONE = .TRUE. | |
| 1210 GO TO 490 | |
| 1211 425 CONTINUE | |
| 1212 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1213 * RWORK(LPHI),RWORK(LPSI)) | |
| 1214 T = TOUT | |
| 1215 IDID = 3 | |
| 1216 DONE = .TRUE. | |
| 1217 GO TO 490 | |
| 1218 430 IF(INFO(3) .EQ. 1) GO TO 440 | |
| 1219 TSTOP=RWORK(LTSTOP) | |
| 1220 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 | |
| 1221 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 | |
| 1222 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 | |
| 1223 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1224 * RWORK(LPHI),RWORK(LPSI)) | |
| 1225 T=TOUT | |
| 1226 IDID = 3 | |
| 1227 DONE = .TRUE. | |
| 1228 GO TO 490 | |
| 1229 440 TSTOP = RWORK(LTSTOP) | |
| 1230 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 | |
| 1231 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 | |
| 1232 IF((TN-T)*H .LE. 0.0D0) GO TO 450 | |
| 1233 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 | |
| 1234 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1235 * RWORK(LPHI),RWORK(LPSI)) | |
| 1236 T = TN | |
| 1237 IDID = 1 | |
| 1238 DONE = .TRUE. | |
| 1239 GO TO 490 | |
| 1240 445 CONTINUE | |
| 1241 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1242 * RWORK(LPHI),RWORK(LPSI)) | |
| 1243 T = TOUT | |
| 1244 IDID = 3 | |
| 1245 DONE = .TRUE. | |
| 1246 GO TO 490 | |
| 1247 450 CONTINUE | |
| 1248 C CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP | |
| 1249 IF(ABS(TN-TSTOP).GT.100.0D0*UROUND* | |
| 1250 * (ABS(TN)+ABS(H)))GO TO 460 | |
| 1251 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1252 * RWORK(LPHI),RWORK(LPSI)) | |
| 1253 IDID=2 | |
| 1254 T=TSTOP | |
| 1255 DONE = .TRUE. | |
| 1256 GO TO 490 | |
| 1257 460 TNEXT=TN+H | |
| 1258 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 | |
| 1259 H=TSTOP-TN | |
| 1260 RWORK(LH)=H | |
| 1261 C | |
| 1262 490 IF (DONE) GO TO 580 | |
| 1263 C | |
| 1264 C------------------------------------------------------- | |
| 1265 C THE NEXT BLOCK CONTAINS THE CALL TO THE | |
| 1266 C ONE-STEP INTEGRATOR DDASTP. | |
| 1267 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS. | |
| 1268 C CHECK FOR TOO MANY STEPS. | |
| 1269 C UPDATE WT. | |
| 1270 C CHECK FOR TOO MUCH ACCURACY REQUESTED. | |
| 1271 C COMPUTE MINIMUM STEPSIZE. | |
| 1272 C------------------------------------------------------- | |
| 1273 C | |
| 1274 500 CONTINUE | |
| 1275 C CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME | |
| 1276 IF (IDID .EQ. -12) GO TO 527 | |
| 1277 C | |
| 1278 C CHECK FOR TOO MANY STEPS | |
| 4429 | 1279 IF((IWORK(LNST)-IWORK(LNSTL)).LT.IWORK(LMXSTP)) |
| 2329 | 1280 * GO TO 510 |
| 1281 IDID=-1 | |
| 1282 GO TO 527 | |
| 1283 C | |
| 1284 C UPDATE WT | |
| 1285 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI), | |
| 1286 * RWORK(LWT),RPAR,IPAR) | |
| 1287 DO 520 I=1,NEQ | |
| 1288 IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520 | |
| 1289 IDID=-3 | |
| 1290 GO TO 527 | |
| 1291 520 CONTINUE | |
| 1292 C | |
| 1293 C TEST FOR TOO MUCH ACCURACY REQUESTED. | |
| 1294 R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)* | |
| 1295 * 100.0D0*UROUND | |
| 1296 IF(R.LE.1.0D0)GO TO 525 | |
| 1297 C MULTIPLY RTOL AND ATOL BY R AND RETURN | |
| 1298 IF(INFO(2).EQ.1)GO TO 523 | |
| 1299 RTOL(1)=R*RTOL(1) | |
| 1300 ATOL(1)=R*ATOL(1) | |
| 1301 IDID=-2 | |
| 1302 GO TO 527 | |
| 1303 523 DO 524 I=1,NEQ | |
| 1304 RTOL(I)=R*RTOL(I) | |
| 1305 524 ATOL(I)=R*ATOL(I) | |
| 1306 IDID=-2 | |
| 1307 GO TO 527 | |
| 1308 525 CONTINUE | |
| 1309 C | |
| 1310 C COMPUTE MINIMUM STEPSIZE | |
| 1311 HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT)) | |
| 1312 C | |
| 1313 C TEST H VS. HMAX | |
| 1314 IF (INFO(7) .EQ. 0) GO TO 526 | |
| 1315 RH = ABS(H)/RWORK(LHMAX) | |
| 1316 IF (RH .GT. 1.0D0) H = H/RH | |
|
19627
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
15271
diff
changeset
|
1317 526 CONTINUE |
| 2329 | 1318 C |
| 1319 CALL DDASTP(TN,Y,YPRIME,NEQ, | |
| 1320 * RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR, | |
| 1321 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), | |
| 1322 * RWORK(LWM),IWORK(LIWM), | |
| 1323 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), | |
| 1324 * RWORK(LPSI),RWORK(LSIGMA), | |
| 1325 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD), | |
| 1326 * RWORK(LS),HMIN,RWORK(LROUND), | |
| 1327 * IWORK(LPHASE),IWORK(LJCALC),IWORK(LK), | |
| 1328 * IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP) | |
| 1329 527 IF(IDID.LT.0)GO TO 600 | |
| 1330 C | |
| 1331 C-------------------------------------------------------- | |
| 1332 C THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN | |
| 1333 C FROM DDASTP (IDID=1). TEST FOR STOP CONDITIONS. | |
| 1334 C-------------------------------------------------------- | |
| 1335 C | |
| 1336 IF(INFO(4).NE.0)GO TO 540 | |
| 1337 IF(INFO(3).NE.0)GO TO 530 | |
| 1338 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 | |
| 1339 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1340 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1341 IDID=3 | |
| 1342 T=TOUT | |
| 1343 GO TO 580 | |
| 1344 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 | |
| 1345 T=TN | |
| 1346 IDID=1 | |
| 1347 GO TO 580 | |
| 1348 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1349 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1350 IDID=3 | |
| 1351 T=TOUT | |
| 1352 GO TO 580 | |
| 1353 540 IF(INFO(3).NE.0)GO TO 550 | |
| 1354 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 | |
| 1355 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1356 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1357 T=TOUT | |
| 1358 IDID=3 | |
| 1359 GO TO 580 | |
| 1360 542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND* | |
| 1361 * (ABS(TN)+ABS(H)))GO TO 545 | |
| 1362 TNEXT=TN+H | |
| 1363 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 | |
| 1364 H=TSTOP-TN | |
| 1365 GO TO 500 | |
| 1366 545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, | |
| 1367 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1368 IDID=2 | |
| 1369 T=TSTOP | |
| 1370 GO TO 580 | |
| 1371 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 | |
| 1372 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552 | |
| 1373 T=TN | |
| 1374 IDID=1 | |
| 1375 GO TO 580 | |
| 1376 552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, | |
| 1377 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1378 IDID=2 | |
| 1379 T=TSTOP | |
| 1380 GO TO 580 | |
| 1381 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1382 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1383 T=TOUT | |
| 1384 IDID=3 | |
| 1385 GO TO 580 | |
| 1386 C | |
| 1387 C-------------------------------------------------------- | |
| 1388 C ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM | |
| 1389 C THIS BLOCK. | |
| 1390 C-------------------------------------------------------- | |
| 1391 C | |
| 1392 580 CONTINUE | |
| 1393 RWORK(LTN)=TN | |
| 1394 RWORK(LH)=H | |
| 1395 RETURN | |
| 1396 C | |
| 1397 C----------------------------------------------------------------------- | |
| 1398 C THIS BLOCK HANDLES ALL UNSUCCESSFUL | |
| 1399 C RETURNS OTHER THAN FOR ILLEGAL INPUT. | |
| 1400 C----------------------------------------------------------------------- | |
| 1401 C | |
| 1402 600 CONTINUE | |
| 1403 ITEMP=-IDID | |
| 1404 GO TO (610,620,630,690,690,640,650,660,670,675, | |
| 1405 * 680,685), ITEMP | |
| 1406 C | |
| 1407 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE | |
| 1408 C REACHING TOUT | |
| 1409 610 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1410 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1411 * 'AT CURRENT T = ' // XERN3 // ' 500 STEPS TAKEN ON THIS ' // | |
| 1412 * 'CALL BEFORE REACHING TOUT', IDID, 1) | |
| 1413 GO TO 690 | |
| 1414 C | |
| 1415 C TOO MUCH ACCURACY FOR MACHINE PRECISION | |
| 1416 620 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1417 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1418 * 'AT T = ' // XERN3 // ' TOO MUCH ACCURACY REQUESTED FOR ' // | |
| 1419 * 'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' // | |
| 1420 * 'APPROPRIATE VALUES', IDID, 1) | |
| 1421 GO TO 690 | |
| 1422 C | |
| 1423 C WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM) | |
| 1424 630 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1425 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1426 * 'AT T = ' // XERN3 // ' SOME ELEMENT OF WT HAS BECOME .LE. ' // | |
| 1427 * '0.0', IDID, 1) | |
| 1428 GO TO 690 | |
| 1429 C | |
| 1430 C ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN | |
| 1431 640 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1432 WRITE (XERN4, '(1P,D15.6)') H | |
| 1433 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1434 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1435 * ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN', | |
| 1436 * IDID, 1) | |
| 1437 GO TO 690 | |
| 1438 C | |
| 1439 C CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN | |
| 1440 650 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1441 WRITE (XERN4, '(1P,D15.6)') H | |
| 1442 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1443 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1444 * ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' // | |
| 1445 * 'ABS(H)=HMIN', IDID, 1) | |
| 1446 GO TO 690 | |
| 1447 C | |
| 1448 C THE ITERATION MATRIX IS SINGULAR | |
| 1449 660 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1450 WRITE (XERN4, '(1P,D15.6)') H | |
| 1451 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1452 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1453 * ' THE ITERATION MATRIX IS SINGULAR', IDID, 1) | |
| 1454 GO TO 690 | |
| 1455 C | |
| 1456 C CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES. | |
| 1457 670 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1458 WRITE (XERN4, '(1P,D15.6)') H | |
| 1459 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1460 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1461 * ' THE CORRECTOR COULD NOT CONVERGE. ALSO, THE ERROR TEST ' // | |
| 1462 * 'FAILED REPEATEDLY.', IDID, 1) | |
| 1463 GO TO 690 | |
| 1464 C | |
| 1465 C CORRECTOR FAILURE BECAUSE IRES = -1 | |
| 1466 675 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1467 WRITE (XERN4, '(1P,D15.6)') H | |
| 1468 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1469 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1470 * ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' // | |
| 1471 * 'TO MINUS ONE', IDID, 1) | |
| 1472 GO TO 690 | |
| 1473 C | |
| 1474 C FAILURE BECAUSE IRES = -2 | |
| 1475 680 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1476 WRITE (XERN4, '(1P,D15.6)') H | |
| 1477 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1478 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1479 * ' IRES WAS EQUAL TO MINUS TWO', IDID, 1) | |
| 1480 GO TO 690 | |
| 1481 C | |
| 1482 C FAILED TO COMPUTE INITIAL YPRIME | |
| 1483 685 WRITE (XERN3, '(1P,D15.6)') TN | |
| 1484 WRITE (XERN4, '(1P,D15.6)') HO | |
| 1485 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1486 * 'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 // | |
| 1487 * ' THE INITIAL YPRIME COULD NOT BE COMPUTED', IDID, 1) | |
| 1488 GO TO 690 | |
| 1489 C | |
| 1490 690 CONTINUE | |
| 1491 INFO(1)=-1 | |
| 1492 T=TN | |
| 1493 RWORK(LTN)=TN | |
| 1494 RWORK(LH)=H | |
| 1495 RETURN | |
| 1496 C | |
| 1497 C----------------------------------------------------------------------- | |
| 1498 C THIS BLOCK HANDLES ALL ERROR RETURNS DUE | |
| 1499 C TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING | |
| 1500 C DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS | |
| 1501 C CALLED. IF THIS HAPPENS TWICE IN | |
| 1502 C SUCCESSION, EXECUTION IS TERMINATED | |
| 1503 C | |
| 1504 C----------------------------------------------------------------------- | |
| 1505 701 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1506 * 'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1) | |
| 1507 GO TO 750 | |
| 1508 C | |
| 1509 702 WRITE (XERN1, '(I8)') NEQ | |
| 1510 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1511 * 'NEQ = ' // XERN1 // ' .LE. 0', 2, 1) | |
| 1512 GO TO 750 | |
| 1513 C | |
| 1514 703 WRITE (XERN1, '(I8)') MXORD | |
| 1515 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1516 * 'MAXORD = ' // XERN1 // ' NOT IN RANGE', 3, 1) | |
| 1517 GO TO 750 | |
| 1518 C | |
| 1519 704 WRITE (XERN1, '(I8)') LENRW | |
| 1520 WRITE (XERN2, '(I8)') LRW | |
| 1521 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1522 * 'RWORK LENGTH NEEDED, LENRW = ' // XERN1 // | |
| 1523 * ', EXCEEDS LRW = ' // XERN2, 4, 1) | |
| 1524 GO TO 750 | |
| 1525 C | |
| 1526 705 WRITE (XERN1, '(I8)') LENIW | |
| 1527 WRITE (XERN2, '(I8)') LIW | |
| 1528 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1529 * 'IWORK LENGTH NEEDED, LENIW = ' // XERN1 // | |
| 1530 * ', EXCEEDS LIW = ' // XERN2, 5, 1) | |
| 1531 GO TO 750 | |
| 1532 C | |
| 1533 706 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1534 * 'SOME ELEMENT OF RTOL IS .LT. 0', 6, 1) | |
| 1535 GO TO 750 | |
| 1536 C | |
| 1537 707 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1538 * 'SOME ELEMENT OF ATOL IS .LT. 0', 7, 1) | |
| 1539 GO TO 750 | |
| 1540 C | |
| 1541 708 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1542 * 'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1) | |
| 1543 GO TO 750 | |
| 1544 C | |
| 1545 709 WRITE (XERN3, '(1P,D15.6)') TSTOP | |
| 1546 WRITE (XERN4, '(1P,D15.6)') TOUT | |
| 1547 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1548 * 'INFO(4) = 1 AND TSTOP = ' // XERN3 // ' BEHIND TOUT = ' // | |
| 1549 * XERN4, 9, 1) | |
| 1550 GO TO 750 | |
| 1551 C | |
| 1552 710 WRITE (XERN3, '(1P,D15.6)') HMAX | |
| 1553 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1554 * 'HMAX = ' // XERN3 // ' .LT. 0.0', 10, 1) | |
| 1555 GO TO 750 | |
| 1556 C | |
| 1557 711 WRITE (XERN3, '(1P,D15.6)') TOUT | |
| 1558 WRITE (XERN4, '(1P,D15.6)') T | |
| 1559 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1560 * 'TOUT = ' // XERN3 // ' BEHIND T = ' // XERN4, 11, 1) | |
| 1561 GO TO 750 | |
| 1562 C | |
| 1563 712 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1564 * 'INFO(8)=1 AND H0=0.0', 12, 1) | |
| 1565 GO TO 750 | |
| 1566 C | |
| 1567 713 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1568 * 'SOME ELEMENT OF WT IS .LE. 0.0', 13, 1) | |
| 1569 GO TO 750 | |
| 1570 C | |
| 1571 714 WRITE (XERN3, '(1P,D15.6)') TOUT | |
| 1572 WRITE (XERN4, '(1P,D15.6)') T | |
| 1573 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1574 * 'TOUT = ' // XERN3 // ' TOO CLOSE TO T = ' // XERN4 // | |
| 1575 * ' TO START INTEGRATION', 14, 1) | |
| 1576 GO TO 750 | |
| 1577 C | |
| 1578 715 WRITE (XERN3, '(1P,D15.6)') TSTOP | |
| 1579 WRITE (XERN4, '(1P,D15.6)') T | |
| 1580 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1581 * 'INFO(4)=1 AND TSTOP = ' // XERN3 // ' BEHIND T = ' // XERN4, | |
| 1582 * 15, 1) | |
| 1583 GO TO 750 | |
| 1584 C | |
| 4429 | 1585 716 WRITE (XERN1, '(I8)') MXSTP |
| 1586 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1587 * 'INFO(12)=1 AND MXSTP = ' // XERN1 // ' ILLEGAL.', 3, 1) | |
| 1588 GO TO 750 | |
| 1589 C | |
| 2329 | 1590 717 WRITE (XERN1, '(I8)') IWORK(LML) |
| 1591 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1592 * 'ML = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ', | |
| 1593 * 17, 1) | |
| 1594 GO TO 750 | |
| 1595 C | |
| 1596 718 WRITE (XERN1, '(I8)') IWORK(LMU) | |
| 1597 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1598 * 'MU = ' // XERN1 // ' ILLEGAL. EITHER .LT. 0 OR .GT. NEQ', | |
| 1599 * 18, 1) | |
| 1600 GO TO 750 | |
| 1601 C | |
| 1602 719 WRITE (XERN3, '(1P,D15.6)') TOUT | |
| 1603 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1604 * 'TOUT = T = ' // XERN3, 19, 1) | |
| 1605 GO TO 750 | |
| 1606 C | |
| 1607 750 IDID=-33 | |
| 1608 IF(INFO(1).EQ.-1) THEN | |
| 1609 CALL XERMSG ('SLATEC', 'DDASSL', | |
| 1610 * 'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' // | |
| 1611 * 'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2) | |
| 1612 ENDIF | |
| 1613 C | |
| 1614 INFO(1)=-1 | |
| 1615 RETURN | |
| 1616 C-----------END OF SUBROUTINE DDASSL------------------------------------ | |
| 1617 END |