Mercurial > octave-nkf
annotate liboctave/cruft/dasrt/ddasrt.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 |
|---|---|
| 4000 | 1 SUBROUTINE DDASRT (RES,NEQ,T,Y,YPRIME,TOUT, |
| 2 * INFO,RTOL,ATOL,IDID,RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC, | |
| 3 * G,NG,JROOT) | |
| 4 C | |
| 5 C***BEGIN PROLOGUE DDASRT | |
| 6 C***DATE WRITTEN 821001 (YYMMDD) | |
| 7 C***REVISION DATE 910624 (YYMMDD) | |
| 8 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC,BACKWARD DIFFERENTIATION FORMULAS | |
| 9 C IMPLICIT DIFFERENTIAL SYSTEMS | |
| 10 C***AUTHOR PETZOLD,LINDA R.,COMPUTING AND MATHEMATICS RESEARCH DIVISION | |
| 11 C LAWRENCE LIVERMORE NATIONAL LABORATORY | |
| 12 C L - 316, P.O. Box 808, | |
| 13 C LIVERMORE, CA. 94550 | |
| 14 C***PURPOSE This code solves a system of differential/algebraic | |
| 15 C equations of the form F(T,Y,YPRIME) = 0. | |
| 16 C***DESCRIPTION | |
| 17 C | |
| 18 C *Usage: | |
| 19 C | |
| 20 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) | |
| 21 C EXTERNAL RES, JAC, G | |
| 22 C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR, NG, | |
| 23 C * JROOT(NG) | |
| 24 C DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL, | |
| 25 C * RWORK(LRW), RPAR | |
| 26 C | |
| 27 C CALL DDASRT (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, | |
| 28 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC) | |
| 29 C | |
| 30 C | |
| 31 C | |
| 32 C *Arguments: | |
| 33 C | |
| 34 C RES:EXT This is a subroutine which you provide to define the | |
| 35 C differential/algebraic system. | |
| 36 C | |
| 37 C NEQ:IN This is the number of equations to be solved. | |
| 38 C | |
| 39 C T:INOUT This is the current value of the independent variable. | |
| 40 C | |
| 41 C Y(*):INOUT This array contains the solution components at T. | |
| 42 C | |
| 43 C YPRIME(*):INOUT This array contains the derivatives of the solution | |
| 44 C components at T. | |
| 45 C | |
| 46 C TOUT:IN This is a point at which a solution is desired. | |
| 47 C | |
| 48 C INFO(N):IN The basic task of the code is to solve the system from T | |
| 49 C to TOUT and return an answer at TOUT. INFO is an integer | |
| 50 C array which is used to communicate exactly how you want | |
| 51 C this task to be carried out. N must be greater than or | |
| 52 C equal to 15. | |
| 53 C | |
| 54 C RTOL,ATOL:INOUT These quantities represent absolute and relative | |
| 55 C error tolerances which you provide to indicate how | |
| 56 C accurately you wish the solution to be computed. | |
| 57 C You may choose them to be both scalars or else | |
| 58 C both vectors. | |
| 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 decide | |
| 62 C 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. | |
| 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. | |
| 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 to | |
| 79 C provide for defining a matrix of partial derivatives | |
| 80 C described below. | |
| 81 C | |
| 82 C G This is the name of the subroutine for defining | |
| 83 C constraint functions, G(T,Y), whose roots are desired | |
| 84 C during the integration. This name must be declared | |
| 85 C external in the calling program. | |
| 86 C | |
| 87 C NG This is the number of constraint functions G(I). | |
| 88 C If there are none, set NG=0, and pass a dummy name | |
| 89 C for G. | |
| 90 C | |
| 91 C JROOT This is an integer array of length NG for output | |
| 92 C of root information. | |
| 93 C | |
| 94 C | |
| 95 C *Description | |
| 96 C | |
| 97 C QUANTITIES WHICH MAY BE ALTERED BY THE CODE ARE | |
| 98 C T,Y(*),YPRIME(*),INFO(1),RTOL,ATOL, | |
| 99 C IDID,RWORK(*) AND IWORK(*). | |
| 100 C | |
| 101 C Subroutine DDASRT uses the backward differentiation formulas of | |
| 102 C orders one through five to solve a system of the above form for Y and | |
| 103 C YPRIME. Values for Y and YPRIME at the initial time must be given as | |
| 104 C input. These values must be consistent, (that is, if T,Y,YPRIME are | |
| 105 C the given initial values, they must satisfy F(T,Y,YPRIME) = 0.). The | |
| 106 C subroutine solves the system from T to TOUT. | |
| 107 C It is easy to continue the solution to get results at additional | |
| 108 C TOUT. This is the interval mode of operation. Intermediate results | |
| 109 C can also be obtained easily by using the intermediate-output | |
| 110 C capability. If DDASRT detects a sign-change in G(T,Y), then | |
| 111 C it will return the intermediate value of T and Y for which | |
| 112 C G(T,Y) = 0. | |
| 113 C | |
| 114 C ---------INPUT-WHAT TO DO ON THE FIRST CALL TO DDASRT--------------- | |
| 115 C | |
| 116 C | |
| 117 C The first call of the code is defined to be the start of each new | |
| 118 C problem. Read through the descriptions of all the following items, | |
| 119 C provide sufficient storage space for designated arrays, set | |
| 120 C appropriate variables for the initialization of the problem, and | |
| 121 C give information about how you want the problem to be solved. | |
| 122 C | |
| 123 C | |
| 124 C RES -- Provide a subroutine of the form | |
| 125 C SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) | |
| 126 C to define the system of differential/algebraic | |
| 127 C equations which is to be solved. For the given values | |
| 128 C of T,Y and YPRIME, the subroutine should | |
| 129 C return the residual of the defferential/algebraic | |
| 130 C system | |
| 131 C DELTA = F(T,Y,YPRIME) | |
| 132 C (DELTA(*) is a vector of length NEQ which is | |
| 133 C output for RES.) | |
| 134 C | |
| 135 C Subroutine RES must not alter T,Y or YPRIME. | |
| 136 C You must declare the name RES in an external | |
| 137 C statement in your program that calls DDASRT. | |
| 138 C You must dimension Y,YPRIME and DELTA in RES. | |
| 139 C | |
| 140 C IRES is an integer flag which is always equal to | |
| 141 C zero on input. Subroutine RES should alter IRES | |
| 142 C only if it encounters an illegal value of Y or | |
| 143 C a stop condition. Set IRES = -1 if an input value | |
| 144 C is illegal, and DDASRT will try to solve the problem | |
| 145 C without getting IRES = -1. If IRES = -2, DDASRT | |
| 146 C will return control to the calling program | |
| 147 C with IDID = -11. | |
| 148 C | |
| 149 C RPAR and IPAR are real and integer parameter arrays which | |
| 150 C you can use for communication between your calling program | |
| 151 C and subroutine RES. They are not altered by DDASRT. If you | |
| 152 C do not need RPAR or IPAR, ignore these parameters by treat- | |
| 153 C ing them as dummy arguments. If you do choose to use them, | |
| 154 C dimension them in your calling program and in RES as arrays | |
| 155 C of appropriate length. | |
| 156 C | |
| 157 C NEQ -- Set it to the number of differential equations. | |
| 158 C (NEQ .GE. 1) | |
| 159 C | |
| 160 C T -- Set it to the initial point of the integration. | |
| 161 C T must be defined as a variable. | |
| 162 C | |
| 163 C Y(*) -- Set this vector to the initial values of the NEQ solution | |
| 164 C components at the initial point. You must dimension Y of | |
| 165 C length at least NEQ in your calling program. | |
| 166 C | |
| 167 C YPRIME(*) -- Set this vector to the initial values of | |
| 168 C the NEQ first derivatives of the solution | |
| 169 C components at the initial point. You | |
| 170 C must dimension YPRIME at least NEQ | |
| 171 C in your calling program. If you do not | |
| 172 C know initial values of some of the solution | |
| 173 C components, see the explanation of INFO(11). | |
| 174 C | |
| 175 C TOUT - Set it to the first point at which a solution | |
| 176 C is desired. You can not take TOUT = T. | |
| 177 C integration either forward in T (TOUT .GT. T) or | |
| 178 C backward in T (TOUT .LT. T) is permitted. | |
| 179 C | |
| 180 C The code advances the solution from T to TOUT using | |
| 181 C step sizes which are automatically selected so as to | |
| 182 C achieve the desired accuracy. If you wish, the code will | |
| 183 C return with the solution and its derivative at | |
| 184 C intermediate steps (intermediate-output mode) so that | |
| 185 C you can monitor them, but you still must provide TOUT in | |
| 186 C accord with the basic aim of the code. | |
| 187 C | |
| 188 C the first step taken by the code is a critical one | |
| 189 C because it must reflect how fast the solution changes near | |
| 190 C the initial point. The code automatically selects an | |
| 191 C initial step size which is practically always suitable for | |
| 192 C the problem. By using the fact that the code will not step | |
| 193 C past TOUT in the first step, you could, if necessary, | |
| 194 C restrict the length of the initial step size. | |
| 195 C | |
| 196 C For some problems it may not be permissable to integrate | |
| 197 C past a point TSTOP because a discontinuity occurs there | |
| 198 C or the solution or its derivative is not defined beyond | |
| 199 C TSTOP. When you have declared a TSTOP point (SEE INFO(4) | |
| 200 C and RWORK(1)), you have told the code not to integrate | |
| 201 C past TSTOP. In this case any TOUT beyond TSTOP is invalid | |
| 202 C input. | |
| 203 C | |
| 204 C INFO(*) - Use the INFO array to give the code more details about | |
| 205 C how you want your problem solved. This array should be | |
| 206 C dimensioned of length 15, though DDASRT uses | |
| 207 C only the first twelve entries. You must respond to all of | |
| 208 C the following items which are arranged as questions. The | |
| 209 C simplest use of the code corresponds to answering all | |
| 210 C questions as yes, i.e. setting all entries of INFO to 0. | |
| 211 C | |
| 212 C INFO(1) - This parameter enables the code to initialize | |
| 213 C itself. You must set it to indicate the start of every | |
| 214 C new problem. | |
| 215 C | |
| 216 C **** Is this the first call for this problem ... | |
| 217 C Yes - Set INFO(1) = 0 | |
| 218 C No - Not applicable here. | |
| 219 C See below for continuation calls. **** | |
| 220 C | |
| 221 C INFO(2) - How much accuracy you want of your solution | |
| 222 C is specified by the error tolerances RTOL and ATOL. | |
| 223 C The simplest use is to take them both to be scalars. | |
| 224 C To obtain more flexibility, they can both be vectors. | |
| 225 C The code must be told your choice. | |
| 226 C | |
| 227 C **** Are both error tolerances RTOL, ATOL scalars ... | |
| 228 C Yes - Set INFO(2) = 0 | |
| 229 C and input scalars for both RTOL and ATOL | |
| 230 C No - Set INFO(2) = 1 | |
| 231 C and input arrays for both RTOL and ATOL **** | |
| 232 C | |
| 233 C INFO(3) - The code integrates from T in the direction | |
| 234 C of TOUT by steps. If you wish, it will return the | |
| 235 C computed solution and derivative at the next | |
| 236 C intermediate step (the intermediate-output mode) or | |
| 237 C TOUT, whichever comes first. This is a good way to | |
| 238 C proceed if you want to see the behavior of the solution. | |
| 239 C If you must have solutions at a great many specific | |
| 240 C TOUT points, this code will compute them efficiently. | |
| 241 C | |
| 242 C **** Do you want the solution only at | |
| 243 C TOUT (and not at the next intermediate step) ... | |
| 244 C Yes - Set INFO(3) = 0 | |
| 245 C No - Set INFO(3) = 1 **** | |
| 246 C | |
| 247 C INFO(4) - To handle solutions at a great many specific | |
| 248 C values TOUT efficiently, this code may integrate past | |
| 249 C TOUT and interpolate to obtain the result at TOUT. | |
| 250 C Sometimes it is not possible to integrate beyond some | |
| 251 C point TSTOP because the equation changes there or it is | |
| 252 C not defined past TSTOP. Then you must tell the code | |
| 253 C not to go past. | |
| 254 C | |
| 255 C **** Can the integration be carried out without any | |
| 256 C restrictions on the independent variable T ... | |
| 257 C Yes - Set INFO(4)=0 | |
| 258 C No - Set INFO(4)=1 | |
| 259 C and define the stopping point TSTOP by | |
| 260 C setting RWORK(1)=TSTOP **** | |
| 261 C | |
| 262 C INFO(5) - To solve differential/algebraic problems it is | |
| 263 C necessary to use a matrix of partial derivatives of the | |
| 264 C system of differential equations. If you do not | |
| 265 C provide a subroutine to evaluate it analytically (see | |
| 266 C description of the item JAC in the call list), it will | |
| 267 C be approximated by numerical differencing in this code. | |
| 268 C although it is less trouble for you to have the code | |
| 269 C compute partial derivatives by numerical differencing, | |
| 270 C the solution will be more reliable if you provide the | |
| 271 C derivatives via JAC. Sometimes numerical differencing | |
| 272 C is cheaper than evaluating derivatives in JAC and | |
| 273 C sometimes it is not - this depends on your problem. | |
| 274 C | |
| 275 C **** Do you want the code to evaluate the partial | |
| 276 C derivatives automatically by numerical differences ... | |
| 277 C Yes - Set INFO(5)=0 | |
| 278 C No - Set INFO(5)=1 | |
| 279 C and provide subroutine JAC for evaluating the | |
| 280 C matrix of partial derivatives **** | |
| 281 C | |
| 282 C INFO(6) - DDASRT will perform much better if the matrix of | |
| 283 C partial derivatives, DG/DY + CJ*DG/DYPRIME, | |
| 284 C (here CJ is a scalar determined by DDASRT) | |
| 285 C is banded and the code is told this. In this | |
| 286 C case, the storage needed will be greatly reduced, | |
| 287 C numerical differencing will be performed much cheaper, | |
| 288 C and a number of important algorithms will execute much | |
| 289 C faster. The differential equation is said to have | |
| 290 C half-bandwidths ML (lower) and MU (upper) if equation i | |
| 291 C involves only unknowns Y(J) with | |
| 292 C I-ML .LE. J .LE. I+MU | |
| 293 C for all I=1,2,...,NEQ. Thus, ML and MU are the widths | |
| 294 C of the lower and upper parts of the band, respectively, | |
| 295 C with the main diagonal being excluded. If you do not | |
| 296 C indicate that the equation has a banded matrix of partial | |
| 297 C derivatives, the code works with a full matrix of NEQ**2 | |
| 298 C elements (stored in the conventional way). Computations | |
| 299 C with banded matrices cost less time and storage than with | |
| 300 C full matrices if 2*ML+MU .LT. NEQ. If you tell the | |
| 301 C code that the matrix of partial derivatives has a banded | |
| 302 C structure and you want to provide subroutine JAC to | |
| 303 C compute the partial derivatives, then you must be careful | |
| 304 C to store the elements of the matrix in the special form | |
| 305 C indicated in the description of JAC. | |
| 306 C | |
| 307 C **** Do you want to solve the problem using a full | |
| 308 C (dense) matrix (and not a special banded | |
| 309 C structure) ... | |
| 310 C Yes - Set INFO(6)=0 | |
| 311 C No - Set INFO(6)=1 | |
| 312 C and provide the lower (ML) and upper (MU) | |
| 313 C bandwidths by setting | |
| 314 C IWORK(1)=ML | |
| 315 C IWORK(2)=MU **** | |
| 316 C | |
| 317 C | |
| 318 C INFO(7) -- You can specify a maximum (absolute value of) | |
| 319 C stepsize, so that the code | |
| 320 C will avoid passing over very | |
| 321 C large regions. | |
| 322 C | |
| 323 C **** Do you want the code to decide | |
| 324 C on its own maximum stepsize? | |
| 325 C Yes - Set INFO(7)=0 | |
| 326 C No - Set INFO(7)=1 | |
| 327 C and define HMAX by setting | |
| 328 C RWORK(2)=HMAX **** | |
| 329 C | |
| 330 C INFO(8) -- Differential/algebraic problems | |
| 331 C may occaisionally suffer from | |
| 332 C severe scaling difficulties on the | |
| 333 C first step. If you know a great deal | |
| 334 C about the scaling of your problem, you can | |
| 335 C help to alleviate this problem by | |
| 336 C specifying an initial stepsize H0. | |
| 337 C | |
| 338 C **** Do you want the code to define | |
| 339 C its own initial stepsize? | |
| 340 C Yes - Set INFO(8)=0 | |
| 341 C No - Set INFO(8)=1 | |
| 342 C and define H0 by setting | |
| 343 C RWORK(3)=H0 **** | |
| 344 C | |
| 345 C INFO(9) -- If storage is a severe problem, | |
| 346 C you can save some locations by | |
| 347 C restricting the maximum order MAXORD. | |
| 348 C the default value is 5. for each | |
| 349 C order decrease below 5, the code | |
| 350 C requires NEQ fewer locations, however | |
| 351 C it is likely to be slower. In any | |
| 352 C case, you must have 1 .LE. MAXORD .LE. 5 | |
| 353 C **** Do you want the maximum order to | |
| 354 C default to 5? | |
| 355 C Yes - Set INFO(9)=0 | |
| 356 C No - Set INFO(9)=1 | |
| 357 C and define MAXORD by setting | |
| 358 C IWORK(3)=MAXORD **** | |
| 359 C | |
| 360 C INFO(10) --If you know that the solutions to your equations | |
| 361 C will always be nonnegative, it may help to set this | |
| 362 C parameter. However, it is probably best to | |
| 363 C try the code without using this option first, | |
| 364 C and only to use this option if that doesn't | |
| 365 C work very well. | |
| 366 C **** Do you want the code to solve the problem without | |
| 367 C invoking any special nonnegativity constraints? | |
| 368 C Yes - Set INFO(10)=0 | |
| 369 C No - Set INFO(10)=1 | |
| 370 C | |
| 371 C INFO(11) --DDASRT normally requires the initial T, | |
| 372 C Y, and YPRIME to be consistent. That is, | |
| 373 C you must have F(T,Y,YPRIME) = 0 at the initial | |
| 374 C time. If you do not know the initial | |
| 375 C derivative precisely, you can let DDASRT try | |
| 376 C to compute it. | |
| 377 C **** Are the initial T, Y, YPRIME consistent? | |
| 378 C Yes - Set INFO(11) = 0 | |
| 379 C No - Set INFO(11) = 1, | |
| 380 C and set YPRIME to an initial approximation | |
| 381 C to YPRIME. (If you have no idea what | |
| 382 C YPRIME should be, set it to zero. Note | |
| 383 C that the initial Y should be such | |
| 384 C that there must exist a YPRIME so that | |
| 385 C F(T,Y,YPRIME) = 0.) | |
| 386 C | |
| 387 C INFO(12) --Maximum number of steps. | |
| 388 C **** Do you want to let DDASRT use the default limit for | |
| 389 C the number of steps? | |
| 4429 | 390 C Yes - Set INFO(12) = 0 |
| 391 C No - Set INFO(12) = 1, | |
| 4000 | 392 C and define the maximum number of steps |
| 4428 | 393 C by setting IWORK(21)=MXSTEP |
| 4000 | 394 C |
| 395 C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL | |
| 396 C error tolerances to tell the code how accurately you | |
| 397 C want the solution to be computed. They must be defined | |
| 398 C as variables because the code may change them. You | |
| 399 C have two choices -- | |
| 400 C Both RTOL and ATOL are scalars. (INFO(2)=0) | |
| 401 C Both RTOL and ATOL are vectors. (INFO(2)=1) | |
| 402 C in either case all components must be non-negative. | |
| 403 C | |
| 404 C The tolerances are used by the code in a local error | |
| 405 C test at each step which requires roughly that | |
| 406 C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL | |
| 407 C for each vector component. | |
| 408 C (More specifically, a root-mean-square norm is used to | |
| 409 C measure the size of vectors, and the error test uses the | |
| 410 C magnitude of the solution at the beginning of the step.) | |
| 411 C | |
| 412 C The true (global) error is the difference between the | |
| 413 C true solution of the initial value problem and the | |
| 414 C computed approximation. Practically all present day | |
| 415 C codes, including this one, control the local error at | |
| 416 C each step and do not even attempt to control the global | |
| 417 C error directly. | |
| 418 C Usually, but not always, the true accuracy of the | |
| 419 C computed Y is comparable to the error tolerances. This | |
| 420 C code will usually, but not always, deliver a more | |
| 421 C accurate solution if you reduce the tolerances and | |
| 422 C integrate again. By comparing two such solutions you | |
| 423 C can get a fairly reliable idea of the true error in the | |
| 424 C solution at the bigger tolerances. | |
| 425 C | |
| 426 C Setting ATOL=0. results in a pure relative error test on | |
| 427 C that component. Setting RTOL=0. results in a pure | |
| 428 C absolute error test on that component. A mixed test | |
| 429 C with non-zero RTOL and ATOL corresponds roughly to a | |
| 430 C relative error test when the solution component is much | |
| 431 C bigger than ATOL and to an absolute error test when the | |
| 432 C solution component is smaller than the threshhold ATOL. | |
| 433 C | |
| 434 C The code will not attempt to compute a solution at an | |
| 435 C accuracy unreasonable for the machine being used. It | |
| 436 C will advise you if you ask for too much accuracy and | |
| 437 C inform you as to the maximum accuracy it believes | |
| 438 C possible. | |
| 439 C | |
| 440 C RWORK(*) -- Dimension this real work array of length LRW in your | |
| 441 C calling program. | |
| 442 C | |
| 443 C LRW -- Set it to the declared length of the RWORK array. | |
| 444 C You must have | |
| 4133 | 445 C LRW .GE. 50+(MAXORD+4)*NEQ+NEQ**2+3*NG |
| 4000 | 446 C for the full (dense) JACOBIAN case (when INFO(6)=0), or |
| 4133 | 447 C LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ+3*NG |
| 4000 | 448 C for the banded user-defined JACOBIAN case |
| 449 C (when INFO(5)=1 and INFO(6)=1), or | |
| 450 C LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ | |
| 4133 | 451 C +2*(NEQ/(ML+MU+1)+1)+3*NG |
| 4000 | 452 C for the banded finite-difference-generated JACOBIAN case |
| 453 C (when INFO(5)=0 and INFO(6)=1) | |
| 454 C | |
| 455 C IWORK(*) -- Dimension this integer work array of length LIW in | |
| 456 C your calling program. | |
| 457 C | |
| 458 C LIW -- Set it to the declared length of the IWORK array. | |
| 4428 | 459 C you must have LIW .GE. 21+NEQ |
| 4000 | 460 C |
| 461 C RPAR, IPAR -- These are parameter arrays, of real and integer | |
| 462 C type, respectively. You can use them for communication | |
| 463 C between your program that calls DDASRT and the | |
| 464 C RES subroutine (and the JAC subroutine). They are not | |
| 465 C altered by DDASRT. If you do not need RPAR or IPAR, | |
| 466 C ignore these parameters by treating them as dummy | |
| 467 C arguments. If you do choose to use them, dimension | |
| 468 C them in your calling program and in RES (and in JAC) | |
| 469 C as arrays of appropriate length. | |
| 470 C | |
| 471 C JAC -- If you have set INFO(5)=0, you can ignore this parameter | |
| 472 C by treating it as a dummy argument. Otherwise, you must | |
| 473 C provide a subroutine of the form | |
| 474 C JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR) | |
| 475 C to define the matrix of partial derivatives | |
| 476 C PD=DG/DY+CJ*DG/DYPRIME | |
| 477 C CJ is a scalar which is input to JAC. | |
| 478 C For the given values of T,Y,YPRIME, the | |
| 479 C subroutine must evaluate the non-zero partial | |
| 480 C derivatives for each equation and each solution | |
| 481 C component, and store these values in the | |
| 482 C matrix PD. The elements of PD are set to zero | |
| 483 C before each call to JAC so only non-zero elements | |
| 484 C need to be defined. | |
| 485 C | |
| 486 C Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ. | |
| 487 C You must declare the name JAC in an | |
| 488 C EXTERNAL STATEMENT in your program that calls | |
| 489 C DDASRT. You must dimension Y, YPRIME and PD | |
| 490 C in JAC. | |
| 491 C | |
| 492 C The way you must store the elements into the PD matrix | |
| 493 C depends on the structure of the matrix which you | |
| 494 C indicated by INFO(6). | |
| 495 C *** INFO(6)=0 -- Full (dense) matrix *** | |
| 496 C Give PD a first dimension of NEQ. | |
| 497 C When you evaluate the (non-zero) partial derivative | |
| 498 C of equation I with respect to variable J, you must | |
| 499 C store it in PD according to | |
| 500 C PD(I,J) = * DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)* | |
| 501 C *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU | |
| 502 C upper diagonal bands (refer to INFO(6) description | |
| 503 C of ML and MU) *** | |
| 504 C Give PD a first dimension of 2*ML+MU+1. | |
| 505 C when you evaluate the (non-zero) partial derivative | |
| 506 C of equation I with respect to variable J, you must | |
| 507 C store it in PD according to | |
| 508 C IROW = I - J + ML + MU + 1 | |
| 509 C PD(IROW,J) = *DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)* | |
| 510 C RPAR and IPAR are real and integer parameter arrays | |
| 511 C which you can use for communication between your calling | |
| 512 C program and your JACOBIAN subroutine JAC. They are not | |
| 513 C altered by DDASRT. If you do not need RPAR or IPAR, | |
| 514 C ignore these parameters by treating them as dummy | |
| 515 C arguments. If you do choose to use them, dimension | |
| 516 C them in your calling program and in JAC as arrays of | |
| 517 C appropriate length. | |
| 518 C | |
| 519 C G -- This is the name of the subroutine for defining constraint | |
| 520 C functions, whose roots are desired during the | |
| 521 C integration. It is to have the form | |
| 522 C SUBROUTINE G(NEQ,T,Y,NG,GOUT,RPAR,IPAR) | |
| 523 C DIMENSION Y(NEQ),GOUT(NG), | |
| 524 C where NEQ, T, Y and NG are INPUT, and the array GOUT is | |
| 525 C output. NEQ, T, and Y have the same meaning as in the | |
| 526 C RES routine, and GOUT is an array of length NG. | |
| 527 C For I=1,...,NG, this routine is to load into GOUT(I) | |
| 528 C the value at (T,Y) of the I-th constraint function G(I). | |
| 529 C DDASRT will find roots of the G(I) of odd multiplicity | |
| 530 C (that is, sign changes) as they occur during | |
| 531 C the integration. G must be declared EXTERNAL in the | |
| 532 C calling program. | |
| 533 C | |
| 534 C CAUTION..because of numerical errors in the functions | |
| 535 C G(I) due to roundoff and integration error, DDASRT | |
| 536 C may return false roots, or return the same root at two | |
| 537 C or more nearly equal values of T. If such false roots | |
| 538 C are suspected, the user should consider smaller error | |
| 539 C tolerances and/or higher precision in the evaluation of | |
| 540 C the G(I). | |
| 541 C | |
| 542 C If a root of some G(I) defines the end of the problem, | |
| 543 C the input to DDASRT should nevertheless allow | |
| 544 C integration to a point slightly past that ROOT, so | |
| 545 C that DDASRT can locate the root by interpolation. | |
| 546 C | |
| 547 C NG -- The number of constraint functions G(I). If there are none, | |
| 548 C set NG = 0, and pass a dummy name for G. | |
| 549 C | |
| 550 C JROOT -- This is an integer array of length NG. It is used only for | |
| 551 C output. On a return where one or more roots have been | |
| 552 C found, JROOT(I)=1 If G(I) has a root at T, | |
| 553 C or JROOT(I)=0 if not. | |
| 554 C | |
| 555 C | |
| 556 C | |
| 557 C OPTIONALLY REPLACEABLE NORM ROUTINE: | |
| 558 C DDASRT uses a weighted norm DDANRM to measure the size | |
| 559 C of vectors such as the estimated error in each step. | |
| 560 C A FUNCTION subprogram | |
| 561 C DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR) | |
| 562 C DIMENSION V(NEQ),WT(NEQ) | |
| 563 C is used to define this norm. Here, V is the vector | |
| 564 C whose norm is to be computed, and WT is a vector of | |
| 565 C weights. A DDANRM routine has been included with DDASRT | |
| 566 C which computes the weighted root-mean-square norm | |
| 567 C given by | |
| 568 C DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2) | |
| 569 C this norm is suitable for most problems. In some | |
| 570 C special cases, it may be more convenient and/or | |
| 571 C efficient to define your own norm by writing a function | |
| 572 C subprogram to be called instead of DDANRM. This should | |
| 573 C ,however, be attempted only after careful thought and | |
| 574 C consideration. | |
| 575 C | |
| 576 C | |
| 577 C------OUTPUT-AFTER ANY RETURN FROM DDASRT---- | |
| 578 C | |
| 579 C The principal aim of the code is to return a computed solution at | |
| 580 C TOUT, although it is also possible to obtain intermediate results | |
| 581 C along the way. To find out whether the code achieved its goal | |
| 582 C or if the integration process was interrupted before the task was | |
| 583 C completed, you must check the IDID parameter. | |
| 584 C | |
| 585 C | |
| 586 C T -- The solution was successfully advanced to the | |
| 587 C output value of T. | |
| 588 C | |
| 589 C Y(*) -- Contains the computed solution approximation at T. | |
| 590 C | |
| 591 C YPRIME(*) -- Contains the computed derivative | |
| 592 C approximation at T. | |
| 593 C | |
| 594 C IDID -- Reports what the code did. | |
| 595 C | |
| 596 C *** Task completed *** | |
| 597 C Reported by positive values of IDID | |
| 598 C | |
| 599 C IDID = 1 -- A step was successfully taken in the | |
| 600 C intermediate-output mode. The code has not | |
| 601 C yet reached TOUT. | |
| 602 C | |
| 603 C IDID = 2 -- The integration to TSTOP was successfully | |
| 604 C completed (T=TSTOP) by stepping exactly to TSTOP. | |
| 605 C | |
| 606 C IDID = 3 -- The integration to TOUT was successfully | |
| 607 C completed (T=TOUT) by stepping past TOUT. | |
| 608 C Y(*) is obtained by interpolation. | |
| 609 C YPRIME(*) is obtained by interpolation. | |
| 610 C | |
| 611 C IDID = 4 -- The integration was successfully completed | |
| 612 C by finding one or more roots of G at T. | |
| 613 C | |
| 614 C *** Task interrupted *** | |
| 615 C Reported by negative values of IDID | |
| 616 C | |
| 617 C IDID = -1 -- A large amount of work has been expended. | |
| 618 C (About INFO(12) steps) | |
| 619 C | |
| 620 C IDID = -2 -- The error tolerances are too stringent. | |
| 621 C | |
| 622 C IDID = -3 -- The local error test cannot be satisfied | |
| 623 C because you specified a zero component in ATOL | |
| 624 C and the corresponding computed solution | |
| 625 C component is zero. Thus, a pure relative error | |
| 626 C test is impossible for this component. | |
| 627 C | |
| 628 C IDID = -6 -- DDASRT had repeated error test | |
| 629 C failures on the last attempted step. | |
| 630 C | |
| 631 C IDID = -7 -- The corrector could not converge. | |
| 632 C | |
| 633 C IDID = -8 -- The matrix of partial derivatives | |
| 634 C is singular. | |
| 635 C | |
| 636 C IDID = -9 -- The corrector could not converge. | |
| 637 C there were repeated error test failures | |
| 638 C in this step. | |
| 639 C | |
| 640 C IDID =-10 -- The corrector could not converge | |
| 641 C because IRES was equal to minus one. | |
| 642 C | |
| 643 C IDID =-11 -- IRES equal to -2 was encountered | |
| 644 C and control is being returned to the | |
| 645 C calling program. | |
| 646 C | |
| 647 C IDID =-12 -- DDASRT failed to compute the initial | |
| 648 C YPRIME. | |
| 649 C | |
| 650 C | |
| 651 C | |
| 652 C IDID = -13,..,-32 -- Not applicable for this code | |
| 653 C | |
| 654 C *** Task terminated *** | |
| 655 C Reported by the value of IDID=-33 | |
| 656 C | |
| 657 C IDID = -33 -- The code has encountered trouble from which | |
| 658 C it cannot recover. A message is printed | |
| 659 C explaining the trouble and control is returned | |
| 660 C to the calling program. For example, this occurs | |
| 661 C when invalid input is detected. | |
| 662 C | |
| 663 C RTOL, ATOL -- These quantities remain unchanged except when | |
| 664 C IDID = -2. In this case, the error tolerances have been | |
| 665 C increased by the code to values which are estimated to | |
| 666 C be appropriate for continuing the integration. However, | |
| 667 C the reported solution at T was obtained using the input | |
| 668 C values of RTOL and ATOL. | |
| 669 C | |
| 670 C RWORK, IWORK -- Contain information which is usually of no | |
| 671 C interest to the user but necessary for subsequent calls. | |
| 672 C However, you may find use for | |
| 673 C | |
| 674 C RWORK(3)--Which contains the step size H to be | |
| 675 C attempted on the next step. | |
| 676 C | |
| 677 C RWORK(4)--Which contains the current value of the | |
| 678 C independent variable, i.e., the farthest point | |
| 679 C integration has reached. This will be different | |
| 680 C from T only when interpolation has been | |
| 681 C performed (IDID=3). | |
| 682 C | |
| 683 C RWORK(7)--Which contains the stepsize used | |
| 684 C on the last successful step. | |
| 685 C | |
| 686 C IWORK(7)--Which contains the order of the method to | |
| 687 C be attempted on the next step. | |
| 688 C | |
| 689 C IWORK(8)--Which contains the order of the method used | |
| 690 C on the last step. | |
| 691 C | |
| 692 C IWORK(11)--Which contains the number of steps taken so | |
| 693 C far. | |
| 694 C | |
| 695 C IWORK(12)--Which contains the number of calls to RES | |
| 696 C so far. | |
| 697 C | |
| 698 C IWORK(13)--Which contains the number of evaluations of | |
| 699 C the matrix of partial derivatives needed so | |
| 700 C far. | |
| 701 C | |
| 702 C IWORK(14)--Which contains the total number | |
| 703 C of error test failures so far. | |
| 704 C | |
| 705 C IWORK(15)--Which contains the total number | |
| 706 C of convergence test failures so far. | |
| 707 C (includes singular iteration matrix | |
| 708 C failures.) | |
| 709 C | |
| 710 C IWORK(16)--Which contains the total number of calls | |
| 711 C to the constraint function g so far | |
| 712 C | |
| 713 C | |
| 714 C | |
| 715 C INPUT -- What to do to continue the integration | |
| 716 C (calls after the first) ** | |
| 717 C | |
| 718 C This code is organized so that subsequent calls to continue the | |
| 719 C integration involve little (if any) additional effort on your | |
| 720 C part. You must monitor the IDID parameter in order to determine | |
| 721 C what to do next. | |
| 722 C | |
| 723 C Recalling that the principal task of the code is to integrate | |
| 724 C from T to TOUT (the interval mode), usually all you will need | |
| 725 C to do is specify a new TOUT upon reaching the current TOUT. | |
| 726 C | |
| 727 C Do not alter any quantity not specifically permitted below, | |
| 728 C in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*) | |
| 729 C or the differential equation in subroutine RES. Any such | |
| 730 C alteration constitutes a new problem and must be treated as such, | |
| 731 C i.e., you must start afresh. | |
| 732 C | |
| 733 C You cannot change from vector to scalar error control or vice | |
| 734 C versa (INFO(2)), but you can change the size of the entries of | |
| 735 C RTOL, ATOL. Increasing a tolerance makes the equation easier | |
| 736 C to integrate. Decreasing a tolerance will make the equation | |
| 737 C harder to integrate and should generally be avoided. | |
| 738 C | |
| 739 C You can switch from the intermediate-output mode to the | |
| 740 C interval mode (INFO(3)) or vice versa at any time. | |
| 741 C | |
| 742 C If it has been necessary to prevent the integration from going | |
| 743 C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the | |
| 744 C code will not integrate to any TOUT beyond the currently | |
| 745 C specified TSTOP. Once TSTOP has been reached you must change | |
| 746 C the value of TSTOP or set INFO(4)=0. You may change INFO(4) | |
| 747 C or TSTOP at any time but you must supply the value of TSTOP in | |
| 748 C RWORK(1) whenever you set INFO(4)=1. | |
| 749 C | |
| 750 C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2) | |
| 751 C unless you are going to restart the code. | |
| 752 C | |
| 753 C *** Following a completed task *** | |
| 754 C If | |
| 755 C IDID = 1, call the code again to continue the integration | |
| 756 C another step in the direction of TOUT. | |
| 757 C | |
| 758 C IDID = 2 or 3, define a new TOUT and call the code again. | |
| 759 C TOUT must be different from T. You cannot change | |
| 760 C the direction of integration without restarting. | |
| 761 C | |
| 762 C IDID = 4, call the code again to continue the integration | |
| 763 C another step in the direction of TOUT. You may | |
| 764 C change the functions in G after a return with IDID=4, | |
| 765 C but the number of constraint functions NG must remain | |
| 766 C the same. If you wish to change | |
| 767 C the functions in RES or in G, then you | |
| 768 C must restart the code. | |
| 769 C | |
| 770 C *** Following an interrupted task *** | |
| 771 C To show the code that you realize the task was | |
| 772 C interrupted and that you want to continue, you | |
| 773 C must take appropriate action and set INFO(1) = 1 | |
| 774 C If | |
| 775 C IDID = -1, The code has reached the step iteration. | |
| 776 C If you want to continue, set INFO(1) = 1 and | |
| 777 C call the code again. See also INFO(12). | |
| 778 C | |
| 779 C IDID = -2, The error tolerances RTOL, ATOL have been | |
| 780 C increased to values the code estimates appropriate | |
| 781 C for continuing. You may want to change them | |
| 782 C yourself. If you are sure you want to continue | |
| 783 C with relaxed error tolerances, set INFO(1)=1 and | |
| 784 C call the code again. | |
| 785 C | |
| 786 C IDID = -3, A solution component is zero and you set the | |
| 787 C corresponding component of ATOL to zero. If you | |
| 788 C are sure you want to continue, you must first | |
| 789 C alter the error criterion to use positive values | |
| 790 C for those components of ATOL corresponding to zero | |
| 791 C solution components, then set INFO(1)=1 and call | |
| 792 C the code again. | |
| 793 C | |
| 794 C IDID = -4,-5 --- Cannot occur with this code. | |
| 795 C | |
| 796 C IDID = -6, Repeated error test failures occurred on the | |
| 797 C last attempted step in DDASRT. A singularity in the | |
| 798 C solution may be present. If you are absolutely | |
| 799 C certain you want to continue, you should restart | |
| 800 C the integration. (Provide initial values of Y and | |
| 801 C YPRIME which are consistent) | |
| 802 C | |
| 803 C IDID = -7, Repeated convergence test failures occurred | |
| 804 C on the last attempted step in DDASRT. An inaccurate | |
| 805 C or ill-conditioned JACOBIAN may be the problem. If | |
| 806 C you are absolutely certain you want to continue, you | |
| 807 C should restart the integration. | |
| 808 C | |
| 809 C IDID = -8, The matrix of partial derivatives is singular. | |
| 810 C Some of your equations may be redundant. | |
| 811 C DDASRT cannot solve the problem as stated. | |
| 812 C It is possible that the redundant equations | |
| 813 C could be removed, and then DDASRT could | |
| 814 C solve the problem. It is also possible | |
| 815 C that a solution to your problem either | |
| 816 C does not exist or is not unique. | |
| 817 C | |
| 818 C IDID = -9, DDASRT had multiple convergence test | |
| 819 C failures, preceeded by multiple error | |
| 820 C test failures, on the last attempted step. | |
| 821 C It is possible that your problem | |
| 822 C is ill-posed, and cannot be solved | |
| 823 C using this code. Or, there may be a | |
| 824 C discontinuity or a singularity in the | |
| 825 C solution. If you are absolutely certain | |
| 826 C you want to continue, you should restart | |
| 827 C the integration. | |
| 828 C | |
| 829 C IDID =-10, DDASRT had multiple convergence test failures | |
| 830 C because IRES was equal to minus one. | |
| 831 C If you are absolutely certain you want | |
| 832 C to continue, you should restart the | |
| 833 C integration. | |
| 834 C | |
| 835 C IDID =-11, IRES=-2 was encountered, and control is being | |
| 836 C returned to the calling program. | |
| 837 C | |
| 838 C IDID =-12, DDASRT failed to compute the initial YPRIME. | |
| 839 C This could happen because the initial | |
| 840 C approximation to YPRIME was not very good, or | |
| 841 C if a YPRIME consistent with the initial Y | |
| 842 C does not exist. The problem could also be caused | |
| 843 C by an inaccurate or singular iteration matrix. | |
| 844 C | |
| 845 C | |
| 846 C | |
| 847 C IDID = -13,..,-32 --- Cannot occur with this code. | |
| 848 C | |
| 849 C *** Following a terminated task *** | |
| 850 C If IDID= -33, you cannot continue the solution of this | |
| 851 C problem. An attempt to do so will result in your | |
| 852 C run being terminated. | |
| 853 C | |
| 854 C --------------------------------------------------------------------- | |
| 855 C | |
| 856 C***REFERENCE | |
|
19627
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
15271
diff
changeset
|
857 C K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical |
| 4000 | 858 C Solution of Initial-Value Problems in Differential-Algebraic |
| 859 C Equations, Elsevier, New York, 1989. | |
| 860 C | |
| 861 C***ROUTINES CALLED DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,DRCHEK,DROOTS, | |
| 4040 | 862 C XERRWD,D1MACH |
| 4000 | 863 C***END PROLOGUE DDASRT |
| 864 C | |
| 865 C**End | |
| 866 C | |
| 867 IMPLICIT DOUBLE PRECISION(A-H,O-Z) | |
| 868 LOGICAL DONE | |
| 869 EXTERNAL RES, JAC, G | |
| 870 DIMENSION Y(*),YPRIME(*) | |
| 871 DIMENSION INFO(15) | |
| 872 DIMENSION RWORK(*),IWORK(*) | |
| 873 DIMENSION RTOL(*),ATOL(*) | |
| 874 DIMENSION RPAR(*),IPAR(*) | |
| 875 CHARACTER MSG*80 | |
| 876 C | |
| 877 C SET POINTERS INTO IWORK | |
| 878 PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11, | |
| 879 * LNRE=12, LNJE=13, LETF=14, LCTF=15, LNGE=16, LNPD=17, | |
| 4428 | 880 * LIRFND=18, LMXSTP=21, LIPVT=22, LJCALC=5, LPHASE=6, LK=7, |
| 4000 | 881 * LKOLD=8, LNS=9, LNSTL=10, LIWM=1) |
| 882 C | |
| 883 C SET RELATIVE OFFSET INTO RWORK | |
| 884 PARAMETER (NPD=1) | |
| 885 C | |
| 886 C SET POINTERS INTO RWORK | |
| 887 PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, | |
| 888 * LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9, | |
| 889 * LALPHA=11, LBETA=17, LGAMMA=23, | |
| 890 * LPSI=29, LSIGMA=35, LT0=41, LTLAST=42, LALPHR=43, LX2=44, | |
| 891 * LDELTA=51) | |
| 892 C | |
| 893 C***FIRST EXECUTABLE STATEMENT DDASRT | |
| 894 IF(INFO(1).NE.0)GO TO 100 | |
| 895 C | |
| 896 C----------------------------------------------------------------------- | |
| 897 C THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY. | |
| 898 C IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS. | |
| 899 C----------------------------------------------------------------------- | |
| 900 C | |
| 901 C FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO | |
| 902 C ARE EITHER ZERO OR ONE. | |
| 903 DO 10 I=2,12 | |
| 904 IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701 | |
| 905 10 CONTINUE | |
| 906 C | |
| 907 IF(NEQ.LE.0)GO TO 702 | |
| 908 C | |
| 909 C CHECK AND COMPUTE MAXIMUM ORDER | |
| 910 MXORD=5 | |
| 911 IF(INFO(9).EQ.0)GO TO 20 | |
| 912 MXORD=IWORK(LMXORD) | |
| 913 IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703 | |
| 914 20 IWORK(LMXORD)=MXORD | |
| 915 C | |
| 916 C COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU. | |
| 917 IF(INFO(6).NE.0)GO TO 40 | |
| 918 LENPD=NEQ**2 | |
| 4133 | 919 LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+3*NG |
| 4000 | 920 IF(INFO(5).NE.0)GO TO 30 |
| 921 IWORK(LMTYPE)=2 | |
| 922 GO TO 60 | |
| 923 30 IWORK(LMTYPE)=1 | |
| 924 GO TO 60 | |
| 925 40 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 | |
| 926 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 | |
| 927 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ | |
| 928 IF(INFO(5).NE.0)GO TO 50 | |
| 929 IWORK(LMTYPE)=5 | |
| 930 MBAND=IWORK(LML)+IWORK(LMU)+1 | |
| 931 MSAVE=(NEQ/MBAND)+1 | |
| 4133 | 932 LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE+3*NG |
| 4000 | 933 GO TO 60 |
| 934 50 IWORK(LMTYPE)=4 | |
| 4133 | 935 LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+3*NG |
| 4000 | 936 C |
| 937 C CHECK LENGTHS OF RWORK AND IWORK | |
| 4428 | 938 60 LENIW=21+NEQ |
| 4000 | 939 IWORK(LNPD)=LENPD |
| 940 IF(LRW.LT.LENRW)GO TO 704 | |
| 941 IF(LIW.LT.LENIW)GO TO 705 | |
| 942 C | |
| 943 C CHECK TO SEE THAT TOUT IS DIFFERENT FROM T | |
| 944 C Also check to see that NG is larger than 0. | |
| 945 IF(TOUT .EQ. T)GO TO 719 | |
| 946 IF(NG .LT. 0) GO TO 730 | |
| 947 C | |
| 948 C CHECK HMAX | |
| 949 IF(INFO(7).EQ.0)GO TO 70 | |
| 950 HMAX=RWORK(LHMAX) | |
| 951 IF(HMAX.LE.0.0D0)GO TO 710 | |
| 952 70 CONTINUE | |
| 953 C | |
| 954 C CHECK AND COMPUTE MAXIMUM STEPS | |
| 955 MXSTP=500 | |
| 956 IF(INFO(12).EQ.0)GO TO 80 | |
| 957 MXSTP=IWORK(LMXSTP) | |
| 4429 | 958 IF(MXSTP.LT.0)GO TO 716 |
| 4000 | 959 80 IWORK(LMXSTP)=MXSTP |
| 960 C | |
| 961 C INITIALIZE COUNTERS | |
| 962 IWORK(LNST)=0 | |
| 963 IWORK(LNRE)=0 | |
| 964 IWORK(LNJE)=0 | |
| 965 IWORK(LNGE)=0 | |
| 966 C | |
| 967 IWORK(LNSTL)=0 | |
| 968 IDID=1 | |
| 969 GO TO 200 | |
| 970 C | |
| 971 C----------------------------------------------------------------------- | |
| 972 C THIS BLOCK IS FOR CONTINUATION CALLS | |
| 973 C ONLY. HERE WE CHECK INFO(1),AND IF THE | |
| 974 C LAST STEP WAS INTERRUPTED WE CHECK WHETHER | |
| 975 C APPROPRIATE ACTION WAS TAKEN. | |
| 976 C----------------------------------------------------------------------- | |
| 977 C | |
| 978 100 CONTINUE | |
| 979 IF(INFO(1).EQ.1)GO TO 110 | |
| 980 IF(INFO(1).NE.-1)GO TO 701 | |
| 981 C IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED | |
| 982 C BY AN ERROR CONDITION FROM DDASTP,AND | |
| 983 C APPROPRIATE ACTION WAS NOT TAKEN. THIS | |
| 984 C IS A FATAL ERROR. | |
| 985 MSG = 'DASRT-- THE LAST STEP TERMINATED WITH A NEGATIVE' | |
| 4040 | 986 CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 987 MSG = 'DASRT-- VALUE (=I1) OF IDID AND NO APPROPRIATE' |
| 4040 | 988 CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0) |
| 4000 | 989 MSG = 'DASRT-- ACTION WAS TAKEN. RUN TERMINATED' |
| 4040 | 990 CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 991 RETURN |
| 992 110 CONTINUE | |
| 993 IWORK(LNSTL)=IWORK(LNST) | |
| 994 C | |
| 995 C----------------------------------------------------------------------- | |
| 996 C THIS BLOCK IS EXECUTED ON ALL CALLS. | |
| 997 C THE ERROR TOLERANCE PARAMETERS ARE | |
| 998 C CHECKED, AND THE WORK ARRAY POINTERS | |
| 999 C ARE SET. | |
| 1000 C----------------------------------------------------------------------- | |
| 1001 C | |
| 1002 200 CONTINUE | |
| 1003 C CHECK RTOL,ATOL | |
| 1004 NZFLG=0 | |
| 1005 RTOLI=RTOL(1) | |
| 1006 ATOLI=ATOL(1) | |
| 1007 DO 210 I=1,NEQ | |
| 1008 IF(INFO(2).EQ.1)RTOLI=RTOL(I) | |
| 1009 IF(INFO(2).EQ.1)ATOLI=ATOL(I) | |
| 1010 IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1 | |
| 1011 IF(RTOLI.LT.0.0D0)GO TO 706 | |
| 1012 IF(ATOLI.LT.0.0D0)GO TO 707 | |
| 1013 210 CONTINUE | |
| 1014 IF(NZFLG.EQ.0)GO TO 708 | |
| 1015 C | |
| 1016 C SET UP RWORK STORAGE.IWORK STORAGE IS FIXED | |
| 1017 C IN DATA STATEMENT. | |
| 1018 LG0=LDELTA+NEQ | |
| 1019 LG1=LG0+NG | |
| 1020 LGX=LG1+NG | |
| 1021 LE=LGX+NG | |
| 1022 LWT=LE+NEQ | |
| 1023 LPHI=LWT+NEQ | |
| 1024 LPD=LPHI+(IWORK(LMXORD)+1)*NEQ | |
| 1025 LWM=LPD | |
| 1026 NTEMP=NPD+IWORK(LNPD) | |
| 1027 IF(INFO(1).EQ.1)GO TO 400 | |
| 1028 C | |
| 1029 C----------------------------------------------------------------------- | |
| 1030 C THIS BLOCK IS EXECUTED ON THE INITIAL CALL | |
| 1031 C ONLY. SET THE INITIAL STEP SIZE, AND | |
| 1032 C THE ERROR WEIGHT VECTOR, AND PHI. | |
| 1033 C COMPUTE INITIAL YPRIME, IF NECESSARY. | |
| 1034 C----------------------------------------------------------------------- | |
| 1035 C | |
| 1036 300 CONTINUE | |
| 1037 TN=T | |
| 1038 IDID=1 | |
| 1039 C | |
| 1040 C SET ERROR WEIGHT VECTOR WT | |
| 1041 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) | |
| 1042 DO 305 I = 1,NEQ | |
| 1043 IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713 | |
| 1044 305 CONTINUE | |
| 1045 C | |
| 1046 C COMPUTE UNIT ROUNDOFF AND HMIN | |
| 1047 UROUND = D1MACH(4) | |
| 1048 RWORK(LROUND) = UROUND | |
| 1049 HMIN = 4.0D0*UROUND*DMAX1(DABS(T),DABS(TOUT)) | |
| 1050 C | |
| 1051 C CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH | |
| 1052 TDIST = DABS(TOUT - T) | |
| 1053 IF(TDIST .LT. HMIN) GO TO 714 | |
| 1054 C | |
| 1055 C CHECK H0, IF THIS WAS INPUT | |
| 1056 IF (INFO(8) .EQ. 0) GO TO 310 | |
| 1057 HO = RWORK(LH) | |
| 1058 IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711 | |
| 1059 IF (HO .EQ. 0.0D0) GO TO 712 | |
| 1060 GO TO 320 | |
| 1061 310 CONTINUE | |
| 1062 C | |
| 1063 C COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER | |
| 1064 C DDASTP OR DDAINI, DEPENDING ON INFO(11) | |
| 1065 HO = 0.001D0*TDIST | |
| 1066 YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR) | |
| 1067 IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM | |
| 1068 HO = DSIGN(HO,TOUT-T) | |
| 1069 C ADJUST HO IF NECESSARY TO MEET HMAX BOUND | |
| 1070 320 IF (INFO(7) .EQ. 0) GO TO 330 | |
| 1071 RH = DABS(HO)/RWORK(LHMAX) | |
| 1072 IF (RH .GT. 1.0D0) HO = HO/RH | |
| 1073 C COMPUTE TSTOP, IF APPLICABLE | |
| 1074 330 IF (INFO(4) .EQ. 0) GO TO 340 | |
| 1075 TSTOP = RWORK(LTSTOP) | |
| 1076 IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715 | |
| 1077 IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T | |
| 1078 IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709 | |
| 1079 C | |
| 1080 C COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE | |
| 1081 340 IF (INFO(11) .EQ. 0) GO TO 350 | |
| 1082 CALL DDAINI(TN,Y,YPRIME,NEQ, | |
| 1083 * RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR, | |
| 1084 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), | |
| 1085 * RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND), | |
| 1086 * INFO(10),NTEMP) | |
| 1087 IF (IDID .LT. 0) GO TO 390 | |
| 1088 C | |
| 1089 C LOAD H WITH H0. STORE H IN RWORK(LH) | |
| 1090 350 H = HO | |
| 1091 RWORK(LH) = H | |
| 1092 C | |
| 1093 C LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2) | |
| 1094 360 ITEMP = LPHI + NEQ | |
| 1095 DO 370 I = 1,NEQ | |
| 1096 RWORK(LPHI + I - 1) = Y(I) | |
| 1097 370 RWORK(ITEMP + I - 1) = H*YPRIME(I) | |
| 1098 C | |
| 1099 C INITIALIZE T0 IN RWORK AND CHECK FOR A ZERO OF G NEAR THE | |
| 1100 C INITIAL T. | |
| 1101 C | |
| 1102 RWORK(LT0) = T | |
| 1103 IWORK(LIRFND) = 0 | |
| 1104 RWORK(LPSI)=H | |
| 1105 RWORK(LPSI+1)=2.0D0*H | |
| 1106 IWORK(LKOLD)=1 | |
| 1107 IF(NG .EQ. 0) GO TO 390 | |
| 1108 CALL DRCHEK(1,G,NG,NEQ,T,TOUT,Y,RWORK(LE),RWORK(LPHI), | |
| 1109 * RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1), | |
| 1110 * RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3), | |
| 1111 * RWORK,IWORK,RPAR,IPAR) | |
| 1112 IF(IRT .NE. 0) GO TO 732 | |
| 1113 C | |
| 1114 C Check for a root in the interval (T0,TN], unless DDASRT | |
| 1115 C did not have to initialize YPRIME. | |
| 1116 C | |
| 1117 IF(NG .EQ. 0 .OR. INFO(11) .EQ. 0) GO TO 390 | |
| 1118 CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI), | |
| 1119 * RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1), | |
| 1120 * RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3), | |
| 1121 * RWORK,IWORK,RPAR,IPAR) | |
| 1122 IF(IRT .NE. 1) GO TO 390 | |
| 1123 IWORK(LIRFND) = 1 | |
| 1124 IDID = 4 | |
| 1125 T = RWORK(LT0) | |
| 1126 GO TO 580 | |
| 1127 C | |
| 1128 390 GO TO 500 | |
| 1129 C | |
| 1130 C------------------------------------------------------- | |
| 1131 C THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS | |
| 1132 C PURPOSE IS TO CHECK STOP CONDITIONS BEFORE | |
| 1133 C TAKING A STEP. | |
| 1134 C ADJUST H IF NECESSARY TO MEET HMAX BOUND | |
| 1135 C------------------------------------------------------- | |
| 1136 C | |
| 1137 400 CONTINUE | |
| 1138 UROUND=RWORK(LROUND) | |
| 1139 DONE = .FALSE. | |
| 1140 TN=RWORK(LTN) | |
| 1141 H=RWORK(LH) | |
| 1142 IF(NG .EQ. 0) GO TO 405 | |
| 1143 C | |
| 1144 C Check for a zero of G near TN. | |
| 1145 C | |
| 1146 CALL DRCHEK(2,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI), | |
| 1147 * RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1), | |
| 1148 * RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3), | |
| 1149 * RWORK,IWORK,RPAR,IPAR) | |
| 1150 IF(IRT .NE. 1) GO TO 405 | |
| 1151 IWORK(LIRFND) = 1 | |
| 1152 IDID = 4 | |
| 1153 T = RWORK(LT0) | |
| 1154 DONE = .TRUE. | |
| 1155 GO TO 490 | |
| 1156 C | |
| 1157 405 CONTINUE | |
| 1158 IF(INFO(7) .EQ. 0) GO TO 410 | |
| 1159 RH = DABS(H)/RWORK(LHMAX) | |
| 1160 IF(RH .GT. 1.0D0) H = H/RH | |
| 1161 410 CONTINUE | |
| 1162 IF(T .EQ. TOUT) GO TO 719 | |
| 1163 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 | |
| 1164 IF(INFO(4) .EQ. 1) GO TO 430 | |
| 1165 IF(INFO(3) .EQ. 1) GO TO 420 | |
| 1166 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 | |
| 1167 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1168 * RWORK(LPHI),RWORK(LPSI)) | |
| 1169 T=TOUT | |
| 1170 IDID = 3 | |
| 1171 DONE = .TRUE. | |
| 1172 GO TO 490 | |
| 1173 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 | |
| 1174 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 | |
| 1175 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1176 * RWORK(LPHI),RWORK(LPSI)) | |
| 1177 T = TN | |
| 1178 IDID = 1 | |
| 1179 DONE = .TRUE. | |
| 1180 GO TO 490 | |
| 1181 425 CONTINUE | |
| 1182 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1183 * RWORK(LPHI),RWORK(LPSI)) | |
| 1184 T = TOUT | |
| 1185 IDID = 3 | |
| 1186 DONE = .TRUE. | |
| 1187 GO TO 490 | |
| 1188 430 IF(INFO(3) .EQ. 1) GO TO 440 | |
| 1189 TSTOP=RWORK(LTSTOP) | |
| 1190 IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 | |
| 1191 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 | |
| 1192 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 | |
| 1193 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1194 * RWORK(LPHI),RWORK(LPSI)) | |
| 1195 T=TOUT | |
| 1196 IDID = 3 | |
| 1197 DONE = .TRUE. | |
| 1198 GO TO 490 | |
| 1199 440 TSTOP = RWORK(LTSTOP) | |
| 1200 IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 | |
| 1201 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 | |
| 1202 IF((TN-T)*H .LE. 0.0D0) GO TO 450 | |
| 1203 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 | |
| 1204 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1205 * RWORK(LPHI),RWORK(LPSI)) | |
| 1206 T = TN | |
| 1207 IDID = 1 | |
| 1208 DONE = .TRUE. | |
| 1209 GO TO 490 | |
| 1210 445 CONTINUE | |
| 1211 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1212 * RWORK(LPHI),RWORK(LPSI)) | |
| 1213 T = TOUT | |
| 1214 IDID = 3 | |
| 1215 DONE = .TRUE. | |
| 1216 GO TO 490 | |
| 1217 450 CONTINUE | |
| 1218 C CHECK WHETHER WE ARE WITH IN ROUNDOFF OF TSTOP | |
| 1219 IF(DABS(TN-TSTOP).GT.100.0D0*UROUND* | |
| 1220 * (DABS(TN)+DABS(H)))GO TO 460 | |
| 1221 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD), | |
| 1222 * RWORK(LPHI),RWORK(LPSI)) | |
| 1223 IDID=2 | |
| 1224 T=TSTOP | |
| 1225 DONE = .TRUE. | |
| 1226 GO TO 490 | |
| 1227 460 TNEXT=TN+H | |
| 1228 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 | |
| 1229 H=TSTOP-TN | |
| 1230 RWORK(LH)=H | |
| 1231 C | |
| 1232 490 IF (DONE) GO TO 590 | |
| 1233 C | |
| 1234 C------------------------------------------------------- | |
| 1235 C THE NEXT BLOCK CONTAINS THE CALL TO THE | |
| 1236 C ONE-STEP INTEGRATOR DDASTP. | |
| 1237 C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS. | |
| 1238 C CHECK FOR TOO MANY STEPS. | |
| 1239 C UPDATE WT. | |
| 1240 C CHECK FOR TOO MUCH ACCURACY REQUESTED. | |
| 1241 C COMPUTE MINIMUM STEPSIZE. | |
| 1242 C------------------------------------------------------- | |
| 1243 C | |
| 1244 500 CONTINUE | |
| 1245 C CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME | |
| 1246 IF (IDID .EQ. -12) GO TO 527 | |
| 1247 C | |
| 1248 C CHECK FOR TOO MANY STEPS | |
| 1249 IF((IWORK(LNST)-IWORK(LNSTL)).LT.IWORK(LMXSTP)) | |
| 1250 * GO TO 510 | |
| 1251 IDID=-1 | |
| 1252 GO TO 527 | |
| 1253 C | |
| 1254 C UPDATE WT | |
| 1255 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI), | |
| 1256 * RWORK(LWT),RPAR,IPAR) | |
| 1257 DO 520 I=1,NEQ | |
| 1258 IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520 | |
| 1259 IDID=-3 | |
| 1260 GO TO 527 | |
| 1261 520 CONTINUE | |
| 1262 C | |
| 1263 C TEST FOR TOO MUCH ACCURACY REQUESTED. | |
| 1264 R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)* | |
| 1265 * 100.0D0*UROUND | |
| 1266 IF(R.LE.1.0D0)GO TO 525 | |
| 1267 C MULTIPLY RTOL AND ATOL BY R AND RETURN | |
| 1268 IF(INFO(2).EQ.1)GO TO 523 | |
| 1269 RTOL(1)=R*RTOL(1) | |
| 1270 ATOL(1)=R*ATOL(1) | |
| 1271 IDID=-2 | |
| 1272 GO TO 527 | |
| 1273 523 DO 524 I=1,NEQ | |
| 1274 RTOL(I)=R*RTOL(I) | |
| 1275 524 ATOL(I)=R*ATOL(I) | |
| 1276 IDID=-2 | |
| 1277 GO TO 527 | |
| 1278 525 CONTINUE | |
| 1279 C | |
| 1280 C COMPUTE MINIMUM STEPSIZE | |
| 1281 HMIN=4.0D0*UROUND*DMAX1(DABS(TN),DABS(TOUT)) | |
| 1282 C | |
| 1283 C TEST H VS. HMAX | |
| 1284 IF (INFO(7) .EQ. 0) GO TO 526 | |
| 1285 RH = ABS(H)/RWORK(LHMAX) | |
| 1286 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
|
1287 526 CONTINUE |
| 4000 | 1288 C |
| 1289 CALL DDASTP(TN,Y,YPRIME,NEQ, | |
| 1290 * RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR, | |
| 1291 * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), | |
| 1292 * RWORK(LWM),IWORK(LIWM), | |
| 1293 * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), | |
| 1294 * RWORK(LPSI),RWORK(LSIGMA), | |
| 1295 * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD), | |
| 1296 * RWORK(LS),HMIN,RWORK(LROUND), | |
| 1297 * IWORK(LPHASE),IWORK(LJCALC),IWORK(LK), | |
| 1298 * IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP) | |
| 1299 527 IF(IDID.LT.0)GO TO 600 | |
| 1300 C | |
| 1301 C-------------------------------------------------------- | |
| 1302 C THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN | |
| 1303 C FROM DDASTP (IDID=1). TEST FOR STOP CONDITIONS. | |
| 1304 C-------------------------------------------------------- | |
| 1305 C | |
| 1306 IF(NG .EQ. 0) GO TO 529 | |
| 1307 C | |
| 1308 C Check for a zero of G near TN. | |
| 1309 C | |
| 1310 CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI), | |
| 1311 * RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1), | |
| 1312 * RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3), | |
| 1313 * RWORK,IWORK,RPAR,IPAR) | |
| 1314 IF(IRT .NE. 1) GO TO 529 | |
| 1315 IWORK(LIRFND) = 1 | |
| 1316 IDID = 4 | |
| 1317 T = RWORK(LT0) | |
| 1318 GO TO 580 | |
| 1319 C | |
| 1320 529 CONTINUE | |
| 1321 IF(INFO(4).NE.0)GO TO 540 | |
| 1322 IF(INFO(3).NE.0)GO TO 530 | |
| 1323 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 | |
| 1324 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1325 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1326 IDID=3 | |
| 1327 T=TOUT | |
| 1328 GO TO 580 | |
| 1329 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 | |
| 1330 T=TN | |
| 1331 IDID=1 | |
| 1332 GO TO 580 | |
| 1333 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1334 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1335 IDID=3 | |
| 1336 T=TOUT | |
| 1337 GO TO 580 | |
| 1338 540 IF(INFO(3).NE.0)GO TO 550 | |
| 1339 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 | |
| 1340 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1341 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1342 T=TOUT | |
| 1343 IDID=3 | |
| 1344 GO TO 580 | |
| 1345 542 IF(DABS(TN-TSTOP).LE.100.0D0*UROUND* | |
| 1346 * (DABS(TN)+DABS(H)))GO TO 545 | |
| 1347 TNEXT=TN+H | |
| 1348 IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 | |
| 1349 H=TSTOP-TN | |
| 1350 GO TO 500 | |
| 1351 545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, | |
| 1352 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1353 IDID=2 | |
| 1354 T=TSTOP | |
| 1355 GO TO 580 | |
| 1356 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 | |
| 1357 IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*(DABS(TN)+DABS(H)))GO TO 552 | |
| 1358 T=TN | |
| 1359 IDID=1 | |
| 1360 GO TO 580 | |
| 1361 552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ, | |
| 1362 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1363 IDID=2 | |
| 1364 T=TSTOP | |
| 1365 GO TO 580 | |
| 1366 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, | |
| 1367 * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) | |
| 1368 T=TOUT | |
| 1369 IDID=3 | |
| 1370 580 CONTINUE | |
| 1371 C | |
| 1372 C-------------------------------------------------------- | |
| 1373 C ALL SUCCESSFUL RETURNS FROM DDASRT ARE MADE FROM | |
| 1374 C THIS BLOCK. | |
| 1375 C-------------------------------------------------------- | |
| 1376 C | |
| 1377 590 CONTINUE | |
| 1378 RWORK(LTN)=TN | |
| 1379 RWORK(LH)=H | |
| 1380 RWORK(LTLAST) = T | |
| 1381 RETURN | |
| 1382 C | |
| 1383 C----------------------------------------------------------------------- | |
| 1384 C THIS BLOCK HANDLES ALL UNSUCCESSFUL | |
| 1385 C RETURNS OTHER THAN FOR ILLEGAL INPUT. | |
| 1386 C----------------------------------------------------------------------- | |
| 1387 C | |
| 1388 600 CONTINUE | |
| 1389 ITEMP=-IDID | |
| 1390 GO TO (610,620,630,690,690,640,650,660,670,675, | |
| 1391 * 680,685), ITEMP | |
| 1392 C | |
| 1393 C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE | |
| 1394 C REACHING TOUT | |
| 1395 610 MSG = 'DASRT-- AT CURRENT T (=R1) 500 STEPS' | |
| 4040 | 1396 CALL XERRWD(MSG,38,610,0,0,0,0,1,TN,0.0D0) |
| 4000 | 1397 MSG = 'DASRT-- TAKEN ON THIS CALL BEFORE REACHING TOUT' |
| 4040 | 1398 CALL XERRWD(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1399 GO TO 690 |
| 1400 C | |
| 1401 C TOO MUCH ACCURACY FOR MACHINE PRECISION | |
| 1402 620 MSG = 'DASRT-- AT T (=R1) TOO MUCH ACCURACY REQUESTED' | |
| 4040 | 1403 CALL XERRWD(MSG,47,620,0,0,0,0,1,TN,0.0D0) |
| 4000 | 1404 MSG = 'DASRT-- FOR PRECISION OF MACHINE. RTOL AND ATOL' |
| 4040 | 1405 CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1406 MSG = 'DASRT-- WERE INCREASED TO APPROPRIATE VALUES' |
| 4040 | 1407 CALL XERRWD(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1408 C |
| 1409 GO TO 690 | |
| 1410 C WT(I) .LE. 0.0D0 FOR SOME I (NOT AT START OF PROBLEM) | |
| 1411 630 MSG = 'DASRT-- AT T (=R1) SOME ELEMENT OF WT' | |
| 4040 | 1412 CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0) |
| 4000 | 1413 MSG = 'DASRT-- HAS BECOME .LE. 0.0' |
| 4040 | 1414 CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1415 GO TO 690 |
| 1416 C | |
| 1417 C ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN | |
| 1418 640 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE' | |
| 4040 | 1419 CALL XERRWD(MSG,44,640,0,0,0,0,2,TN,H) |
| 4000 | 1420 MSG='DASRT-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN' |
| 4040 | 1421 CALL XERRWD(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1422 GO TO 690 |
| 1423 C | |
| 1424 C CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN | |
| 1425 650 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE' | |
| 4040 | 1426 CALL XERRWD(MSG,44,650,0,0,0,0,2,TN,H) |
| 4000 | 1427 MSG = 'DASRT-- CORRECTOR FAILED TO CONVERGE REPEATEDLY' |
| 4040 | 1428 CALL XERRWD(MSG,48,651,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1429 MSG = 'DASRT-- OR WITH ABS(H)=HMIN' |
| 4040 | 1430 CALL XERRWD(MSG,28,652,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1431 GO TO 690 |
| 1432 C | |
| 1433 C THE ITERATION MATRIX IS SINGULAR | |
| 1434 660 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE' | |
| 4040 | 1435 CALL XERRWD(MSG,44,660,0,0,0,0,2,TN,H) |
| 4000 | 1436 MSG = 'DASRT-- ITERATION MATRIX IS SINGULAR' |
| 4040 | 1437 CALL XERRWD(MSG,37,661,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1438 GO TO 690 |
| 1439 C | |
| 1440 C CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES. | |
| 1441 670 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE' | |
| 4040 | 1442 CALL XERRWD(MSG,44,670,0,0,0,0,2,TN,H) |
| 4000 | 1443 MSG = 'DASRT-- CORRECTOR COULD NOT CONVERGE. ALSO, THE' |
| 4040 | 1444 CALL XERRWD(MSG,49,671,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1445 MSG = 'DASRT-- ERROR TEST FAILED REPEATEDLY.' |
| 4040 | 1446 CALL XERRWD(MSG,38,672,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1447 GO TO 690 |
| 1448 C | |
| 1449 C CORRECTOR FAILURE BECAUSE IRES = -1 | |
| 1450 675 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE' | |
| 4040 | 1451 CALL XERRWD(MSG,44,675,0,0,0,0,2,TN,H) |
| 4000 | 1452 MSG = 'DASRT-- CORRECTOR COULD NOT CONVERGE BECAUSE' |
| 4040 | 1453 CALL XERRWD(MSG,45,676,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1454 MSG = 'DASRT-- IRES WAS EQUAL TO MINUS ONE' |
| 4040 | 1455 CALL XERRWD(MSG,36,677,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1456 GO TO 690 |
| 1457 C | |
| 1458 C FAILURE BECAUSE IRES = -2 | |
| 1459 680 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2)' | |
| 4040 | 1460 CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H) |
| 4000 | 1461 MSG = 'DASRT-- IRES WAS EQUAL TO MINUS TWO' |
| 4040 | 1462 CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1463 GO TO 690 |
| 1464 C | |
| 1465 C FAILED TO COMPUTE INITIAL YPRIME | |
| 1466 685 MSG = 'DASRT-- AT T (=R1) AND STEPSIZE H (=R2) THE' | |
| 4040 | 1467 CALL XERRWD(MSG,44,685,0,0,0,0,2,TN,HO) |
| 4000 | 1468 MSG = 'DASRT-- INITIAL YPRIME COULD NOT BE COMPUTED' |
| 4040 | 1469 CALL XERRWD(MSG,45,686,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1470 GO TO 690 |
| 1471 690 CONTINUE | |
| 1472 INFO(1)=-1 | |
| 1473 T=TN | |
| 1474 RWORK(LTN)=TN | |
| 1475 RWORK(LH)=H | |
| 1476 RETURN | |
| 1477 C----------------------------------------------------------------------- | |
| 1478 C THIS BLOCK HANDLES ALL ERROR RETURNS DUE | |
| 1479 C TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING | |
| 1480 C DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS | |
| 1481 C CALLED. IF THIS HAPPENS TWICE IN | |
| 1482 C SUCCESSION, EXECUTION IS TERMINATED | |
| 1483 C | |
| 1484 C----------------------------------------------------------------------- | |
| 1485 701 MSG = 'DASRT-- SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE' | |
| 4040 | 1486 CALL XERRWD(MSG,55,1,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1487 GO TO 750 |
| 1488 702 MSG = 'DASRT-- NEQ (=I1) .LE. 0' | |
| 4040 | 1489 CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0) |
| 4000 | 1490 GO TO 750 |
| 1491 703 MSG = 'DASRT-- MAXORD (=I1) NOT IN RANGE' | |
| 4040 | 1492 CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0) |
| 4000 | 1493 GO TO 750 |
| 1494 704 MSG='DASRT-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)' | |
| 4040 | 1495 CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0) |
| 4000 | 1496 GO TO 750 |
| 1497 705 MSG='DASRT-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)' | |
| 4040 | 1498 CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0) |
| 4000 | 1499 GO TO 750 |
| 1500 706 MSG = 'DASRT-- SOME ELEMENT OF RTOL IS .LT. 0' | |
| 4040 | 1501 CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1502 GO TO 750 |
| 1503 707 MSG = 'DASRT-- SOME ELEMENT OF ATOL IS .LT. 0' | |
| 4040 | 1504 CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1505 GO TO 750 |
| 1506 708 MSG = 'DASRT-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO' | |
| 4040 | 1507 CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1508 GO TO 750 |
| 1509 709 MSG='DASRT-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)' | |
| 4040 | 1510 CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT) |
| 4000 | 1511 GO TO 750 |
| 1512 710 MSG = 'DASRT-- HMAX (=R1) .LT. 0.0' | |
| 4040 | 1513 CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0) |
| 4000 | 1514 GO TO 750 |
| 1515 711 MSG = 'DASRT-- TOUT (=R1) BEHIND T (=R2)' | |
| 4040 | 1516 CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T) |
| 4000 | 1517 GO TO 750 |
| 1518 712 MSG = 'DASRT-- INFO(8)=1 AND H0=0.0' | |
| 4040 | 1519 CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1520 GO TO 750 |
| 1521 713 MSG = 'DASRT-- SOME ELEMENT OF WT IS .LE. 0.0' | |
| 4040 | 1522 CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1523 GO TO 750 |
| 1524 714 MSG='DASRT-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION' | |
| 4040 | 1525 CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T) |
| 4000 | 1526 GO TO 750 |
| 1527 715 MSG = 'DASRT-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)' | |
| 4040 | 1528 CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T) |
| 4000 | 1529 GO TO 750 |
| 1530 716 MSG = 'DASRT-- INFO(12)=1 AND MXSTP (=I1) .LT. 0' | |
| 4040 | 1531 CALL XERRWD(MSG,42,16,0,1,IWORK(LMXSTP),0,0,0.0D0,0.0D0) |
| 4000 | 1532 GO TO 750 |
| 1533 717 MSG = 'DASRT-- ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ' | |
| 4040 | 1534 CALL XERRWD(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0) |
| 4000 | 1535 GO TO 750 |
| 1536 718 MSG = 'DASRT-- MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ' | |
| 4040 | 1537 CALL XERRWD(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0) |
| 4000 | 1538 GO TO 750 |
| 1539 719 MSG = 'DASRT-- TOUT (=R1) IS EQUAL TO T (=R2)' | |
| 4040 | 1540 CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T) |
| 4000 | 1541 GO TO 750 |
| 1542 730 MSG = 'DASRT-- NG (=I1) .LT. 0' | |
| 4040 | 1543 CALL XERRWD(MSG,24,30,1,1,NG,0,0,0.0D0,0.0D0) |
| 4000 | 1544 GO TO 750 |
| 1545 732 MSG = 'DASRT-- ONE OR MORE COMPONENTS OF G HAS A ROOT' | |
| 4040 | 1546 CALL XERRWD(MSG,47,32,1,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1547 MSG = ' TOO NEAR TO THE INITIAL POINT' |
| 4040 | 1548 CALL XERRWD(MSG,38,32,1,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1549 750 IF(INFO(1).EQ.-1) GO TO 760 |
| 1550 INFO(1)=-1 | |
| 1551 IDID=-33 | |
| 1552 RETURN | |
| 1553 760 MSG = 'DASRT-- REPEATED OCCURRENCES OF ILLEGAL INPUT' | |
| 4040 | 1554 CALL XERRWD(MSG,46,801,0,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1555 770 MSG = 'DASRT-- RUN TERMINATED. APPARENT INFINITE LOOP' |
| 4040 | 1556 CALL XERRWD(MSG,47,802,1,0,0,0,0,0.0D0,0.0D0) |
| 4000 | 1557 RETURN |
| 1558 C-----------END OF SUBROUTINE DDASRT------------------------------------ | |
| 1559 END |