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