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
DASRT.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2002-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 "DASRT.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 (*dasrt_fcn_ptr) (const double&, const double*,
38  const double*, double*,
39  octave_idx_type&, double*,
41 
42 typedef octave_idx_type (*dasrt_jac_ptr) (const double&, const double*,
43  const double*, double*,
44  const double&, double*,
46 
48  const double&, const double*,
49  const octave_idx_type&,
50  double*, double*,
52 
53 extern "C"
54 {
55  F77_RET_T
57  double&, double*, double*, const double&,
58  octave_idx_type*, const double*,
59  const double*, octave_idx_type&, double*,
60  const octave_idx_type&, octave_idx_type*,
61  const octave_idx_type&, double*,
62  octave_idx_type*, dasrt_jac_ptr,
63  dasrt_constr_ptr, const octave_idx_type&,
64  octave_idx_type*);
65 }
66 
70 
72 
73 static octave_idx_type
74 ddasrt_f (const double& t, const double *state, const double *deriv,
75  double *delta, octave_idx_type& ires, double *, octave_idx_type *)
76 {
77  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
78 
79  ColumnVector tmp_state (nn);
80  ColumnVector tmp_deriv (nn);
81 
82  for (octave_idx_type i = 0; i < nn; i++)
83  {
84  tmp_state(i) = state[i];
85  tmp_deriv(i) = deriv[i];
86  }
87 
88  ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires);
89 
90  if (tmp_fval.length () == 0)
91  ires = -2;
92  else
93  {
94  for (octave_idx_type i = 0; i < nn; i++)
95  delta[i] = tmp_fval(i);
96  }
97 
98  END_INTERRUPT_WITH_EXCEPTIONS;
99 
100  return 0;
101 }
102 
104 ddasrt_j (const double& time, const double *state, const double *deriv,
105  double *pd, const double& cj, double *, octave_idx_type *)
106 {
107  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
108 
109  // FIXME: would be nice to avoid copying the data.
110 
111  ColumnVector tmp_state (nn);
112  ColumnVector tmp_deriv (nn);
113 
114  for (octave_idx_type i = 0; i < nn; i++)
115  {
116  tmp_deriv.elem (i) = deriv[i];
117  tmp_state.elem (i) = state[i];
118  }
119 
120  Matrix tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj);
121 
122  for (octave_idx_type j = 0; j < nn; j++)
123  for (octave_idx_type i = 0; i < nn; i++)
124  pd[nn * j + i] = tmp_pd.elem (i, j);
125 
126  END_INTERRUPT_WITH_EXCEPTIONS;
127 
128  return 0;
129 }
130 
131 static octave_idx_type
132 ddasrt_g (const octave_idx_type& neq, const double& t, const double *state,
133  const octave_idx_type& ng, double *gout, double *, octave_idx_type *)
134 {
135  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
136 
137  octave_idx_type n = neq;
138 
139  ColumnVector tmp_state (n);
140  for (octave_idx_type i = 0; i < n; i++)
141  tmp_state(i) = state[i];
142 
143  ColumnVector tmp_fval = (*user_csub) (tmp_state, t);
144 
145  for (octave_idx_type i = 0; i < ng; i++)
146  gout[i] = tmp_fval(i);
147 
148  END_INTERRUPT_WITH_EXCEPTIONS;
149 
150  return 0;
151 }
152 
153 void
154 DASRT::integrate (double tout)
155 {
156  DASRT_result retval;
157 
158  // I suppose this is the safe thing to do. If this is the first
159  // call, or if anything about the problem has changed, we should
160  // start completely fresh.
161 
162  if (! initialized || restart
164  {
165  integration_error = false;
166 
167  initialized = true;
168 
169  info.resize (dim_vector (15, 1));
170 
171  for (octave_idx_type i = 0; i < 15; i++)
172  info(i) = 0;
173 
174  octave_idx_type n = size ();
175 
176  nn = n;
177 
178  // DAERTFunc
179 
181 
182  if (user_csub)
183  {
184  ColumnVector tmp = (*user_csub) (x, t);
185  ng = tmp.length ();
186  }
187  else
188  ng = 0;
189 
190  octave_idx_type maxord = maximum_order ();
191  if (maxord >= 0)
192  {
193  if (maxord > 0 && maxord < 6)
194  {
195  info(8) = 1;
196  iwork(2) = maxord;
197  }
198  else
199  {
200  (*current_liboctave_error_handler)
201  ("dassl: invalid value for maximum order");
202  integration_error = true;
203  return;
204  }
205  }
206 
207  liw = 21 + n;
208  lrw = 50 + 9*n + n*n + 3*ng;
209 
210  iwork.resize (dim_vector (liw, 1));
211  rwork.resize (dim_vector (lrw, 1));
212 
213  info(0) = 0;
214 
215  if (stop_time_set)
216  {
217  info(3) = 1;
218  rwork(0) = stop_time;
219  }
220  else
221  info(3) = 0;
222 
223  restart = false;
224 
225  // DAEFunc
226 
229 
230  if (user_fsub)
231  {
232  octave_idx_type ires = 0;
233 
234  ColumnVector fval = (*user_fsub) (x, xdot, t, ires);
235 
236  if (fval.length () != x.length ())
237  {
238  (*current_liboctave_error_handler)
239  ("dasrt: inconsistent sizes for state and residual vectors");
240 
241  integration_error = true;
242  return;
243  }
244  }
245  else
246  {
247  (*current_liboctave_error_handler)
248  ("dasrt: no user supplied RHS subroutine!");
249 
250  integration_error = true;
251  return;
252  }
253 
254  info(4) = user_jsub ? 1 : 0;
255 
256  DAEFunc::reset = false;
257 
258  jroot.resize (dim_vector (ng, 1), 1);
259 
260  DAERTFunc::reset = false;
261 
262  // DASRT_options
263 
264  double mss = maximum_step_size ();
265  if (mss >= 0.0)
266  {
267  rwork(1) = mss;
268  info(6) = 1;
269  }
270  else
271  info(6) = 0;
272 
273  double iss = initial_step_size ();
274  if (iss >= 0.0)
275  {
276  rwork(2) = iss;
277  info(7) = 1;
278  }
279  else
280  info(7) = 0;
281 
282  if (step_limit () >= 0)
283  {
284  info(11) = 1;
285  iwork(20) = step_limit ();
286  }
287  else
288  info(11) = 0;
289 
292 
293  octave_idx_type abs_tol_len = abs_tol.length ();
294  octave_idx_type rel_tol_len = rel_tol.length ();
295 
296  if (abs_tol_len == 1 && rel_tol_len == 1)
297  {
298  info.elem (1) = 0;
299  }
300  else if (abs_tol_len == n && rel_tol_len == n)
301  {
302  info.elem (1) = 1;
303  }
304  else
305  {
306  (*current_liboctave_error_handler)
307  ("dasrt: inconsistent sizes for tolerance arrays");
308 
309  integration_error = true;
310  return;
311  }
312 
313  DASRT_options::reset = false;
314  }
315 
316  double *px = x.fortran_vec ();
317  double *pxdot = xdot.fortran_vec ();
318 
319  octave_idx_type *pinfo = info.fortran_vec ();
320 
321  double *prel_tol = rel_tol.fortran_vec ();
322  double *pabs_tol = abs_tol.fortran_vec ();
323 
324  double *prwork = rwork.fortran_vec ();
325  octave_idx_type *piwork = iwork.fortran_vec ();
326 
327  octave_idx_type *pjroot = jroot.fortran_vec ();
328 
329  double *dummy = 0;
330  octave_idx_type *idummy = 0;
331 
332  F77_XFCN (ddasrt, DDASRT, (ddasrt_f, nn, t, px, pxdot, tout, pinfo,
333  prel_tol, pabs_tol, istate, prwork, lrw,
334  piwork, liw, dummy, idummy, ddasrt_j,
335  ddasrt_g, ng, pjroot));
336 
337  switch (istate)
338  {
339  case 1: // A step was successfully taken in intermediate-output
340  // mode. The code has not yet reached TOUT.
341  case 2: // The integration to TOUT was successfully completed
342  // (T=TOUT) by stepping exactly to TOUT.
343  case 3: // The integration to TOUT was successfully completed
344  // (T=TOUT) by stepping past TOUT. Y(*) is obtained by
345  // interpolation. YPRIME(*) is obtained by interpolation.
346  t = tout;
347  break;
348 
349  case 4: // The integration was successfully completed
350  // by finding one or more roots of G at T.
351  break;
352 
353  case -1: // A large amount of work has been expended.
354  case -2: // The error tolerances are too stringent.
355  case -3: // The local error test cannot be satisfied because you
356  // specified a zero component in ATOL and the
357  // corresponding computed solution component is zero.
358  // Thus, a pure relative error test is impossible for
359  // this component.
360  case -6: // DDASRT had repeated error test failures on the last
361  // attempted step.
362  case -7: // The corrector could not converge.
363  case -8: // The matrix of partial derivatives is singular.
364  case -9: // The corrector could not converge. There were repeated
365  // error test failures in this step.
366  case -10: // The corrector could not converge because IRES was
367  // equal to minus one.
368  case -11: // IRES equal to -2 was encountered and control is being
369  // returned to the calling program.
370  case -12: // DASSL failed to compute the initial YPRIME.
371  case -33: // The code has encountered trouble from which it cannot
372  // recover. A message is printed explaining the trouble
373  // and control is returned to the calling program. For
374  // example, this occurs when invalid input is detected.
375  integration_error = true;
376  break;
377 
378  default:
379  integration_error = true;
380  (*current_liboctave_error_handler)
381  ("unrecognized value of istate (= %d) returned from ddasrt",
382  istate);
383  break;
384  }
385 }
386 
389 {
390  DASRT_result retval;
391 
392  Matrix x_out;
393  Matrix xdot_out;
394  ColumnVector t_out = tout;
395 
396  octave_idx_type n_out = tout.capacity ();
397  octave_idx_type n = size ();
398 
399  if (n_out > 0 && n > 0)
400  {
401  x_out.resize (n_out, n);
402  xdot_out.resize (n_out, n);
403 
404  for (octave_idx_type i = 0; i < n; i++)
405  {
406  x_out(0,i) = x(i);
407  xdot_out(0,i) = xdot(i);
408  }
409 
410  for (octave_idx_type j = 1; j < n_out; j++)
411  {
412  integrate (tout(j));
413 
414  if (integration_error)
415  {
416  retval = DASRT_result (x_out, xdot_out, t_out);
417  return retval;
418  }
419 
420  if (istate == 4)
421  t_out(j) = t;
422  else
423  t_out(j) = tout(j);
424 
425  for (octave_idx_type i = 0; i < n; i++)
426  {
427  x_out(j,i) = x(i);
428  xdot_out(j,i) = xdot(i);
429  }
430 
431  if (istate == 4)
432  {
433  x_out.resize (j+1, n);
434  xdot_out.resize (j+1, n);
435  t_out.resize (j+1);
436  break;
437  }
438  }
439  }
440 
441  retval = DASRT_result (x_out, xdot_out, t_out);
442 
443  return retval;
444 }
445 
447 DASRT::integrate (const ColumnVector& tout, const ColumnVector& tcrit)
448 {
449  DASRT_result retval;
450 
451  Matrix x_out;
452  Matrix xdot_out;
453  ColumnVector t_outs = tout;
454 
455  octave_idx_type n_out = tout.capacity ();
456  octave_idx_type n = size ();
457 
458  if (n_out > 0 && n > 0)
459  {
460  x_out.resize (n_out, n);
461  xdot_out.resize (n_out, n);
462 
463  octave_idx_type n_crit = tcrit.capacity ();
464 
465  if (n_crit > 0)
466  {
467  octave_idx_type i_crit = 0;
468  octave_idx_type i_out = 1;
469  double next_crit = tcrit(0);
470  double next_out;
471  while (i_out < n_out)
472  {
473  bool do_restart = false;
474 
475  next_out = tout(i_out);
476  if (i_crit < n_crit)
477  next_crit = tcrit(i_crit);
478 
479  octave_idx_type save_output;
480  double t_out;
481 
482  if (next_crit == next_out)
483  {
484  set_stop_time (next_crit);
485  t_out = next_out;
486  save_output = 1;
487  i_out++;
488  i_crit++;
489  do_restart = true;
490  }
491  else if (next_crit < next_out)
492  {
493  if (i_crit < n_crit)
494  {
495  set_stop_time (next_crit);
496  t_out = next_crit;
497  save_output = 0;
498  i_crit++;
499  do_restart = true;
500  }
501  else
502  {
503  clear_stop_time ();
504  t_out = next_out;
505  save_output = 1;
506  i_out++;
507  }
508  }
509  else
510  {
511  set_stop_time (next_crit);
512  t_out = next_out;
513  save_output = 1;
514  i_out++;
515  }
516 
517  integrate (t_out);
518 
519  if (integration_error)
520  {
521  retval = DASRT_result (x_out, xdot_out, t_outs);
522  return retval;
523  }
524 
525  if (istate == 4)
526  t_out = t;
527 
528  if (save_output)
529  {
530  for (octave_idx_type i = 0; i < n; i++)
531  {
532  x_out(i_out-1,i) = x(i);
533  xdot_out(i_out-1,i) = xdot(i);
534  }
535 
536  t_outs(i_out-1) = t_out;
537 
538  if (istate == 4)
539  {
540  x_out.resize (i_out, n);
541  xdot_out.resize (i_out, n);
542  t_outs.resize (i_out);
543  i_out = n_out;
544  }
545  }
546 
547  if (do_restart)
548  force_restart ();
549  }
550 
551  retval = DASRT_result (x_out, xdot_out, t_outs);
552  }
553  else
554  {
555  retval = integrate (tout);
556 
557  if (integration_error)
558  return retval;
559  }
560  }
561 
562  return retval;
563 }
564 
565 std::string
567 {
568  std::string retval;
569 
570  std::ostringstream buf;
571  buf << t;
572  std::string t_curr = buf.str ();
573 
574  switch (istate)
575  {
576  case 1:
577  retval = "a step was successfully taken in intermediate-output mode.";
578  break;
579 
580  case 2:
581  retval = "integration completed by stepping exactly to TOUT";
582  break;
583 
584  case 3:
585  retval = "integration to tout completed by stepping past TOUT";
586  break;
587 
588  case 4:
589  retval = "integration completed by finding one or more roots of G at T";
590  break;
591 
592  case -1:
593  retval = std::string ("a large amount of work has been expended (t =")
594  + t_curr + ")";
595  break;
596 
597  case -2:
598  retval = "the error tolerances are too stringent";
599  break;
600 
601  case -3:
602  retval = std::string ("error weight became zero during problem. (t = ")
603  + t_curr
604  + "; solution component i vanished, and atol or atol(i) == 0)";
605  break;
606 
607  case -6:
608  retval = std::string ("repeated error test failures on the last attempted step (t = ")
609  + t_curr + ")";
610  break;
611 
612  case -7:
613  retval = std::string ("the corrector could not converge (t = ")
614  + t_curr + ")";
615  break;
616 
617  case -8:
618  retval = std::string ("the matrix of partial derivatives is singular (t = ")
619  + t_curr + ")";
620  break;
621 
622  case -9:
623  retval = std::string ("the corrector could not converge (t = ")
624  + t_curr + "; repeated test failures)";
625  break;
626 
627  case -10:
628  retval = std::string ("corrector could not converge because IRES was -1 (t = ")
629  + t_curr + ")";
630  break;
631 
632  case -11:
633  retval = std::string ("return requested in user-supplied function (t = ")
634  + t_curr + ")";
635  break;
636 
637  case -12:
638  retval = "failed to compute consistent initial conditions";
639  break;
640 
641  case -33:
642  retval = "unrecoverable error (see printed message)";
643  break;
644 
645  default:
646  retval = "unknown error state";
647  break;
648  }
649 
650  return retval;
651 }