Mercurial > hg > octave-nkf
diff liboctave/DASRT.cc @ 4049:a35a3c5d4740
[project @ 2002-08-16 08:54:31 by jwe]
author | jwe |
---|---|
date | Fri, 16 Aug 2002 08:54:31 +0000 |
parents | 6fae69a1796e |
children | 6481f41a79f3 |
line wrap: on
line diff
--- a/liboctave/DASRT.cc +++ b/liboctave/DASRT.cc @@ -48,10 +48,6 @@ #include "f77-fcn.h" #include "lo-error.h" -#ifndef F77_FUNC -#define F77_FUNC(x, X) F77_FCN (x, X) -#endif - typedef int (*dasrt_fcn_ptr) (const double&, const double*, const double*, double*, int&, double*, int*); @@ -142,102 +138,63 @@ return 0; } -DASRT::DASRT (void) - : DAERT () -{ - initialized = false; - - sanity_checked = false; - - info.resize (30, 0); - - ng = 0; - - liw = 0; - lrw = 0; -} - -DASRT::DASRT (const ColumnVector& state, double time, DAERTFunc& f) - : DAERT (state, time, f) -{ - n = size (); - - initialized = false; - - liw = 20 + n; - lrw = 50 + 9*n + n*n; - - sanity_checked = false; - - info.resize (15, 0); - - DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); - - if (tmp_csub) - { - ColumnVector tmp = tmp_csub (state, time); - ng = tmp.length (); - } - else - ng = 0; - - jroot.resize (ng, 1); -} - -DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv, - double time, DAERTFunc& f) - : DAERT (state, deriv, time, f) -{ - n = size (); - - initialized = false; - - sanity_checked = false; - - info.resize (30, 0); - - DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); - - if (tmp_csub) - { - ColumnVector tmp = tmp_csub (state, time); - ng = tmp.length (); - } - else - ng = 0; - - liw = 20 + n + 1000; - lrw = 50 + 9*n + n*n + 1000; - - jroot.resize (ng, 1); -} - void DASRT::integrate (double tout) { DASRT_result retval; - if (! initialized) + // I suppose this is the safe thing to do. If this is the first + // call, or if anything about the problem has changed, we should + // start completely fresh. + + if (! initialized || restart + || DAEFunc::reset || DAERTFunc::reset || DASRT_options::reset) { - info(0) = 0; - integration_error = false; - user_fsub = DAEFunc::function (); - user_jsub = DAEFunc::jacobian_function (); - user_csub = DAERTFunc::constraint_function (); + initialized = true; + + info.resize (15); + + for (int i = 0; i < 15; i++) + info(i) = 0; + + pinfo = info.fortran_vec (); + + int n = size (); + + nn = n; - if (user_jsub) - info(4) = 1; + liw = 20 + n; + lrw = 50 + 9*n + n*n; + + iwork.resize (liw); + rwork.resize (lrw); + + info(0) = 0; + + if (stop_time_set) + { + info(3) = 1; + rwork(0) = stop_time; + } else - info(4) = 0; + info(3) = 0; px = x.fortran_vec (); pxdot = xdot.fortran_vec (); - nn = n; + piwork = iwork.fortran_vec (); + prwork = rwork.fortran_vec (); + + restart = false; - if (! sanity_checked) + // DAEFunc + + user_fsub = DAEFunc::function (); + user_jsub = DAEFunc::jacobian_function (); + + if (user_fsub) { int ires = 0; @@ -246,20 +203,85 @@ if (fval.length () != x.length ()) { (*current_liboctave_error_handler) - ("dassl: inconsistent sizes for state and residual vectors"); + ("dasrt: inconsistent sizes for state and residual vectors"); integration_error = true; return; } + } + else + { + (*current_liboctave_error_handler) + ("dasrt: no user supplied RHS subroutine!"); - sanity_checked = true; + integration_error = true; + return; + } + + info(4) = user_jsub ? 1 : 0; + + DAEFunc::reset = false; + + // DAERTFunc + + user_csub = DAERTFunc::constraint_function (); + + if (user_csub) + { + ColumnVector tmp = (*user_csub) (x, t); + ng = tmp.length (); } - - if (iwork.length () != liw) - iwork.resize (liw); + else + ng = 0; + + jroot.resize (ng, 1); + + pjroot = jroot.fortran_vec (); + + DAERTFunc::reset = false; + + // DASRT_options + + if (maximum_step_size () >= 0.0) + { + rwork(1) = maximum_step_size (); + info(6) = 1; + } + else + info(6) = 0; - if (rwork.length () != lrw) - rwork.resize (lrw); + if (initial_step_size () >= 0.0) + { + rwork(2) = initial_step_size (); + info(7) = 1; + } + else + info(7) = 0; + + int maxord = maximum_order (); + if (maxord >= 0) + { + if (maxord > 0 && maxord < 6) + { + info(8) = 1; + iwork(2) = maxord; + } + else + { + (*current_liboctave_error_handler) + ("dassl: invalid value for maximum order"); + integration_error = true; + return; + } + } + + if (step_limit () >= 0) + { + info(11) = 1; + iwork(18) = step_limit (); + } + else + info(11) = 0; abs_tol = absolute_tolerance (); rel_tol = relative_tolerance (); @@ -278,65 +300,22 @@ else { (*current_liboctave_error_handler) - ("dassl: inconsistent sizes for tolerance arrays"); + ("dasrt: inconsistent sizes for tolerance arrays"); integration_error = true; return; } - if (initial_step_size () >= 0.0) - { - rwork(2) = initial_step_size (); - info(7) = 1; - } - else - info(7) = 0; - - if (step_limit () >= 0) - { - info(11) = 1; - iwork(18) = step_limit (); - } - else - info(11) = 0; - - if (maximum_step_size () >= 0.0) - { - rwork(1) = maximum_step_size (); - info(6) = 1; - } - else - info(6) = 0; - - pinfo = info.fortran_vec (); - piwork = iwork.fortran_vec (); pabs_tol = abs_tol.fortran_vec (); prel_tol = rel_tol.fortran_vec (); - prwork = rwork.fortran_vec (); - pjroot = jroot.fortran_vec (); - info(5) = 0; - info(8) = 0; - initialized = true; + DASRT_options::reset = false; } - if (restart) - { - info(0) = 0; + static double *dummy = 0; + static int *idummy = 0; - if (stop_time_set) - { - info(3) = 1; - rwork(0) = stop_time; - } - else - info(3) = 0; - } - - double *dummy = 0; - int *idummy = 0; - - F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, px, pxdot, tout, pinfo, + F77_XFCN (ddasrt, DASRT, (ddasrt_f, nn, t, px, pxdot, tout, pinfo, prel_tol, pabs_tol, istate, prwork, lrw, piwork, liw, dummy, idummy, ddasrt_j, ddasrt_g, ng, pjroot)); @@ -344,7 +323,7 @@ if (f77_exception_encountered) { integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in dassl"); + (*current_liboctave_error_handler) ("unrecoverable error in dasrt"); } else { @@ -409,6 +388,7 @@ ColumnVector t_out = tout; int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) { @@ -467,6 +447,7 @@ ColumnVector t_outs = tout; int n_out = tout.capacity (); + int n = size (); if (n_out > 0 && n > 0) {