GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
DASSL.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1993-2018 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software: you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <https://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include <sstream>
28 
29 #include "DASSL.h"
30 #include "dMatrix.h"
31 #include "f77-fcn.h"
32 #include "lo-error.h"
33 #include "quit.h"
34 
35 typedef F77_INT (*dassl_fcn_ptr) (const double&, const double*,
36  const double*, double*, F77_INT&,
37  double*, F77_INT*);
38 
39 typedef F77_INT (*dassl_jac_ptr) (const double&, const double*,
40  const double*, double*, const double&,
41  double*, F77_INT*);
42 
43 extern "C"
44 {
45  F77_RET_T
46  F77_FUNC (ddassl, DDASSL) (dassl_fcn_ptr, const F77_INT&, F77_DBLE&,
47  F77_DBLE*, F77_DBLE*, F77_DBLE&, const F77_INT*,
48  const F77_DBLE*, const F77_DBLE*, F77_INT&,
49  F77_DBLE*, const F77_INT&, F77_INT*,
50  const F77_INT&, const F77_DBLE*, const F77_INT*,
52 }
53 
56 
57 static F77_INT nn;
58 
59 static F77_INT
60 ddassl_f (const double& time, const double *state, const double *deriv,
61  double *delta, F77_INT& ires, double *, F77_INT *)
62 {
63  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
64 
65  // FIXME: would be nice to avoid copying the data.
66 
67  ColumnVector tmp_deriv (nn);
68  ColumnVector tmp_state (nn);
69  ColumnVector tmp_delta (nn);
70 
71  for (F77_INT i = 0; i < nn; i++)
72  {
73  tmp_deriv.elem (i) = deriv[i];
74  tmp_state.elem (i) = state[i];
75  }
76 
77  octave_idx_type tmp_ires = ires;
78 
79  tmp_delta = user_fun (tmp_state, tmp_deriv, time, tmp_ires);
80 
81  ires = octave::to_f77_int (tmp_ires);
82 
83  if (ires >= 0)
84  {
85  if (tmp_delta.isempty ())
86  ires = -2;
87  else
88  {
89  for (F77_INT i = 0; i < nn; i++)
90  delta[i] = tmp_delta.elem (i);
91  }
92  }
93 
94  END_INTERRUPT_WITH_EXCEPTIONS;
95 
96  return 0;
97 }
98 
99 static F77_INT
100 ddassl_j (const double& time, const double *state, const double *deriv,
101  double *pd, const double& cj, double *, F77_INT *)
102 {
103  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
104 
105  // FIXME: would be nice to avoid copying the data.
106 
107  ColumnVector tmp_state (nn);
108  ColumnVector tmp_deriv (nn);
109 
110  for (F77_INT i = 0; i < nn; i++)
111  {
112  tmp_deriv.elem (i) = deriv[i];
113  tmp_state.elem (i) = state[i];
114  }
115 
116  Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
117 
118  for (F77_INT j = 0; j < nn; j++)
119  for (F77_INT i = 0; i < nn; i++)
120  pd[nn * j + i] = tmp_pd.elem (i, j);
121 
122  END_INTERRUPT_WITH_EXCEPTIONS;
123 
124  return 0;
125 }
126 
128 DASSL::do_integrate (double tout)
129 {
131 
132  if (! initialized || restart || DAEFunc::reset || DASSL_options::reset)
133  {
134  integration_error = false;
135 
136  initialized = true;
137 
138  info.resize (dim_vector (15, 1));
139 
140  for (F77_INT i = 0; i < 15; i++)
141  info(i) = 0;
142 
143  F77_INT n = octave::to_f77_int (size ());
144 
145  liw = 21 + n;
146  lrw = 40 + 9*n + n*n;
147 
148  nn = n;
149 
150  iwork.resize (dim_vector (liw, 1));
151  rwork.resize (dim_vector (lrw, 1));
152 
153  info(0) = 0;
154 
155  if (stop_time_set)
156  {
157  rwork(0) = stop_time;
158  info(3) = 1;
159  }
160  else
161  info(3) = 0;
162 
163  restart = false;
164 
165  // DAEFunc
166 
169 
170  if (user_fun)
171  {
172  octave_idx_type ires = 0;
173 
174  ColumnVector res = (*user_fun) (x, xdot, t, ires);
175 
176  if (res.numel () != x.numel ())
177  {
178  (*current_liboctave_error_handler)
179  ("dassl: inconsistent sizes for state and residual vectors");
180 
181  integration_error = true;
182  return retval;
183  }
184  }
185  else
186  {
187  (*current_liboctave_error_handler)
188  ("dassl: no user supplied RHS subroutine!");
189 
190  integration_error = true;
191  return retval;
192  }
193 
194  info(4) = (user_jac ? 1 : 0);
195 
196  DAEFunc::reset = false;
197 
198  // DASSL_options
199 
200  double hmax = maximum_step_size ();
201  if (hmax >= 0.0)
202  {
203  rwork(1) = hmax;
204  info(6) = 1;
205  }
206  else
207  info(6) = 0;
208 
209  double h0 = initial_step_size ();
210  if (h0 >= 0.0)
211  {
212  rwork(2) = h0;
213  info(7) = 1;
214  }
215  else
216  info(7) = 0;
217 
218  F77_INT sl = octave::to_f77_int (step_limit ());
219 
220  if (sl >= 0)
221  {
222  info(11) = 1;
223  iwork(20) = sl;
224  }
225  else
226  info(11) = 0;
227 
228  F77_INT maxord = octave::to_f77_int (maximum_order ());
229  if (maxord >= 0)
230  {
231  if (maxord > 0 && maxord < 6)
232  {
233  info(8) = 1;
234  iwork(2) = maxord;
235  }
236  else
237  {
238  (*current_liboctave_error_handler)
239  ("dassl: invalid value for maximum order: %d", maxord);
240  integration_error = true;
241  return retval;
242  }
243  }
244 
245  F77_INT enc = octave::to_f77_int (enforce_nonnegativity_constraints ());
246  info(9) = (enc ? 1 : 0);
247 
248  F77_INT ccic = octave::to_f77_int (compute_consistent_initial_condition ());
249  info(10) = (ccic ? 1 : 0);
250 
251  abs_tol = absolute_tolerance ();
252  rel_tol = relative_tolerance ();
253 
254  F77_INT abs_tol_len = octave::to_f77_int (abs_tol.numel ());
255  F77_INT rel_tol_len = octave::to_f77_int (rel_tol.numel ());
256 
257  if (abs_tol_len == 1 && rel_tol_len == 1)
258  {
259  info(1) = 0;
260  }
261  else if (abs_tol_len == n && rel_tol_len == n)
262  {
263  info(1) = 1;
264  }
265  else
266  {
267  (*current_liboctave_error_handler)
268  ("dassl: inconsistent sizes for tolerance arrays");
269 
270  integration_error = true;
271  return retval;
272  }
273 
274  DASSL_options::reset = false;
275  }
276 
277  double *px = x.fortran_vec ();
278  double *pxdot = xdot.fortran_vec ();
279 
280  F77_INT *pinfo = info.fortran_vec ();
281 
282  double *prel_tol = rel_tol.fortran_vec ();
283  double *pabs_tol = abs_tol.fortran_vec ();
284 
285  double *prwork = rwork.fortran_vec ();
286  F77_INT *piwork = iwork.fortran_vec ();
287 
288  double *dummy = nullptr;
289  F77_INT *idummy = nullptr;
290 
291  F77_INT tmp_istate = octave::to_f77_int (istate);
292 
293  F77_XFCN (ddassl, DDASSL, (ddassl_f, nn, t, px, pxdot, tout, pinfo,
294  prel_tol, pabs_tol, tmp_istate, prwork, lrw,
295  piwork, liw, dummy, idummy, ddassl_j));
296 
297  istate = tmp_istate;
298 
299  switch (istate)
300  {
301  case 1: // A step was successfully taken in intermediate-output
302  // mode. The code has not yet reached TOUT.
303  case 2: // The integration to TSTOP was successfully completed
304  // (T=TSTOP) by stepping exactly to TSTOP.
305  case 3: // The integration to TOUT was successfully completed
306  // (T=TOUT) by stepping past TOUT. Y(*) is obtained by
307  // interpolation. YPRIME(*) is obtained by interpolation.
308  retval = x;
309  t = tout;
310  break;
311 
312  case -1: // A large amount of work has been expended. (~500 steps).
313  case -2: // The error tolerances are too stringent.
314  case -3: // The local error test cannot be satisfied because you
315  // specified a zero component in ATOL and the
316  // corresponding computed solution component is zero.
317  // Thus, a pure relative error test is impossible for
318  // this component.
319  case -6: // DDASSL had repeated error test failures on the last
320  // attempted step.
321  case -7: // The corrector could not converge.
322  case -8: // The matrix of partial derivatives is singular.
323  case -9: // The corrector could not converge. There were repeated
324  // error test failures in this step.
325  case -10: // The corrector could not converge because IRES was
326  // equal to minus one.
327  case -11: // IRES equal to -2 was encountered and control is being
328  // returned to the calling program.
329  case -12: // DDASSL failed to compute the initial YPRIME.
330  case -33: // The code has encountered trouble from which it cannot
331  // recover. A message is printed explaining the trouble
332  // and control is returned to the calling program. For
333  // example, this occurs when invalid input is detected.
334  integration_error = true;
335  break;
336 
337  default:
338  integration_error = true;
339  (*current_liboctave_error_handler)
340  ("unrecognized value of istate (= %d) returned from ddassl",
341  istate);
342  break;
343  }
344 
345  return retval;
346 }
347 
348 Matrix
350 {
351  Matrix dummy;
352  return integrate (tout, dummy);
353 }
354 
355 Matrix
356 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
357 {
358  Matrix retval;
359 
360  octave_idx_type n_out = tout.numel ();
361  F77_INT n = octave::to_f77_int (size ());
362 
363  if (n_out > 0 && n > 0)
364  {
365  retval.resize (n_out, n);
366  xdot_out.resize (n_out, n);
367 
368  for (F77_INT i = 0; i < n; i++)
369  {
370  retval.elem (0, i) = x.elem (i);
371  xdot_out.elem (0, i) = xdot.elem (i);
372  }
373 
374  for (octave_idx_type j = 1; j < n_out; j++)
375  {
376  ColumnVector x_next = do_integrate (tout.elem (j));
377 
378  if (integration_error)
379  return retval;
380 
381  for (F77_INT i = 0; i < n; i++)
382  {
383  retval.elem (j, i) = x_next.elem (i);
384  xdot_out.elem (j, i) = xdot.elem (i);
385  }
386  }
387  }
388 
389  return retval;
390 }
391 
392 Matrix
393 DASSL::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
394 {
395  Matrix dummy;
396  return integrate (tout, dummy, tcrit);
397 }
398 
399 Matrix
400 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out,
401  const ColumnVector& tcrit)
402 {
403  Matrix retval;
404 
405  octave_idx_type n_out = tout.numel ();
406  F77_INT n = octave::to_f77_int (size ());
407 
408  if (n_out > 0 && n > 0)
409  {
410  retval.resize (n_out, n);
411  xdot_out.resize (n_out, n);
412 
413  for (F77_INT i = 0; i < n; i++)
414  {
415  retval.elem (0, i) = x.elem (i);
416  xdot_out.elem (0, i) = xdot.elem (i);
417  }
418 
419  octave_idx_type n_crit = tcrit.numel ();
420 
421  if (n_crit > 0)
422  {
423  octave_idx_type i_crit = 0;
424  octave_idx_type i_out = 1;
425  double next_crit = tcrit.elem (0);
426  double next_out;
427  while (i_out < n_out)
428  {
429  bool do_restart = false;
430 
431  next_out = tout.elem (i_out);
432  if (i_crit < n_crit)
433  next_crit = tcrit.elem (i_crit);
434 
435  bool save_output;
436  double t_out;
437 
438  if (next_crit == next_out)
439  {
440  set_stop_time (next_crit);
441  t_out = next_out;
442  save_output = true;
443  i_out++;
444  i_crit++;
445  do_restart = true;
446  }
447  else if (next_crit < next_out)
448  {
449  if (i_crit < n_crit)
450  {
451  set_stop_time (next_crit);
452  t_out = next_crit;
453  save_output = false;
454  i_crit++;
455  do_restart = true;
456  }
457  else
458  {
459  clear_stop_time ();
460  t_out = next_out;
461  save_output = true;
462  i_out++;
463  }
464  }
465  else
466  {
467  set_stop_time (next_crit);
468  t_out = next_out;
469  save_output = true;
470  i_out++;
471  }
472 
473  ColumnVector x_next = do_integrate (t_out);
474 
475  if (integration_error)
476  return retval;
477 
478  if (save_output)
479  {
480  for (F77_INT i = 0; i < n; i++)
481  {
482  retval.elem (i_out-1, i) = x_next.elem (i);
483  xdot_out.elem (i_out-1, i) = xdot.elem (i);
484  }
485  }
486 
487  if (do_restart)
488  force_restart ();
489  }
490  }
491  else
492  {
493  retval = integrate (tout, xdot_out);
494 
495  if (integration_error)
496  return retval;
497  }
498  }
499 
500  return retval;
501 }
502 
505 {
507 
508  std::ostringstream buf;
509  buf << t;
510  std::string t_curr = buf.str ();
511 
512  switch (istate)
513  {
514  case 1:
515  retval = "a step was successfully taken in intermediate-output mode.";
516  break;
517 
518  case 2:
519  retval = "integration completed by stepping exactly to TOUT";
520  break;
521 
522  case 3:
523  retval = "integration to tout completed by stepping past TOUT";
524  break;
525 
526  case -1:
527  retval = "a large amount of work has been expended (t =" + t_curr + ')';
528  break;
529 
530  case -2:
531  retval = "the error tolerances are too stringent";
532  break;
533 
534  case -3:
535  retval = "error weight became zero during problem. (t = " + t_curr +
536  "; solution component i vanished, and atol or atol(i) == 0)";
537  break;
538 
539  case -6:
540  retval = "repeated error test failures on the last attempted step (t = "
541  + t_curr + ')';
542  break;
543 
544  case -7:
545  retval = "the corrector could not converge (t = " + t_curr + ')';
546  break;
547 
548  case -8:
549  retval = "the matrix of partial derivatives is singular (t = " + t_curr +
550  ')';
551  break;
552 
553  case -9:
554  retval = "the corrector could not converge (t = " + t_curr +
555  "; repeated test failures)";
556  break;
557 
558  case -10:
559  retval = "corrector could not converge because IRES was -1 (t = "
560  + t_curr + ')';
561  break;
562 
563  case -11:
564  retval = "return requested in user-supplied function (t = " + t_curr +
565  ')';
566  break;
567 
568  case -12:
569  retval = "failed to compute consistent initial conditions";
570  break;
571 
572  case -33:
573  retval = "unrecoverable error (see printed message)";
574  break;
575 
576  default:
577  retval = "unknown error state";
578  break;
579  }
580 
581  return retval;
582 }
ColumnVector do_integrate(double t)
Definition: DASSL.cc:128
ColumnVector xdot
Definition: base-dae.h:77
static DAEFunc::DAEJacFunc user_jac
Definition: DASSL.cc:55
subroutine DDASSL(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL, IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
Definition: ddassl.f:3
bool initialized
Definition: DASSL.h:72
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:148
static DAEFunc::DAERHSFunc user_fun
Definition: DASSL.cc:54
Matrix integrate(const ColumnVector &tout, Matrix &xdot_out)
Definition: DASSL.cc:356
bool isempty(void) const
Definition: Array.h:565
Array< octave_f77_int_type > info
Definition: DASSL.h:77
const T * fortran_vec(void) const
Definition: Array.h:584
octave_f77_int_type liw
Definition: DASSL.h:74
std::string error_message(void) const
Definition: DASSL.cc:504
bool integration_error
Definition: base-de.h:115
F77_INT(* dassl_jac_ptr)(const double &, const double *, const double *, double *, const double &, double *, F77_INT *)
Definition: DASSL.cc:39
octave_idx_type istate
Definition: base-de.h:117
T & elem(octave_idx_type n)
Definition: Array.h:488
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:41
octave_f77_int_type lrw
Definition: DASSL.h:75
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)
octave_value resize(const dim_vector &dv, bool fill=false) const
Definition: ov.h:511
static F77_INT ddassl_f(const double &time, const double *state, const double *deriv, double *delta, F77_INT &ires, double *, F77_INT *)
Definition: DASSL.cc:60
bool reset
Definition: DAEFunc.h:101
octave_idx_type size(void) const
Definition: base-de.h:76
Matrix(* DAEJacFunc)(const ColumnVector &x, const ColumnVector &xdot, double t, double cj)
Definition: DAEFunc.h:44
virtual void force_restart(void)
Definition: base-de.h:95
Array< double > rel_tol
Definition: DASSL.h:83
void set_stop_time(double tt)
Definition: base-de.h:82
double t
Definition: base-de.h:107
DAEJacFunc jacobian_function(void) const
Definition: DAEFunc.h:82
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
Definition: Array.cc:1010
octave_value retval
Definition: data.cc:6246
Definition: dMatrix.h:36
void clear_stop_time(void)
Definition: base-de.h:89
static uint32_t state[624]
Definition: randmtzig.cc:183
Array< double > abs_tol
Definition: DASSL.h:82
double stop_time
Definition: base-de.h:109
ColumnVector x
Definition: base-de.h:105
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
DAERHSFunc function(void) const
Definition: DAEFunc.h:73
ColumnVector(* DAERHSFunc)(const ColumnVector &x, const ColumnVector &xdot, double t, octave_idx_type &ires)
Definition: DAEFunc.h:36
for i
Definition: data.cc:5264
Array< double > rwork
Definition: DASSL.h:80
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:366
bool restart
Definition: base-de.h:113
F77_INT(* dassl_fcn_ptr)(const double &, const double *, const double *, double *, F77_INT &, double *, F77_INT *)
Definition: DASSL.cc:35
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:888
F77_RET_T F77_FUNC(ddassl, DDASSL)(dassl_fcn_ptr
static F77_INT nn
Definition: DASSL.cc:57
bool stop_time_set
Definition: base-de.h:111
Array< octave_f77_int_type > iwork
Definition: DASSL.h:78
static F77_INT ddassl_j(const double &time, const double *state, const double *deriv, double *pd, const double &cj, double *, F77_INT *)
Definition: DASSL.cc:100