Mercurial > hg > octave-nkf

comparison liboctave/LSODE.cc @ 10314:07ebe522dac2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave C++ sources
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:23:32 -0500
parents 4c0cdbe0acca
children 12884915a8e4
comparison
equal deleted inserted replaced
10313:f3b65e1ae355 10314:07ebe522dac2
34 #include "lo-error.h" 34 #include "lo-error.h"
35 #include "lo-math.h" 35 #include "lo-math.h"
36 #include "quit.h" 36 #include "quit.h"
37 37
38 typedef octave_idx_type (*lsode_fcn_ptr) (const octave_idx_type&, const double&, double*, 38 typedef octave_idx_type (*lsode_fcn_ptr) (const octave_idx_type&, const double&, double*,
39 double*, octave_idx_type&); 39 double*, octave_idx_type&);
40 40
41 typedef octave_idx_type (*lsode_jac_ptr) (const octave_idx_type&, const double&, double*, 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 42 const octave_idx_type&, const octave_idx_type&, double*, const
43 octave_idx_type&); 43 octave_idx_type&);
44 44
45 extern "C" 45 extern "C"
46 { 46 {
47 F77_RET_T 47 F77_RET_T
48 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, octave_idx_type&, double*, double&, 48 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, octave_idx_type&, double*, double&,
49 double&, octave_idx_type&, double&, const double*, octave_idx_type&, 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&, 50 octave_idx_type&, octave_idx_type&, double*, octave_idx_type&, octave_idx_type*, octave_idx_type&,
51 lsode_jac_ptr, octave_idx_type&); 51 lsode_jac_ptr, octave_idx_type&);
52 } 52 }
53 53
54 static ODEFunc::ODERHSFunc user_fun; 54 static ODEFunc::ODERHSFunc user_fun;
55 static ODEFunc::ODEJacFunc user_jac; 55 static ODEFunc::ODEJacFunc user_jac;
56 static ColumnVector *tmp_x; 56 static ColumnVector *tmp_x;
57 57
58 static octave_idx_type 58 static octave_idx_type
59 lsode_f (const octave_idx_type& neq, const double& time, double *, 59 lsode_f (const octave_idx_type& neq, const double& time, double *,
60 double *deriv, octave_idx_type& ierr) 60 double *deriv, octave_idx_type& ierr)
61 { 61 {
62 BEGIN_INTERRUPT_WITH_EXCEPTIONS; 62 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
63 63
64 ColumnVector tmp_deriv; 64 ColumnVector tmp_deriv;
65 65
72 if (tmp_deriv.length () == 0) 72 if (tmp_deriv.length () == 0)
73 ierr = -1; 73 ierr = -1;
74 else 74 else
75 { 75 {
76 for (octave_idx_type i = 0; i < neq; i++) 76 for (octave_idx_type i = 0; i < neq; i++)
77 deriv [i] = tmp_deriv.elem (i); 77 deriv [i] = tmp_deriv.elem (i);
78 } 78 }
79 79
80 END_INTERRUPT_WITH_EXCEPTIONS; 80 END_INTERRUPT_WITH_EXCEPTIONS;
81 81
82 return 0; 82 return 0;
83 } 83 }
84 84
85 static octave_idx_type 85 static octave_idx_type
86 lsode_j (const octave_idx_type& neq, const double& time, double *, 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) 87 const octave_idx_type&, const octave_idx_type&, double *pd, const octave_idx_type& nrowpd)
88 { 88 {
89 BEGIN_INTERRUPT_WITH_EXCEPTIONS; 89 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
90 90
91 Matrix tmp_jac (neq, neq); 91 Matrix tmp_jac (neq, neq);
92 92
125 nn = n; 125 nn = n;
126 126
127 octave_idx_type max_maxord = 0; 127 octave_idx_type max_maxord = 0;
128 128
129 if (integration_method () == "stiff") 129 if (integration_method () == "stiff")
130 { 130 {
131 max_maxord = 5; 131 max_maxord = 5;
132 132
133 if (jac) 133 if (jac)
134 method_flag = 21; 134 method_flag = 21;
135 else 135 else
136 method_flag = 22; 136 method_flag = 22;
137 137
138 liw = 20 + n; 138 liw = 20 + n;
139 lrw = 22 + n * (9 + n); 139 lrw = 22 + n * (9 + n);
140 } 140 }
141 else 141 else
142 { 142 {
143 max_maxord = 12; 143 max_maxord = 12;
144 144
145 method_flag = 10; 145 method_flag = 10;
146 146
147 liw = 20; 147 liw = 20;
148 lrw = 22 + 16 * n; 148 lrw = 22 + 16 * n;
149 } 149 }
150 150
151 maxord = maximum_order (); 151 maxord = maximum_order ();
152 152
153 iwork.resize (liw); 153 iwork.resize (liw);
154 154
155 for (octave_idx_type i = 4; i < 9; i++) 155 for (octave_idx_type i = 4; i < 9; i++)
156 iwork(i) = 0; 156 iwork(i) = 0;
157 157
158 rwork.resize (lrw); 158 rwork.resize (lrw);
159 159
160 for (octave_idx_type i = 4; i < 9; i++) 160 for (octave_idx_type i = 4; i < 9; i++)
161 rwork(i) = 0; 161 rwork(i) = 0;
162 162
163 if (maxord >= 0) 163 if (maxord >= 0)
164 { 164 {
165 if (maxord > 0 && maxord <= max_maxord) 165 if (maxord > 0 && maxord <= max_maxord)
166 { 166 {
167 iwork(4) = maxord; 167 iwork(4) = maxord;
168 iopt = 1; 168 iopt = 1;
169 } 169 }
170 else 170 else
171 { 171 {
172 (*current_liboctave_error_handler) 172 (*current_liboctave_error_handler)
173 ("lsode: invalid value for maximum order"); 173 ("lsode: invalid value for maximum order");
174 integration_error = true; 174 integration_error = true;
175 return retval; 175 return retval;
176 } 176 }
177 } 177 }
178 178
179 if (stop_time_set) 179 if (stop_time_set)
180 { 180 {
181 itask = 4; 181 itask = 4;
182 rwork(0) = stop_time; 182 rwork(0) = stop_time;
183 iopt = 1; 183 iopt = 1;
184 } 184 }
185 else 185 else
186 { 186 {
187 itask = 1; 187 itask = 1;
188 } 188 }
189 189
190 px = x.fortran_vec (); 190 px = x.fortran_vec ();
191 191
192 piwork = iwork.fortran_vec (); 192 piwork = iwork.fortran_vec ();
193 prwork = rwork.fortran_vec (); 193 prwork = rwork.fortran_vec ();
206 user_jac = jacobian_function (); 206 user_jac = jacobian_function ();
207 207
208 ColumnVector xdot = (*user_fun) (x, t); 208 ColumnVector xdot = (*user_fun) (x, t);
209 209
210 if (x.length () != xdot.length ()) 210 if (x.length () != xdot.length ())
211 { 211 {
212 (*current_liboctave_error_handler) 212 (*current_liboctave_error_handler)
213 ("lsode: inconsistent sizes for state and derivative vectors"); 213 ("lsode: inconsistent sizes for state and derivative vectors");
214 214
215 integration_error = true; 215 integration_error = true;
216 return retval; 216 return retval;
217 } 217 }
218 218
219 ODEFunc::reset = false; 219 ODEFunc::reset = false;
220 220
221 // LSODE_options 221 // LSODE_options
222 222
224 abs_tol = absolute_tolerance (); 224 abs_tol = absolute_tolerance ();
225 225
226 octave_idx_type abs_tol_len = abs_tol.length (); 226 octave_idx_type abs_tol_len = abs_tol.length ();
227 227
228 if (abs_tol_len == 1) 228 if (abs_tol_len == 1)
229 itol = 1; 229 itol = 1;
230 else if (abs_tol_len == n) 230 else if (abs_tol_len == n)
231 itol = 2; 231 itol = 2;
232 else 232 else
233 { 233 {
234 (*current_liboctave_error_handler) 234 (*current_liboctave_error_handler)
235 ("lsode: inconsistent sizes for state and absolute tolerance vectors"); 235 ("lsode: inconsistent sizes for state and absolute tolerance vectors");
236 236
237 integration_error = true; 237 integration_error = true;
238 return retval; 238 return retval;
239 } 239 }
240 240
241 double iss = initial_step_size (); 241 double iss = initial_step_size ();
242 if (iss >= 0.0) 242 if (iss >= 0.0)
243 { 243 {
244 rwork(4) = iss; 244 rwork(4) = iss;
245 iopt = 1; 245 iopt = 1;
246 } 246 }
247 247
248 double maxss = maximum_step_size (); 248 double maxss = maximum_step_size ();
249 if (maxss >= 0.0) 249 if (maxss >= 0.0)
250 { 250 {
251 rwork(5) = maxss; 251 rwork(5) = maxss;
252 iopt = 1; 252 iopt = 1;
253 } 253 }
254 254
255 double minss = minimum_step_size (); 255 double minss = minimum_step_size ();
256 if (minss >= 0.0) 256 if (minss >= 0.0)
257 { 257 {
258 rwork(6) = minss; 258 rwork(6) = minss;
259 iopt = 1; 259 iopt = 1;
260 } 260 }
261 261
262 octave_idx_type sl = step_limit (); 262 octave_idx_type sl = step_limit ();
263 if (sl > 0) 263 if (sl > 0)
264 { 264 {
265 iwork(5) = sl; 265 iwork(5) = sl;
266 iopt = 1; 266 iopt = 1;
267 } 267 }
268 268
269 pabs_tol = abs_tol.fortran_vec (); 269 pabs_tol = abs_tol.fortran_vec ();
270 270
271 LSODE_options::reset = false; 271 LSODE_options::reset = false;
272 } 272 }
273 273
274 F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol, 274 F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol,
275 pabs_tol, itask, istate, iopt, prwork, lrw, 275 pabs_tol, itask, istate, iopt, prwork, lrw,
276 piwork, liw, lsode_j, method_flag)); 276 piwork, liw, lsode_j, method_flag));
277 277
278 switch (istate) 278 switch (istate)
279 { 279 {
280 case 1: // prior to initial integration step. 280 case 1: // prior to initial integration step.
281 case 2: // lsode was successful. 281 case 2: // lsode was successful.
286 case -1: // excess work done on this call (perhaps wrong mf). 286 case -1: // excess work done on this call (perhaps wrong mf).
287 case -2: // excess accuracy requested (tolerances too small). 287 case -2: // excess accuracy requested (tolerances too small).
288 case -3: // invalid input detected (see printed message). 288 case -3: // invalid input detected (see printed message).
289 case -4: // repeated error test failures (check all inputs). 289 case -4: // repeated error test failures (check all inputs).
290 case -5: // repeated convergence failures (perhaps bad jacobian 290 case -5: // repeated convergence failures (perhaps bad jacobian
291 // supplied or wrong choice of mf or tolerances). 291 // supplied or wrong choice of mf or tolerances).
292 case -6: // error weight became zero during problem. (solution 292 case -6: // error weight became zero during problem. (solution
293 // component i vanished, and atol or atol(i) = 0.) 293 // component i vanished, and atol or atol(i) = 0.)
294 case -13: // return requested in user-supplied function. 294 case -13: // return requested in user-supplied function.
295 integration_error = true; 295 integration_error = true;
296 break; 296 break;
297 297
298 default: 298 default:
299 integration_error = true; 299 integration_error = true;
300 (*current_liboctave_error_handler) 300 (*current_liboctave_error_handler)
301 ("unrecognized value of istate (= %d) returned from lsode", 301 ("unrecognized value of istate (= %d) returned from lsode",
302 istate); 302 istate);
303 break; 303 break;
304 } 304 }
305 305
306 return retval; 306 return retval;
307 } 307 }
322 break; 322 break;
323 323
324 case 2: 324 case 2:
325 retval = "successful exit"; 325 retval = "successful exit";
326 break; 326 break;
327 327
328 case 3: 328 case 3:
329 retval = "prior to continuation call with modified parameters"; 329 retval = "prior to continuation call with modified parameters";
330 break; 330 break;
331 331
332 case -1: 332 case -1:
333 retval = std::string ("excess work on this call (t = ") 333 retval = std::string ("excess work on this call (t = ")
334 + t_curr + "; perhaps wrong integration method)"; 334 + t_curr + "; perhaps wrong integration method)";
335 break; 335 break;
336 336
337 case -2: 337 case -2:
338 retval = "excess accuracy requested (tolerances too small)"; 338 retval = "excess accuracy requested (tolerances too small)";
339 break; 339 break;
342 retval = "invalid input detected (see printed message)"; 342 retval = "invalid input detected (see printed message)";
343 break; 343 break;
344 344
345 case -4: 345 case -4:
346 retval = std::string ("repeated error test failures (t = ") 346 retval = std::string ("repeated error test failures (t = ")
347 + t_curr + "check all inputs)"; 347 + t_curr + "check all inputs)";
348 break; 348 break;
349 349
350 case -5: 350 case -5:
351 retval = std::string ("repeated convergence failures (t = ") 351 retval = std::string ("repeated convergence failures (t = ")
352 + t_curr 352 + t_curr
353 + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)"; 353 + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)";
354 break; 354 break;
355 355
356 case -6: 356 case -6:
357 retval = std::string ("error weight became zero during problem. (t = ") 357 retval = std::string ("error weight became zero during problem. (t = ")
358 + t_curr 358 + t_curr
359 + "; solution component i vanished, and atol or atol(i) == 0)"; 359 + "; solution component i vanished, and atol or atol(i) == 0)";
360 break; 360 break;
361 361
362 case -13: 362 case -13:
363 retval = "return requested in user-supplied function (t = " 363 retval = "return requested in user-supplied function (t = "
364 + t_curr + ")"; 364 + t_curr + ")";
365 break; 365 break;
366 366
367 default: 367 default:
368 retval = "unknown error state"; 368 retval = "unknown error state";
369 break; 369 break;
383 if (n_out > 0 && n > 0) 383 if (n_out > 0 && n > 0)
384 { 384 {
385 retval.resize (n_out, n); 385 retval.resize (n_out, n);
386 386
387 for (octave_idx_type i = 0; i < n; i++) 387 for (octave_idx_type i = 0; i < n; i++)
388 retval.elem (0, i) = x.elem (i); 388 retval.elem (0, i) = x.elem (i);
389 389
390 for (octave_idx_type j = 1; j < n_out; j++) 390 for (octave_idx_type j = 1; j < n_out; j++)
391 { 391 {
392 ColumnVector x_next = do_integrate (tout.elem (j)); 392 ColumnVector x_next = do_integrate (tout.elem (j));
393 393
394 if (integration_error) 394 if (integration_error)
395 return retval; 395 return retval;
396 396
397 for (octave_idx_type i = 0; i < n; i++) 397 for (octave_idx_type i = 0; i < n; i++)
398 retval.elem (j, i) = x_next.elem (i); 398 retval.elem (j, i) = x_next.elem (i);
399 } 399 }
400 } 400 }
401 401
402 return retval; 402 return retval;
403 } 403 }
404 404
413 if (n_out > 0 && n > 0) 413 if (n_out > 0 && n > 0)
414 { 414 {
415 retval.resize (n_out, n); 415 retval.resize (n_out, n);
416 416
417 for (octave_idx_type i = 0; i < n; i++) 417 for (octave_idx_type i = 0; i < n; i++)
418 retval.elem (0, i) = x.elem (i); 418 retval.elem (0, i) = x.elem (i);
419 419
420 octave_idx_type n_crit = tcrit.capacity (); 420 octave_idx_type n_crit = tcrit.capacity ();
421 421
422 if (n_crit > 0) 422 if (n_crit > 0)
423 { 423 {
424 octave_idx_type i_crit = 0; 424 octave_idx_type i_crit = 0;
425 octave_idx_type i_out = 1; 425 octave_idx_type i_out = 1;
426 double next_crit = tcrit.elem (0); 426 double next_crit = tcrit.elem (0);
427 double next_out; 427 double next_out;
428 while (i_out < n_out) 428 while (i_out < n_out)
429 { 429 {
430 bool do_restart = false; 430 bool do_restart = false;
431 431
432 next_out = tout.elem (i_out); 432 next_out = tout.elem (i_out);
433 if (i_crit < n_crit) 433 if (i_crit < n_crit)
434 next_crit = tcrit.elem (i_crit); 434 next_crit = tcrit.elem (i_crit);
435 435
436 octave_idx_type save_output; 436 octave_idx_type save_output;
437 double t_out; 437 double t_out;
438 438
439 if (next_crit == next_out) 439 if (next_crit == next_out)
440 { 440 {
441 set_stop_time (next_crit); 441 set_stop_time (next_crit);
442 t_out = next_out; 442 t_out = next_out;
443 save_output = 1; 443 save_output = 1;
444 i_out++; 444 i_out++;
445 i_crit++; 445 i_crit++;
446 do_restart = true; 446 do_restart = true;
447 } 447 }
448 else if (next_crit < next_out) 448 else if (next_crit < next_out)
449 { 449 {
450 if (i_crit < n_crit) 450 if (i_crit < n_crit)
451 { 451 {
452 set_stop_time (next_crit); 452 set_stop_time (next_crit);
453 t_out = next_crit; 453 t_out = next_crit;
454 save_output = 0; 454 save_output = 0;
455 i_crit++; 455 i_crit++;
456 do_restart = true; 456 do_restart = true;
457 } 457 }
458 else 458 else
459 { 459 {
460 clear_stop_time (); 460 clear_stop_time ();
461 t_out = next_out; 461 t_out = next_out;
462 save_output = 1; 462 save_output = 1;
463 i_out++; 463 i_out++;
464 } 464 }
465 } 465 }
466 else 466 else
467 { 467 {
468 set_stop_time (next_crit); 468 set_stop_time (next_crit);
469 t_out = next_out; 469 t_out = next_out;
470 save_output = 1; 470 save_output = 1;
471 i_out++; 471 i_out++;
472 } 472 }
473 473
474 ColumnVector x_next = do_integrate (t_out); 474 ColumnVector x_next = do_integrate (t_out);
475 475
476 if (integration_error) 476 if (integration_error)
477 return retval; 477 return retval;
478 478
479 if (save_output) 479 if (save_output)
480 { 480 {
481 for (octave_idx_type i = 0; i < n; i++) 481 for (octave_idx_type i = 0; i < n; i++)
482 retval.elem (i_out-1, i) = x_next.elem (i); 482 retval.elem (i_out-1, i) = x_next.elem (i);
483 } 483 }
484 484
485 if (do_restart) 485 if (do_restart)
486 force_restart (); 486 force_restart ();
487 } 487 }
488 } 488 }
489 else 489 else
490 { 490 {
491 retval = do_integrate (tout); 491 retval = do_integrate (tout);
492 492
493 if (integration_error) 493 if (integration_error)
494 return retval; 494 return retval;
495 } 495 }
496 } 496 }
497 497
498 return retval; 498 return retval;
499 } 499 }