Mercurial > hg > octave-max
changeset 3994:a41827ec5677
[project @ 2002-07-16 23:57:09 by jwe]
| author | jwe |
|---|---|
| date | Tue, 16 Jul 2002 23:57:09 +0000 |
| parents | f23bc69132cc |
| children | ee0304212be0 |
| files | liboctave/DAERT.h liboctave/DASRT.cc liboctave/DASRT.h liboctave/DASSL.cc src/DLD-FUNCTIONS/dasrt.cc |
| diffstat | 5 files changed, 103 insertions(+), 199 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/DAERT.h +++ b/liboctave/DAERT.h @@ -35,6 +35,9 @@ DAERT (void) : base_diff_alg_eqn (), DAERTFunc () { } + DAERT (const ColumnVector& x, double t, DAERTFunc& f) + : base_diff_alg_eqn (x, t), DAERTFunc (f) { } + DAERT (const ColumnVector& x, const ColumnVector& xdot, double t, DAERTFunc& f) : base_diff_alg_eqn (x, xdot, t), DAERTFunc (f) { }
--- a/liboctave/DASRT.cc +++ b/liboctave/DASRT.cc @@ -88,7 +88,7 @@ for (int i = 0; i < nn; i++) tmp_deriv(i) = deriv[i]; - ColumnVector tmp_fval = user_fsub (tmp_state, tmp_deriv, t, ires); + ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires); if (tmp_fval.length () == 0) ires = -2; @@ -101,12 +101,6 @@ return 0; } -//typedef int (*efptr) (const double& t, const int& n, double *state, -// double *ework, double *rpar, int *ipar, -// const int& ieform, int& ires); - -//static efptr e_fun; - int ddasrt_j (const double& time, const double *state, const double *deriv, double *pd, const double& cj, double *, int *) @@ -122,7 +116,7 @@ tmp_state.elem (i) = state [i]; } - Matrix tmp_pd = user_jsub (tmp_state, tmp_deriv, time, cj); + Matrix tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++) @@ -141,7 +135,7 @@ for (int i = 0; i < n; i++) tmp_state(i) = state[i]; - ColumnVector tmp_fval = user_csub (tmp_state, t); + ColumnVector tmp_fval = (*user_csub) (tmp_state, t); for (int i = 0; i < ng; i++) gout[i] = tmp_fval(i); @@ -149,7 +143,6 @@ return 0; } - DASRT::DASRT (void) : DAERT () { @@ -163,13 +156,43 @@ info.resize (30, 0); - npar = 0; ng = 0; liw = 0; lrw = 0; } +DASRT::DASRT (const ColumnVector& state, double time, DAERTFunc& f) + : DAERT (state, time, f) +{ + n = size (); + + initialized = false; + restart = false; + + stop_time_set = false; + stop_time = 0.0; + + 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) @@ -185,9 +208,6 @@ sanity_checked = false; info.resize (30, 0); - jroot.resize (ng, 1); - - npar = 0; DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); @@ -199,72 +219,10 @@ else ng = 0; - rpar.resize (npar+1); - ipar.resize (npar+1); - - info(11) = npar; - - // Also store it here, for communication with user-supplied - // subroutines. - ipar(0) = npar; - - y.resize (n, 1, 0.0); - ydot.resize (n, 1, 0.0); -} - -void -DASRT::init_work_size (int info_zero) -{ - double t; - double *py = y.fortran_vec (); - double *pydot = ydot.fortran_vec (); - double rel_tol = relative_tolerance (); - double abs_tol = absolute_tolerance (); - int *pinfo = info.fortran_vec (); - double *prpar = rpar.fortran_vec (); - int *pipar = ipar.fortran_vec (); - int *pjroot = jroot.fortran_vec (); - int idid; - - // We do not have to lie. - rwork.resize (5000+9*n+n*n, 0.0); - iwork.resize (n+20, 0); + liw = 20 + n + 1000; + lrw = 50 + 9*n + n*n + 1000; - liw = n+20; - lrw = 5000+9*n+n*n; - - double *prwork = rwork.fortran_vec (); - int *piwork = iwork.fortran_vec (); - - F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo, - &rel_tol, &abs_tol, idid, prwork, lrw, - piwork, liw, prpar, pipar, ddasrt_j, - ddasrt_g, ng, pjroot); - - int iwadd = iwork(18); - - if (iwadd > 0) - liw += iwadd; - - info(0) = 0; - - iwork.resize (liw, 0); - - piwork = iwork.fortran_vec (); - - F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo, - &rel_tol, &abs_tol, idid, prwork, lrw, - piwork, liw, prpar, pipar, ddasrt_j, - ddasrt_g, ng, pjroot); - - int rwadd = iwork(19); - - if (rwadd > 0) - lrw += rwadd; - - rwork.resize (lrw, 0.0); - - info(0) = info_zero; + jroot.resize (ng, 1); } void @@ -296,16 +254,8 @@ { info(0) = 0; - for (int i = 0; i < n; i++) - { - y(i,0) = x(i); - ydot(i,0) = xdot(i); - } - integration_error = false; - nn = n; - user_fsub = DAEFunc::function (); user_jsub = DAEFunc::jacobian_function (); user_csub = DAERTFunc::constraint_function (); @@ -315,11 +265,16 @@ else info(4) = 0; + px = x.fortran_vec (); + pxdot = xdot.fortran_vec (); + + nn = n; + if (! sanity_checked) { int ires = 0; - ColumnVector fval = user_fsub (x, xdot, t, ires); + ColumnVector fval = (*user_fsub) (x, xdot, t, ires); if (fval.length () != x.length ()) { @@ -333,8 +288,6 @@ sanity_checked = true; } - init_work_size (info(0)); - if (iwork.length () != liw) iwork.resize (liw); @@ -368,13 +321,9 @@ else info(6) = 0; - py = y.fortran_vec (); - pydot = ydot.fortran_vec (); pinfo = info.fortran_vec (); piwork = iwork.fortran_vec (); prwork = rwork.fortran_vec (); - prpar = rpar.fortran_vec (); - pipar = ipar.fortran_vec (); pjroot = jroot.fortran_vec (); info(5) = 0; @@ -395,9 +344,12 @@ info(3) = 0; } - F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, py, pydot, tout, pinfo, + double *dummy = 0; + int *idummy = 0; + + F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, px, pxdot, tout, pinfo, &rel_tol, &abs_tol, idid, prwork, lrw, - piwork, liw, prpar, pipar, ddasrt_j, + piwork, liw, dummy, idummy, ddasrt_j, ddasrt_g, ng, pjroot)); if (f77_exception_encountered) @@ -420,14 +372,10 @@ case 5: // The integration to TSTOP was successfully completed // (T=TSTOP) by stepping to TSTOP within the // tolerance. Must restart to continue. - for (int i = 0; i < n; i++) - x(i) = y(i,0); t = tout; break; case 4: // We've hit the stopping condition. - for (int i = 0; i < n; i++) - x(i) = y(i,0); break; case -1: // A large amount of work has been expended. (~500 steps). @@ -468,7 +416,7 @@ Matrix x_out; Matrix xdot_out; - ColumnVector t_out; + ColumnVector t_out = tout; int oldj = 0; @@ -478,11 +426,17 @@ { x_out.resize (n_out, n); xdot_out.resize (n_out, n); - t_out.resize (n_out); - for (int j = 0; j < n_out; j++) + for (int i = 0; i < n; i++) + { + x_out(0,i) = x(i); + xdot_out(0,i) = xdot(i); + } + + for (int j = 1; j < n_out; j++) { integrate (tout(j)); + if (integration_error) { retval = DASRT_result (x_out, xdot_out, t_out); @@ -496,17 +450,16 @@ for (int i = 0; i < n; i++) { - x_out(j,i) = y(i,0); - xdot_out(j,i) = ydot(i,0); + x_out(j,i) = x(i); + xdot_out(j,i) = xdot(i); } if (idid == 4) { - oldj = j; - j = n_out; - x_out.resize (oldj+1, n); - xdot_out.resize (oldj+1, n); - t_out.resize (oldj+1); + x_out.resize (j+1, n); + xdot_out.resize (j+1, n); + t_out.resize (j+1); + break; } } } @@ -523,7 +476,7 @@ Matrix x_out; Matrix xdot_out; - ColumnVector t_outs; + ColumnVector t_outs = tout; int n_out = tout.capacity (); @@ -531,14 +484,13 @@ { x_out.resize (n_out, n); xdot_out.resize (n_out, n); - t_outs.resize (n_out); int n_crit = tcrit.capacity (); if (n_crit > 0) { int i_crit = 0; - int i_out = 0; + int i_out = 1; double next_crit = tcrit(0); double next_out; while (i_out < n_out) @@ -602,11 +554,13 @@ { for (int i = 0; i < n; i++) { - x_out(i_out-1,i) = y(i,0); - xdot_out(i_out-1,i) = ydot(i,0); + x_out(i_out-1,i) = x(i); + xdot_out(i_out-1,i) = xdot(i); } + t_outs(i_out-1) = t_out; - if (idid ==4) + + if (idid == 4) { x_out.resize (i_out, n); xdot_out.resize (i_out, n);
--- a/liboctave/DASRT.h +++ b/liboctave/DASRT.h @@ -158,6 +158,8 @@ DASRT (void); + DASRT (const ColumnVector& state, double time, DAERTFunc& f); + DASRT (const ColumnVector& state, const ColumnVector& deriv, double time, DAERTFunc& f); @@ -189,38 +191,26 @@ int liw; int lrw; int idid; - int ieform; - int lun; int n; - int npar; int ng; Array<int> info; Array<int> iwork; - Array<int> ipar; Array<int> jroot; Array<double> rwork; - Array<double> rpar; - - Matrix y; - Matrix ydot; double abs_tol; double rel_tol; - double *py; - double *pydot; + double *px; + double *pxdot; int *pinfo; int *piwork; double *prwork; - double *prpar; - int *pipar; int *pjroot; - void init_work_size (int); - void integrate (double t); };
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -56,13 +56,13 @@ DASSL::DASSL (void) : DAE () { - stop_time_set = 0; + stop_time_set = false; stop_time = 0.0; liw = 0; lrw = 0; - sanity_checked = 0; + sanity_checked = false; info.resize (15); @@ -75,27 +75,24 @@ { n = size (); - stop_time_set = 0; + stop_time_set = false; stop_time = 0.0; liw = 20 + n; lrw = 40 + 9*n + n*n; - sanity_checked = 0; + sanity_checked = false; - info.resize (15); - - for (int i = 0; i < 15; i++) - info.elem (i) = 0; + info.resize (15, 0); } DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv, - double time, DAEFunc& f) + double time, DAEFunc& f) : DAE (state, deriv, time, f) { n = size (); - stop_time_set = 0; + stop_time_set = false; stop_time = 0.0; DAEFunc::set_function (f.function ()); @@ -104,7 +101,7 @@ liw = 20 + n; lrw = 40 + 9*n + n*n; - sanity_checked = 0; + sanity_checked = false; info.resize (15); @@ -116,20 +113,20 @@ DASSL::force_restart (void) { restart = 1; - integration_error = 0; + integration_error = false; } void DASSL::set_stop_time (double tt) { - stop_time_set = 1; + stop_time_set = true; stop_time = tt; } void DASSL::clear_stop_time (void) { - stop_time_set = 0; + stop_time_set = false; } int @@ -205,21 +202,21 @@ if (rwork.length () != lrw) rwork.resize (lrw); - integration_error = 0; + integration_error = false; - if (DAEFunc::jacobian_function ()) + user_fun = DAEFunc::fun; + user_jac = DAEFunc::jac; + + if (user_jac) info.elem (4) = 1; else info.elem (4) = 0; - double *px = x.fortran_vec (); + double *px = x.fortran_vec (); double *pxdot = xdot.fortran_vec (); nn = n; - user_fun = DAEFunc::fun; - user_jac = DAEFunc::jac; - if (! sanity_checked) { int ires = 0; @@ -231,11 +228,11 @@ (*current_liboctave_error_handler) ("dassl: inconsistent sizes for state and residual vectors"); - integration_error = 1; + integration_error = true; return retval; } - sanity_checked = 1; + sanity_checked = true; } if (stop_time_set) @@ -280,7 +277,7 @@ if (f77_exception_encountered) { - integration_error = 1; + integration_error = true; (*current_liboctave_error_handler) ("unrecoverable error in dassl"); } else @@ -294,7 +291,6 @@ case 3: // The integration to TOUT was successfully completed // (T=TOUT) by stepping past TOUT. Y(*) is obtained by // interpolation. YPRIME(*) is obtained by interpolation. - retval = x; t = tout; break; @@ -322,7 +318,7 @@ // and control is returned to the calling program. For // example, this occurs when invalid input is detected. default: - integration_error = 1; + integration_error = true; break; } }
--- a/src/DLD-FUNCTIONS/dasrt.cc +++ b/src/DLD-FUNCTIONS/dasrt.cc @@ -410,20 +410,13 @@ DASRT_ABORT2 ("expecting time derivative of state vector as argument %d", argp); - ColumnVector old_out_times (args(argp++).vector_value ()); + ColumnVector out_times (args(argp++).vector_value ()); if (error_state) DASRT_ABORT2 ("expecting output time vector as %s argument %d", argp); - double tzero = old_out_times (0); - - int ol = old_out_times.length (); - - ColumnVector out_times (ol-1, 0.0); - - for (int i = 1; i < ol; i++) - out_times(i-1) = old_out_times(i); + double tzero = out_times (0); ColumnVector crit_times; @@ -459,41 +452,9 @@ if (! error_state) { - ColumnVector old_output_times = output.times (); - - Matrix old_output_deriv = output.deriv (); - Matrix old_output_state = output.state (); - - int lstuff = old_output_times.length (); - int lstate = state.length (); - - ColumnVector output_times (lstuff+1, 0.0); - - Matrix output_deriv (lstuff+1, lstate, 0.0); - Matrix output_state (lstuff+1, lstate, 0.0); - - output_times(0) = tzero; - - for (int i = 0; i < lstate; i++) - { - output_deriv(0,i) = stateprime(i); - output_state(0,i) = state(i); - } - - for (int i = 0; i < lstuff; i++) - { - output_times(i+1) = old_output_times(i); - - for (int j = 0; j < lstate; j++) - { - output_deriv(i+1,j) = old_output_deriv(i,j); - output_state(i+1,j) = old_output_state(i,j); - } - } - - retval(2) = output_times; - retval(1) = output_deriv; - retval(0) = output_state; + retval(2) = output.times (); + retval(1) = output.deriv (); + retval(0) = output.state (); } else {