Mercurial > hg > octave-lyh
annotate liboctave/LSODE.cc @ 8807:401d54a83690
use 'invalid', not 'illegal'
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 18 Feb 2009 14:32:53 -0500 |
parents | 29980c6b8604 |
children | eb63fbe60fab |
rev | line source |
---|---|
3 | 1 /* |
2 | |
7017 | 3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2002, |
4 2003, 2004, 2005, 2006, 2007 John W. Eaton | |
3 | 5 |
6 This file is part of Octave. | |
7 | |
8 Octave is free software; you can redistribute it and/or modify it | |
9 under the terms of the GNU General Public License as published by the | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
3 | 12 |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
17 | |
18 You should have received a copy of the GNU General Public License | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
3 | 21 |
22 */ | |
23 | |
238 | 24 #ifdef HAVE_CONFIG_H |
1192 | 25 #include <config.h> |
3 | 26 #endif |
27 | |
1367 | 28 #include <cfloat> |
29 | |
5765 | 30 #include <sstream> |
31 | |
1842 | 32 #include "LSODE.h" |
1847 | 33 #include "f77-fcn.h" |
227 | 34 #include "lo-error.h" |
7231 | 35 #include "lo-math.h" |
4180 | 36 #include "quit.h" |
3 | 37 |
5275 | 38 typedef octave_idx_type (*lsode_fcn_ptr) (const octave_idx_type&, const double&, double*, |
39 double*, octave_idx_type&); | |
3507 | 40 |
5275 | 41 typedef octave_idx_type (*lsode_jac_ptr) (const octave_idx_type&, const double&, double*, |
42 const octave_idx_type&, const octave_idx_type&, double*, const | |
43 octave_idx_type&); | |
3507 | 44 |
3 | 45 extern "C" |
4552 | 46 { |
47 F77_RET_T | |
5275 | 48 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, octave_idx_type&, double*, double&, |
49 double&, octave_idx_type&, double&, const double*, octave_idx_type&, | |
50 octave_idx_type&, octave_idx_type&, double*, octave_idx_type&, octave_idx_type*, octave_idx_type&, | |
51 lsode_jac_ptr, octave_idx_type&); | |
4552 | 52 } |
3 | 53 |
532 | 54 static ODEFunc::ODERHSFunc user_fun; |
55 static ODEFunc::ODEJacFunc user_jac; | |
3 | 56 static ColumnVector *tmp_x; |
57 | |
5275 | 58 static octave_idx_type |
59 lsode_f (const octave_idx_type& neq, const double& time, double *, | |
60 double *deriv, octave_idx_type& ierr) | |
3 | 61 { |
4180 | 62 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
63 | |
2343 | 64 ColumnVector tmp_deriv; |
3 | 65 |
1360 | 66 // NOTE: this won't work if LSODE passes copies of the state vector. |
67 // In that case we have to create a temporary vector object | |
68 // and copy. | |
69 | |
1251 | 70 tmp_deriv = (*user_fun) (*tmp_x, time); |
3 | 71 |
258 | 72 if (tmp_deriv.length () == 0) |
1251 | 73 ierr = -1; |
258 | 74 else |
75 { | |
5275 | 76 for (octave_idx_type i = 0; i < neq; i++) |
258 | 77 deriv [i] = tmp_deriv.elem (i); |
78 } | |
3 | 79 |
4180 | 80 END_INTERRUPT_WITH_EXCEPTIONS; |
81 | |
3 | 82 return 0; |
83 } | |
84 | |
5275 | 85 static octave_idx_type |
86 lsode_j (const octave_idx_type& neq, const double& time, double *, | |
87 const octave_idx_type&, const octave_idx_type&, double *pd, const octave_idx_type& nrowpd) | |
3 | 88 { |
4180 | 89 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
90 | |
1251 | 91 Matrix tmp_jac (neq, neq); |
3 | 92 |
1360 | 93 // NOTE: this won't work if LSODE passes copies of the state vector. |
94 // In that case we have to create a temporary vector object | |
95 // and copy. | |
96 | |
1251 | 97 tmp_jac = (*user_jac) (*tmp_x, time); |
3 | 98 |
5275 | 99 for (octave_idx_type j = 0; j < neq; j++) |
100 for (octave_idx_type i = 0; i < neq; i++) | |
1251 | 101 pd [nrowpd * j + i] = tmp_jac (i, j); |
3 | 102 |
4180 | 103 END_INTERRUPT_WITH_EXCEPTIONS; |
104 | |
3 | 105 return 0; |
106 } | |
107 | |
108 ColumnVector | |
1842 | 109 LSODE::do_integrate (double tout) |
3 | 110 { |
1945 | 111 ColumnVector retval; |
112 | |
5275 | 113 static octave_idx_type nn = 0; |
4049 | 114 |
115 if (! initialized || restart || ODEFunc::reset || LSODE_options::reset) | |
1945 | 116 { |
4049 | 117 integration_error = false; |
1945 | 118 |
4049 | 119 initialized = true; |
120 | |
121 istate = 1; | |
122 | |
5275 | 123 octave_idx_type n = size (); |
4049 | 124 |
125 nn = n; | |
3955 | 126 |
5275 | 127 octave_idx_type max_maxord = 0; |
4231 | 128 |
4049 | 129 if (integration_method () == "stiff") |
130 { | |
4231 | 131 max_maxord = 5; |
132 | |
4049 | 133 if (jac) |
134 method_flag = 21; | |
135 else | |
136 method_flag = 22; | |
3955 | 137 |
5553 | 138 liw = 20 + n; |
139 lrw = 22 + n * (9 + n); | |
4049 | 140 } |
141 else | |
142 { | |
4231 | 143 max_maxord = 12; |
144 | |
4049 | 145 method_flag = 10; |
3955 | 146 |
4049 | 147 liw = 20; |
148 lrw = 22 + 16 * n; | |
149 } | |
150 | |
4231 | 151 maxord = maximum_order (); |
152 | |
5552 | 153 iwork.resize (liw); |
154 | |
155 for (octave_idx_type i = 4; i < 9; i++) | |
156 iwork(i) = 0; | |
157 | |
158 rwork.resize (lrw); | |
159 | |
160 for (octave_idx_type i = 4; i < 9; i++) | |
161 rwork(i) = 0; | |
162 | |
4231 | 163 if (maxord >= 0) |
164 { | |
165 if (maxord > 0 && maxord <= max_maxord) | |
166 { | |
167 iwork(4) = maxord; | |
168 iopt = 1; | |
169 } | |
170 else | |
171 { | |
172 (*current_liboctave_error_handler) | |
173 ("lsode: invalid value for maximum order"); | |
174 integration_error = true; | |
175 return retval; | |
176 } | |
177 } | |
178 | |
4049 | 179 if (stop_time_set) |
180 { | |
181 itask = 4; | |
182 rwork(0) = stop_time; | |
183 iopt = 1; | |
184 } | |
185 else | |
186 { | |
187 itask = 1; | |
188 } | |
258 | 189 |
4049 | 190 px = x.fortran_vec (); |
3 | 191 |
4049 | 192 piwork = iwork.fortran_vec (); |
193 prwork = rwork.fortran_vec (); | |
194 | |
195 restart = false; | |
196 | |
197 // ODEFunc | |
3 | 198 |
4049 | 199 // NOTE: this won't work if LSODE passes copies of the state vector. |
200 // In that case we have to create a temporary vector object | |
201 // and copy. | |
3 | 202 |
4049 | 203 tmp_x = &x; |
204 | |
205 user_fun = function (); | |
206 user_jac = jacobian_function (); | |
207 | |
2343 | 208 ColumnVector xdot = (*user_fun) (x, t); |
209 | |
210 if (x.length () != xdot.length ()) | |
211 { | |
212 (*current_liboctave_error_handler) | |
213 ("lsode: inconsistent sizes for state and derivative vectors"); | |
214 | |
3995 | 215 integration_error = true; |
2343 | 216 return retval; |
217 } | |
218 | |
4049 | 219 ODEFunc::reset = false; |
220 | |
221 // LSODE_options | |
222 | |
223 rel_tol = relative_tolerance (); | |
224 abs_tol = absolute_tolerance (); | |
225 | |
5275 | 226 octave_idx_type abs_tol_len = abs_tol.length (); |
4049 | 227 |
228 if (abs_tol_len == 1) | |
229 itol = 1; | |
230 else if (abs_tol_len == n) | |
231 itol = 2; | |
232 else | |
233 { | |
234 (*current_liboctave_error_handler) | |
235 ("lsode: inconsistent sizes for state and absolute tolerance vectors"); | |
236 | |
237 integration_error = true; | |
238 return retval; | |
239 } | |
2343 | 240 |
4049 | 241 double iss = initial_step_size (); |
242 if (iss >= 0.0) | |
243 { | |
244 rwork(4) = iss; | |
245 iopt = 1; | |
246 } | |
247 | |
248 double maxss = maximum_step_size (); | |
249 if (maxss >= 0.0) | |
250 { | |
251 rwork(5) = maxss; | |
252 iopt = 1; | |
253 } | |
254 | |
255 double minss = minimum_step_size (); | |
256 if (minss >= 0.0) | |
257 { | |
258 rwork(6) = minss; | |
259 iopt = 1; | |
260 } | |
261 | |
5275 | 262 octave_idx_type sl = step_limit (); |
4049 | 263 if (sl > 0) |
264 { | |
265 iwork(5) = sl; | |
266 iopt = 1; | |
267 } | |
268 | |
269 pabs_tol = abs_tol.fortran_vec (); | |
270 | |
271 LSODE_options::reset = false; | |
3 | 272 } |
273 | |
4583 | 274 F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol, |
275 pabs_tol, itask, istate, iopt, prwork, lrw, | |
276 piwork, liw, lsode_j, method_flag)); | |
3 | 277 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
278 switch (istate) |
3 | 279 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
280 case 1: // prior to initial integration step. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
281 case 2: // lsode was successful. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
282 retval = x; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
283 t = tout; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
284 break; |
3996 | 285 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
286 case -1: // excess work done on this call (perhaps wrong mf). |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
287 case -2: // excess accuracy requested (tolerances too small). |
8807
401d54a83690
use 'invalid', not 'illegal'
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
288 case -3: // invalid input detected (see printed message). |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
289 case -4: // repeated error test failures (check all inputs). |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
290 case -5: // repeated convergence failures (perhaps bad jacobian |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
291 // supplied or wrong choice of mf or tolerances). |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
292 case -6: // error weight became zero during problem. (solution |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
293 // component i vanished, and atol or atol(i) = 0.) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
294 case -13: // return requested in user-supplied function. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
295 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
296 break; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
297 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
298 default: |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
299 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
300 (*current_liboctave_error_handler) |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
301 ("unrecognized value of istate (= %d) returned from lsode", |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
302 istate); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
303 break; |
3 | 304 } |
305 | |
1945 | 306 return retval; |
3 | 307 } |
308 | |
3959 | 309 std::string |
310 LSODE::error_message (void) const | |
311 { | |
312 std::string retval; | |
313 | |
5765 | 314 std::ostringstream buf; |
315 buf << t; | |
316 std::string t_curr = buf.str (); | |
4042 | 317 |
3959 | 318 switch (istate) |
319 { | |
320 case 1: | |
3996 | 321 retval = "prior to initial integration step"; |
3959 | 322 break; |
323 | |
324 case 2: | |
325 retval = "successful exit"; | |
326 break; | |
327 | |
328 case 3: | |
329 retval = "prior to continuation call with modified parameters"; | |
330 break; | |
331 | |
3996 | 332 case -1: |
4042 | 333 retval = std::string ("excess work on this call (t = ") |
4043 | 334 + t_curr + "; perhaps wrong integration method)"; |
3996 | 335 break; |
336 | |
337 case -2: | |
338 retval = "excess accuracy requested (tolerances too small)"; | |
339 break; | |
340 | |
341 case -3: | |
342 retval = "invalid input detected (see printed message)"; | |
343 break; | |
344 | |
345 case -4: | |
4042 | 346 retval = std::string ("repeated error test failures (t = ") |
4043 | 347 + t_curr + "check all inputs)"; |
3996 | 348 break; |
349 | |
350 case -5: | |
4042 | 351 retval = std::string ("repeated convergence failures (t = ") |
4043 | 352 + t_curr |
4042 | 353 + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)"; |
3996 | 354 break; |
355 | |
356 case -6: | |
4042 | 357 retval = std::string ("error weight became zero during problem. (t = ") |
4043 | 358 + t_curr |
4042 | 359 + "; solution component i vanished, and atol or atol(i) == 0)"; |
3996 | 360 break; |
361 | |
362 case -13: | |
4042 | 363 retval = "return requested in user-supplied function (t = " |
4043 | 364 + t_curr + ")"; |
3996 | 365 break; |
366 | |
3959 | 367 default: |
368 retval = "unknown error state"; | |
369 break; | |
370 } | |
371 | |
372 return retval; | |
373 } | |
374 | |
3 | 375 Matrix |
1842 | 376 LSODE::do_integrate (const ColumnVector& tout) |
3 | 377 { |
378 Matrix retval; | |
4049 | 379 |
5275 | 380 octave_idx_type n_out = tout.capacity (); |
381 octave_idx_type n = size (); | |
3 | 382 |
383 if (n_out > 0 && n > 0) | |
384 { | |
385 retval.resize (n_out, n); | |
386 | |
5275 | 387 for (octave_idx_type i = 0; i < n; i++) |
3 | 388 retval.elem (0, i) = x.elem (i); |
389 | |
5275 | 390 for (octave_idx_type j = 1; j < n_out; j++) |
3 | 391 { |
1842 | 392 ColumnVector x_next = do_integrate (tout.elem (j)); |
258 | 393 |
394 if (integration_error) | |
395 return retval; | |
396 | |
5275 | 397 for (octave_idx_type i = 0; i < n; i++) |
3 | 398 retval.elem (j, i) = x_next.elem (i); |
399 } | |
400 } | |
401 | |
402 return retval; | |
403 } | |
404 | |
405 Matrix | |
3511 | 406 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit) |
3 | 407 { |
408 Matrix retval; | |
4049 | 409 |
5275 | 410 octave_idx_type n_out = tout.capacity (); |
411 octave_idx_type n = size (); | |
3 | 412 |
413 if (n_out > 0 && n > 0) | |
414 { | |
415 retval.resize (n_out, n); | |
416 | |
5275 | 417 for (octave_idx_type i = 0; i < n; i++) |
3 | 418 retval.elem (0, i) = x.elem (i); |
419 | |
5275 | 420 octave_idx_type n_crit = tcrit.capacity (); |
3 | 421 |
422 if (n_crit > 0) | |
423 { | |
5275 | 424 octave_idx_type i_crit = 0; |
425 octave_idx_type i_out = 1; | |
3 | 426 double next_crit = tcrit.elem (0); |
427 double next_out; | |
428 while (i_out < n_out) | |
429 { | |
3995 | 430 bool do_restart = false; |
3 | 431 |
432 next_out = tout.elem (i_out); | |
433 if (i_crit < n_crit) | |
434 next_crit = tcrit.elem (i_crit); | |
435 | |
5275 | 436 octave_idx_type save_output; |
3 | 437 double t_out; |
438 | |
439 if (next_crit == next_out) | |
440 { | |
441 set_stop_time (next_crit); | |
442 t_out = next_out; | |
443 save_output = 1; | |
444 i_out++; | |
445 i_crit++; | |
3995 | 446 do_restart = true; |
3 | 447 } |
448 else if (next_crit < next_out) | |
449 { | |
450 if (i_crit < n_crit) | |
451 { | |
452 set_stop_time (next_crit); | |
453 t_out = next_crit; | |
454 save_output = 0; | |
455 i_crit++; | |
3995 | 456 do_restart = true; |
3 | 457 } |
458 else | |
459 { | |
460 clear_stop_time (); | |
461 t_out = next_out; | |
462 save_output = 1; | |
463 i_out++; | |
464 } | |
465 } | |
466 else | |
467 { | |
468 set_stop_time (next_crit); | |
469 t_out = next_out; | |
470 save_output = 1; | |
471 i_out++; | |
472 } | |
473 | |
1842 | 474 ColumnVector x_next = do_integrate (t_out); |
3 | 475 |
258 | 476 if (integration_error) |
477 return retval; | |
478 | |
3 | 479 if (save_output) |
480 { | |
5275 | 481 for (octave_idx_type i = 0; i < n; i++) |
3 | 482 retval.elem (i_out-1, i) = x_next.elem (i); |
483 } | |
484 | |
485 if (do_restart) | |
486 force_restart (); | |
487 } | |
488 } | |
489 else | |
258 | 490 { |
1842 | 491 retval = do_integrate (tout); |
258 | 492 |
493 if (integration_error) | |
494 return retval; | |
495 } | |
3 | 496 } |
497 | |
498 return retval; | |
499 } | |
500 | |
501 /* | |
502 ;;; Local Variables: *** | |
503 ;;; mode: C++ *** | |
504 ;;; End: *** | |
505 */ |