# HG changeset patch # User jwe # Date 1026848213 0 # Node ID 53b4eab68976aa6553c9aa787248d482ee8af243 # Parent 48d2bc4a3729ce094810ec0453b5a7f33e7782ff [project @ 2002-07-16 19:36:52 by jwe] diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,5 +1,8 @@ 2002-07-16 John W. Eaton + * DASRT.h (DASRT::set_ng, DASRT::get_ng): Delete + * DASRT.cc (DASRT::DASRT): Set ng here. + * DAEFunc.h (DAEFunc): Jacobian function now follows format of DASSL. * DASSL.cc (ddassl_j): Make it work. * DASPK.cc (ddaspk_j): Likewise. diff --git a/liboctave/DASRT.cc b/liboctave/DASRT.cc --- a/liboctave/DASRT.cc +++ b/liboctave/DASRT.cc @@ -163,13 +163,14 @@ info.resize (30, 0); npar = 0; + ng = 0; liw = 0; lrw = 0; } -DASRT::DASRT (const int& ng, const ColumnVector& state, - const ColumnVector& deriv, double time, DAERTFunc& f) +DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv, + double time, DAERTFunc& f) : DAERT (state, deriv, time, f) { n = size (); @@ -180,8 +181,6 @@ stop_time_set = false; stop_time = 0.0; - DAERTFunc::operator = (f); - sanity_checked = false; info.resize (30, 0); @@ -189,6 +188,16 @@ npar = 0; + DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function (); + + if (tmp_csub) + { + ColumnVector tmp = tmp_csub (state, time); + ng = tmp.length (); + } + else + ng = 0; + rpar.resize (npar+1); ipar.resize (npar+1); @@ -226,7 +235,6 @@ 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, @@ -234,7 +242,6 @@ int iwadd = iwork(18); - if (iwadd > 0) liw += iwadd; @@ -251,14 +258,12 @@ int rwadd = iwork(19); - if (rwadd > 0) lrw += rwadd; rwork.resize (lrw, 0.0); info(0) = info_zero; - } void @@ -276,17 +281,6 @@ } void -DASRT::set_ng (int the_ng) -{ - ng = the_ng; -} - -int DASRT::get_ng (void) -{ - return ng; -} - -void DASRT::clear_stop_time (void) { stop_time_set = false; @@ -302,10 +296,10 @@ info(0) = 0; for (int i = 0; i < n; i++) - { - y(i,0) = x(i); - ydot(i,0) = xdot(i); - } + { + y(i,0) = x(i); + ydot(i,0) = xdot(i); + } integration_error = false; @@ -335,13 +329,9 @@ sanity_checked = true; } - init_work_size (info(0)); - - - if (iwork.length () != liw) iwork.resize (liw); @@ -351,7 +341,6 @@ abs_tol = absolute_tolerance (); rel_tol = relative_tolerance (); - if (initial_step_size () >= 0.0) { rwork(2) = initial_step_size (); @@ -376,7 +365,6 @@ else info(6) = 0; - py = y.fortran_vec (); pydot = ydot.fortran_vec (); pinfo = info.fortran_vec (); @@ -404,9 +392,6 @@ info(3) = 0; } - - - F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, py, pydot, tout, pinfo, &rel_tol, &abs_tol, idid, prwork, lrw, piwork, liw, prpar, pipar, ddasrt_j, @@ -486,7 +471,6 @@ int n_out = tout.capacity (); - if (n_out > 0 && n > 0) { x_out.resize (n_out, n); @@ -501,25 +485,26 @@ retval = DASRT_result (x_out, xdot_out, t_out); return retval; } + if (idid == 4) t_out(j) = t; else t_out(j) = tout(j); - for (int i = 0; i < n; i++) { x_out(j,i) = y(i,0); xdot_out(j,i) = ydot(i,0); } - 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); - } + + 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); + } } } @@ -606,6 +591,7 @@ retval = DASRT_result (x_out, xdot_out, t_outs); return retval; } + if (idid == 4) t_out = t; @@ -624,7 +610,6 @@ t_outs.resize (i_out); i_out = n_out; } - } if (do_restart) diff --git a/liboctave/DASRT.h b/liboctave/DASRT.h --- a/liboctave/DASRT.h +++ b/liboctave/DASRT.h @@ -158,7 +158,7 @@ DASRT (void); - DASRT (const int& ng, const ColumnVector& x, const ColumnVector& xdot, + DASRT (const ColumnVector& state, const ColumnVector& deriv, double time, DAERTFunc& f); ~DASRT (void) { } @@ -167,8 +167,6 @@ void set_stop_time (double t); void clear_stop_time (void); - void set_ng (int the_ng); - int get_ng (void); DASRT_result integrate (const ColumnVector& tout); diff --git a/liboctave/DASSL.h b/liboctave/DASSL.h --- a/liboctave/DASSL.h +++ b/liboctave/DASSL.h @@ -113,9 +113,9 @@ DASSL (void); - DASSL (const ColumnVector& x, double time, DAEFunc& f); + DASSL (const ColumnVector& state, double time, DAEFunc& f); - DASSL (const ColumnVector& x, const ColumnVector& xdot, + DASSL (const ColumnVector& state, const ColumnVector& xdot, double time, DAEFunc& f); ~DASSL (void) { } diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,5 +1,7 @@ 2002-07-16 John W. Eaton + * DLD-FUNCTIONS/dasrt.cc (Fdasrt): No need to find ng here. + * DLD-FUNCTIONS/dassl.cc (dassl_user_jacobian): New function. (Fdassl): Handle Jacobian function. diff --git a/src/DLD-FUNCTIONS/dasrt.cc b/src/DLD-FUNCTIONS/dasrt.cc --- a/src/DLD-FUNCTIONS/dasrt.cc +++ b/src/DLD-FUNCTIONS/dasrt.cc @@ -327,7 +327,6 @@ DASRT_ABORT1 ("invalid recursive call"); int argp = 0; - int ng = 0; int nargin = args.length (); @@ -393,76 +392,25 @@ DAERTFunc func (dasrt_user_f); argp++; - double t0; - ColumnVector x0; if (args(1).is_string ()) { dasrt_cf = is_valid_function (args(1), "dasrt", true); if (! dasrt_cf) - DASRT_ABORT1 ("expecting function name as argument 2"); - else - { - octave_value_list blah; - octave_value_list inputs; - ColumnVector the_ts; - x0 = ColumnVector (args(2).vector_value ()); - - if (error_state) - DASRT_ABORT1 ("bad x0"); - - the_ts = ColumnVector (args(4).vector_value ()); - - if (error_state) - DASRT_ABORT1 ("bad tout"); - - t0 = the_ts(0); + DASRT_ABORT1 ("expecting function name as argument 2"); - inputs(0) = x0; - inputs(1) = t0; - - string g_name = args(1).string_value(); - - if (error_state) - DASRT_ABORT1 ("bad function name of stopping condition"); - - blah = feval (g_name, inputs, 2); - ColumnVector gout = ColumnVector (blah(0).vector_value ()); - - ng = gout.length(); - int testg; - for (testg = 0; testg < ng; testg++) - { - if (gout(testg) == 0) - DASRT_ABORT1 ("stopping condition satisfied at initial point"); - } - } argp++; func.set_constraint_function (dasrt_user_cf); } - else - { - // Now this second argument is not a string. It has to be x0. - // we call the dummy g function and set ng = 1; - - x0 = ColumnVector (args(1).vector_value ()); - - if (error_state) - DASRT_ABORT1 ("bad x0"); - - func.set_constraint_function (dasrt_dumb_cf); - - ng = 1; - } ColumnVector state (args(argp++).vector_value ()); if (error_state) DASRT_ABORT2 ("expecting state vector as argument %d", argp); - ColumnVector stateprime (args(argp++).vector_value()); + ColumnVector stateprime (args(argp++).vector_value ()); if (error_state) DASRT_ABORT2 @@ -477,12 +425,11 @@ double tzero = old_out_times (0); int ol = old_out_times.length (); - int ijk; ColumnVector out_times (ol-1, 0.0); - for (ijk = 1; ijk < ol; ijk++) - out_times(ijk-1) = old_out_times(ijk); + for (int i = 1; i < ol; i++) + out_times(i-1) = old_out_times(i); ColumnVector crit_times; @@ -506,10 +453,9 @@ DASRT_result output; - DASRT dae = DASRT (ng, state, stateprime, tzero, func); + DASRT dae = DASRT (state, stateprime, tzero, func); dae.copy (dasrt_opts); - dae.set_ng (ng); if (error_state) DASRT_ABORT1 ("something is wrong"); @@ -527,7 +473,7 @@ Matrix old_output_state = output.state (); int lstuff = old_output_times.length (); - int lstate = x0.length (); + int lstate = state.length (); ColumnVector output_times (lstuff+1, 0.0); @@ -536,20 +482,20 @@ output_times(0) = tzero; - for (ijk = 0; ijk < lstate; ijk++) + for (int i = 0; i < lstate; i++) { - output_deriv(0,ijk) = stateprime(ijk); - output_state(0,ijk) = state(ijk); + output_deriv(0,i) = stateprime(i); + output_state(0,i) = state(i); } - for (ijk = 0; ijk < lstuff; ijk++) + for (int i = 0; i < lstuff; i++) { - output_times(ijk+1) = old_output_times(ijk); + output_times(i+1) = old_output_times(i); - for (int lmnop=0; lmnop < lstate; lmnop++) + for (int j = 0; j < lstate; j++) { - output_deriv(ijk+1,lmnop) = old_output_deriv(ijk,lmnop); - output_state(ijk+1,lmnop) = old_output_state(ijk,lmnop); + output_deriv(i+1,j) = old_output_deriv(i,j); + output_state(i+1,j) = old_output_state(i,j); } }