GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
LSODE.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1993-2023 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
49  F77_DBLE&, F77_INT&, F77_DBLE&, const F77_DBLE *,
52  F77_INT&);
53 }
54 
59 
60 static F77_INT
61 lsode_f (const F77_INT& neq, const double& time, double *, double *deriv,
62  F77_INT& ierr)
63 {
64  ColumnVector tmp_deriv;
65 
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 
70  tmp_deriv = (*user_fcn) (*tmp_x, time);
71 
72  if (tmp_deriv.isempty ())
73  ierr = -1;
74  else
75  {
76  for (F77_INT i = 0; i < neq; i++)
77  deriv[i] = tmp_deriv.elem (i);
78  }
79 
80  return 0;
81 }
82 
83 static F77_INT
84 lsode_j (const F77_INT& neq, const double& time, double *,
85  const F77_INT& ml, const F77_INT& mu,
86  double *pd, const F77_INT& nrowpd)
87 {
88  Matrix tmp_jac (neq, neq);
89 
90  // NOTE: this won't work if LSODE passes copies of the state vector.
91  // In that case we have to create a temporary vector object
92  // and copy.
93 
94  tmp_jac = (*user_jac) (*tmp_x, time);
95 
97  for (F77_INT j = 0; j < neq; j++)
98  for (F77_INT i = 0; i < neq; i++)
99  pd[nrowpd * j + i] = tmp_jac (i, j);
100  else
101  // upper left ends of subdiagonals in tmp_jac
102  for (F77_INT i = 0, j = mu; i <= ml; j == 0 ? i++ : j--)
103  // walk down the subdiagonal in tmp_jac
104  for (F77_INT k = i, l = j; k < neq && l < neq; k++, l++)
105  pd[nrowpd * l + k + mu - l] = tmp_jac (k, l);
106 
107  return 0;
108 }
109 
111 LSODE::do_integrate (double tout)
112 {
113  ColumnVector retval;
114 
115  static F77_INT nn = 0;
116 
118  || LSODE_options::m_reset)
119  {
120  m_integration_error = false;
121 
122  m_initialized = true;
123 
124  m_istate = 1;
125 
126  F77_INT n = octave::to_f77_int (size ());
127 
128  nn = n;
129 
130  octave_idx_type max_maxord = 0;
131 
132  user_jac_ignore_ml_mu = true;
133 
135 
136  m_iwork(0) = lower_jacobian_subdiagonals (); // 'ML' in dlsode.f
137 
138  m_iwork(1) = upper_jacobian_subdiagonals (); // 'MU' in dlsode.f
139 
140  if (integration_method () == "stiff")
141  {
142  max_maxord = 5;
143 
144  if (m_jac)
145  {
146  if (jacobian_type () == "banded")
147  {
148  m_method_flag = 24;
149  user_jac_ignore_ml_mu = false;
150  }
151  else
152  m_method_flag = 21;
153  }
154  else
155  {
156  if (jacobian_type () == "full")
157  m_method_flag = 22;
158  else if (jacobian_type () == "banded")
159  m_method_flag = 25;
160  else if (jacobian_type () == "diagonal")
161  m_method_flag = 23;
162  else
163  {
164  // should be prevented by lsode_options
165  (*current_liboctave_error_handler)
166  ("lsode: internal error, wrong jacobian type");
167  m_integration_error = true;
168  return retval;
169  }
170  }
171 
172  m_liw = 20 + n;
173  m_lrw = 22 + n * (9 + n);
174  }
175  else
176  {
177  max_maxord = 12;
178 
179  m_method_flag = 10;
180 
181  m_liw = 20;
182  m_lrw = 22 + 16 * n;
183  }
184 
186 
187  for (F77_INT i = 4; i < 9; i++)
188  m_iwork(i) = 0;
189 
191 
192  for (F77_INT i = 4; i < 9; i++)
193  m_rwork(i) = 0;
194 
195  octave_idx_type maxord = maximum_order ();
196 
197  if (maxord >= 0)
198  {
199  if (maxord > 0 && maxord <= max_maxord)
200  {
201  m_iwork(4) = octave::to_f77_int (maxord);
202  m_iopt = 1;
203  }
204  else
205  {
206  // FIXME: Should this be a warning?
207  (*current_liboctave_error_handler)
208  ("lsode: invalid value for maximum order");
209  m_integration_error = true;
210  return retval;
211  }
212  }
213 
214  if (m_stop_time_set)
215  {
216  m_itask = 4;
217  m_rwork(0) = m_stop_time;
218  m_iopt = 1;
219  }
220  else
221  {
222  m_itask = 1;
223  }
224 
225  m_restart = false;
226 
227  // ODEFunc
228 
229  // NOTE: this won't work if LSODE passes copies of the state vector.
230  // In that case we have to create a temporary vector object
231  // and copy.
232 
233  tmp_x = &m_x;
234 
235  user_fcn = function ();
237 
238  ColumnVector m_xdot = (*user_fcn) (m_x, m_t);
239 
240  if (m_x.numel () != m_xdot.numel ())
241  {
242  // FIXME: Should this be a warning?
243  (*current_liboctave_error_handler)
244  ("lsode: inconsistent sizes for state and derivative vectors");
245 
246  m_integration_error = true;
247  return retval;
248  }
249 
250  ODEFunc::m_reset = false;
251 
252  // LSODE_options
253 
254  m_rel_tol = relative_tolerance ();
255  m_abs_tol = absolute_tolerance ();
256 
257  F77_INT abs_tol_len = octave::to_f77_int (m_abs_tol.numel ());
258 
259  if (abs_tol_len == 1)
260  m_itol = 1;
261  else if (abs_tol_len == n)
262  m_itol = 2;
263  else
264  {
265  // FIXME: Should this be a warning?
266  (*current_liboctave_error_handler)
267  ("lsode: inconsistent sizes for state and absolute tolerance vectors");
268 
269  m_integration_error = true;
270  return retval;
271  }
272 
273  double iss = initial_step_size ();
274  if (iss >= 0.0)
275  {
276  m_rwork(4) = iss;
277  m_iopt = 1;
278  }
279 
280  double maxss = maximum_step_size ();
281  if (maxss >= 0.0)
282  {
283  m_rwork(5) = maxss;
284  m_iopt = 1;
285  }
286 
287  double minss = minimum_step_size ();
288  if (minss >= 0.0)
289  {
290  m_rwork(6) = minss;
291  m_iopt = 1;
292  }
293 
294  F77_INT sl = octave::to_f77_int (step_limit ());
295  if (sl > 0)
296  {
297  m_iwork(5) = sl;
298  m_iopt = 1;
299  }
300 
301  LSODE_options::m_reset = false;
302  }
303 
304  double *px = m_x.fortran_vec ();
305 
306  double *pabs_tol = m_abs_tol.fortran_vec ();
307 
308  F77_INT *piwork = m_iwork.fortran_vec ();
309  double *prwork = m_rwork.fortran_vec ();
310 
311  F77_INT tmp_istate = octave::to_f77_int (m_istate);
312 
313  F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, m_t, tout, m_itol, m_rel_tol,
314  pabs_tol, m_itask, tmp_istate, m_iopt, prwork,
315  m_lrw, piwork, m_liw, lsode_j, m_method_flag));
316 
317  m_istate = tmp_istate;
318 
319  switch (m_istate)
320  {
321  case 1: // prior to initial integration step.
322  case 2: // lsode was successful.
323  retval = m_x;
324  m_t = tout;
325  break;
326 
327  case -1: // excess work done on this call (perhaps wrong mf).
328  case -2: // excess accuracy requested (tolerances too small).
329  case -3: // invalid input detected (see printed message).
330  case -4: // repeated error test failures (check all inputs).
331  case -5: // repeated convergence failures (perhaps bad Jacobian
332  // supplied or wrong choice of mf or tolerances).
333  case -6: // error weight became zero during problem. (solution
334  // component i vanished, and atol or atol(i) = 0.)
335  case -13: // return requested in user-supplied function.
336  m_integration_error = true;
337  break;
338 
339  default:
340  m_integration_error = true;
341  (*current_liboctave_error_handler)
342  ("unrecognized value of istate (= %" OCTAVE_IDX_TYPE_FORMAT ") "
343  "returned from lsode", m_istate);
344  break;
345  }
346 
347  return retval;
348 }
349 
350 std::string
352 {
353  std::string retval;
354 
355  std::ostringstream buf;
356  buf << m_t;
357  std::string t_curr = buf.str ();
358 
359  switch (m_istate)
360  {
361  case 1:
362  retval = "prior to initial integration step";
363  break;
364 
365  case 2:
366  retval = "successful exit";
367  break;
368 
369  case 3:
370  retval = "prior to continuation call with modified parameters";
371  break;
372 
373  case -1:
374  retval = "excess work on this call (t = " + t_curr +
375  "; perhaps wrong integration method)";
376  break;
377 
378  case -2:
379  retval = "excess accuracy requested (tolerances too small)";
380  break;
381 
382  case -3:
383  retval = "invalid input detected (see printed message)";
384  break;
385 
386  case -4:
387  retval = "repeated error test failures (t = " + t_curr +
388  "; check all inputs)";
389  break;
390 
391  case -5:
392  retval = "repeated convergence failures (t = " + t_curr +
393  "; perhaps bad Jacobian supplied or wrong choice of integration method or tolerances)";
394  break;
395 
396  case -6:
397  retval = "error weight became zero during problem. (t = " + t_curr +
398  "; solution component i vanished, and atol or atol(i) == 0)";
399  break;
400 
401  case -13:
402  retval = "return requested in user-supplied function (t = "
403  + t_curr + ')';
404  break;
405 
406  default:
407  retval = "unknown error state";
408  break;
409  }
410 
411  return retval;
412 }
413 
414 Matrix
416 {
417  Matrix retval;
418 
419  octave_idx_type n_out = tout.numel ();
420  F77_INT n = octave::to_f77_int (size ());
421 
422  if (n_out > 0 && n > 0)
423  {
424  retval.resize (n_out, n);
425 
426  for (F77_INT i = 0; i < n; i++)
427  retval.elem (0, i) = m_x.elem (i);
428 
429  for (octave_idx_type j = 1; j < n_out; j++)
430  {
431  ColumnVector x_next = do_integrate (tout.elem (j));
432 
434  return retval;
435 
436  for (F77_INT i = 0; i < n; i++)
437  retval.elem (j, i) = x_next.elem (i);
438  }
439  }
440 
441  return retval;
442 }
443 
444 Matrix
445 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
446 {
447  Matrix retval;
448 
449  octave_idx_type n_out = tout.numel ();
450  F77_INT n = octave::to_f77_int (size ());
451 
452  if (n_out > 0 && n > 0)
453  {
454  retval.resize (n_out, n);
455 
456  for (F77_INT i = 0; i < n; i++)
457  retval.elem (0, i) = m_x.elem (i);
458 
459  octave_idx_type n_crit = tcrit.numel ();
460 
461  if (n_crit > 0)
462  {
463  octave_idx_type i_crit = 0;
464  octave_idx_type i_out = 1;
465  double next_crit = tcrit.elem (0);
466  double next_out;
467  while (i_out < n_out)
468  {
469  bool do_restart = false;
470 
471  next_out = tout.elem (i_out);
472  if (i_crit < n_crit)
473  next_crit = tcrit.elem (i_crit);
474 
475  bool save_output = false;
476  double t_out;
477 
478  if (next_crit == next_out)
479  {
480  set_stop_time (next_crit);
481  t_out = next_out;
482  save_output = true;
483  i_out++;
484  i_crit++;
485  do_restart = true;
486  }
487  else if (next_crit < next_out)
488  {
489  if (i_crit < n_crit)
490  {
491  set_stop_time (next_crit);
492  t_out = next_crit;
493  save_output = false;
494  i_crit++;
495  do_restart = true;
496  }
497  else
498  {
499  clear_stop_time ();
500  t_out = next_out;
501  save_output = true;
502  i_out++;
503  }
504  }
505  else
506  {
507  set_stop_time (next_crit);
508  t_out = next_out;
509  save_output = true;
510  i_out++;
511  }
512 
513  ColumnVector x_next = do_integrate (t_out);
514 
516  return retval;
517 
518  if (save_output)
519  {
520  for (F77_INT i = 0; i < n; i++)
521  retval.elem (i_out-1, i) = x_next.elem (i);
522  }
523 
524  if (do_restart)
525  force_restart ();
526  }
527  }
528  else
529  {
530  retval = do_integrate (tout);
531 
533  return retval;
534  }
535  }
536 
537  return retval;
538 }
static F77_INT nn
Definition: DASPK.cc:66
F77_INT(* lsode_jac_ptr)(const F77_INT &, const double &, double *, const F77_INT &, const F77_INT &, double *, const F77_INT &)
Definition: LSODE.cc:41
static F77_INT lsode_f(const F77_INT &neq, const double &time, double *, double *deriv, F77_INT &ierr)
Definition: LSODE.cc:61
static bool user_jac_ignore_ml_mu
Definition: LSODE.cc:58
F77_INT(* lsode_fcn_ptr)(const F77_INT &, const double &, double *, double *, F77_INT &)
Definition: LSODE.cc:38
F77_RET_T F77_INT F77_DBLE F77_DBLE F77_DBLE F77_INT F77_DBLE const F77_DBLE F77_INT F77_INT F77_INT F77_DBLE F77_INT F77_INT F77_INT F77_INT &static ODEFunc::ODERHSFunc user_fcn
Definition: LSODE.cc:52
static ODEFunc::ODEJacFunc user_jac
Definition: LSODE.cc:56
F77_RET_T F77_FUNC(dlsode, DLSODE)(lsode_fcn_ptr
static F77_INT lsode_j(const F77_INT &neq, const double &time, double *, const F77_INT &ml, const F77_INT &mu, double *pd, const F77_INT &nrowpd)
Definition: LSODE.cc:84
static ColumnVector * tmp_x
Definition: LSODE.cc:57
OCTARRAY_OVERRIDABLE_FUNC_API bool isempty(void) const
Size of the specified dimension.
Definition: Array.h:651
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:414
OCTARRAY_API void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array-base.cc:1032
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array-base.cc:1766
OCTARRAY_OVERRIDABLE_FUNC_API T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
Array< double > m_rwork
Definition: LSODE.h:75
Array< double > m_abs_tol
Definition: LSODE.h:79
octave_f77_int_type m_lrw
Definition: LSODE.h:72
octave_f77_int_type m_liw
Definition: LSODE.h:71
double m_rel_tol
Definition: LSODE.h:77
std::string error_message(void) const
Definition: LSODE.cc:351
ColumnVector do_integrate(double t)
Definition: LSODE.cc:111
octave_f77_int_type m_itol
Definition: LSODE.h:69
octave_f77_int_type m_method_flag
Definition: LSODE.h:66
bool m_initialized
Definition: LSODE.h:64
octave_f77_int_type m_iopt
Definition: LSODE.h:68
octave_f77_int_type m_itask
Definition: LSODE.h:67
Array< octave_f77_int_type > m_iwork
Definition: LSODE.h:74
Definition: dMatrix.h:42
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
ODEJacFunc m_jac
Definition: ODEFunc.h:87
ColumnVector(* ODERHSFunc)(const ColumnVector &, double)
Definition: ODEFunc.h:38
Matrix(* ODEJacFunc)(const ColumnVector &, double)
Definition: ODEFunc.h:39
bool m_reset
Definition: ODEFunc.h:94
ODEJacFunc jacobian_function(void) const
Definition: ODEFunc.h:75
virtual void force_restart(void)
Definition: base-de.h:98
double m_stop_time
Definition: base-de.h:112
bool m_restart
Definition: base-de.h:116
double m_t
Definition: base-de.h:110
octave_idx_type size(void) const
Definition: base-de.h:79
octave_idx_type m_istate
Definition: base-de.h:120
bool m_stop_time_set
Definition: base-de.h:114
bool m_integration_error
Definition: base-de.h:118
ColumnVector m_x
Definition: base-de.h:108
void set_stop_time(double tt)
Definition: base-de.h:85
void clear_stop_time(void)
Definition: base-de.h:92
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
Definition: oct-time.h:52
subroutine dlsode(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
Definition: dlsode.f:3
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
double F77_DBLE
Definition: f77-fcn.h:302
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
F77_RET_T(F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, const F77_INT &, const F77_INT &, const F77_INT &, F77_INT &, F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_DBLE *, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, const F77_INT &, F77_DBLE *, F77_INT *, F77_INT &F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL)
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE const F77_INT F77_INT & ierr
octave_idx_type n
Definition: mx-inlines.cc:753