comparison liboctave/LSODE.cc @ 2343:d7592de300ea

[project @ 1996-07-24 21:42:44 by jwe]
author jwe
date Wed, 24 Jul 1996 21:44:51 +0000
parents 1b57120c997b
children 8b262e771614
comparison
equal deleted inserted replaced
2342:95e511896bf5 2343:d7592de300ea
69 itask = 1; 69 itask = 1;
70 iopt = 0; 70 iopt = 0;
71 71
72 liw = 20 + n; 72 liw = 20 + n;
73 lrw = 22 + n * (9 + n); 73 lrw = 22 + n * (9 + n);
74
75 sanity_checked = 0;
74 } 76 }
75 77
76 LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f) 78 LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f)
77 : ODE (state, time, f), LSODE_options () 79 : ODE (state, time, f), LSODE_options ()
78 { 80 {
89 itask = 1; 91 itask = 1;
90 iopt = 0; 92 iopt = 0;
91 93
92 liw = 20 + n; 94 liw = 20 + n;
93 lrw = 22 + n * (9 + n); 95 lrw = 22 + n * (9 + n);
96
97 sanity_checked = 0;
94 } 98 }
95 99
96 void 100 void
97 LSODE::force_restart (void) 101 LSODE::force_restart (void)
98 { 102 {
114 118
115 int 119 int
116 lsode_f (const int& neq, const double& time, double *, 120 lsode_f (const int& neq, const double& time, double *,
117 double *deriv, int& ierr) 121 double *deriv, int& ierr)
118 { 122 {
119 ColumnVector tmp_deriv (neq); 123 ColumnVector tmp_deriv;
120 124
121 // NOTE: this won't work if LSODE passes copies of the state vector. 125 // NOTE: this won't work if LSODE passes copies of the state vector.
122 // In that case we have to create a temporary vector object 126 // In that case we have to create a temporary vector object
123 // and copy. 127 // and copy.
124 128
195 // and copy. 199 // and copy.
196 200
197 tmp_x = &x; 201 tmp_x = &x;
198 user_fun = function (); 202 user_fun = function ();
199 user_jac = jacobian_function (); 203 user_jac = jacobian_function ();
204
205 if (! sanity_checked)
206 {
207 ColumnVector xdot = (*user_fun) (x, t);
208
209 if (x.length () != xdot.length ())
210 {
211 (*current_liboctave_error_handler)
212 ("lsode: inconsistent sizes for state and derivative vectors");
213
214 integration_error = 1;
215 return retval;
216 }
217
218 sanity_checked = 1;
219 }
200 220
201 // Try 5000 steps before giving up. 221 // Try 5000 steps before giving up.
202 222
203 iwork.elem (5) = 5000; 223 iwork.elem (5) = 5000;
204 224