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