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
DASPK.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-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 "DASPK.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 (*daspk_fcn_ptr) (const double&, const double*,
38  const double*, const double&,
39  double*, octave_idx_type&,
40  double*, octave_idx_type*);
41 
42 typedef octave_idx_type (*daspk_jac_ptr) (const double&, const double*,
43  const double*, double*,
44  const double&, double*,
46 
48  const double&, const double*,
49  const double*, const double*,
50  const double&, const double*,
51  double*, octave_idx_type*,
52  double*, const double&,
53  octave_idx_type&, double*,
55 
56 extern "C"
57 {
58  F77_RET_T
59  F77_FUNC (ddaspk, DDASPK) (daspk_fcn_ptr, const octave_idx_type&,
60  double&, double*, double*, double&,
61  const octave_idx_type*, const double*,
62  const double*, octave_idx_type&,
63  double*, const octave_idx_type&,
64  octave_idx_type*, const octave_idx_type&,
65  const double*, const octave_idx_type*,
67 }
68 
72 
73 static octave_idx_type
74 ddaspk_f (const double& time, const double *state, const double *deriv,
75  const double&, double *delta, octave_idx_type& ires, double *,
77 {
78  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
79 
80  ColumnVector tmp_deriv (nn);
81  ColumnVector tmp_state (nn);
82  ColumnVector tmp_delta (nn);
83 
84  for (octave_idx_type i = 0; i < nn; i++)
85  {
86  tmp_deriv.elem (i) = deriv[i];
87  tmp_state.elem (i) = state[i];
88  }
89 
90  tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires);
91 
92  if (ires >= 0)
93  {
94  if (tmp_delta.length () == 0)
95  ires = -2;
96  else
97  {
98  for (octave_idx_type i = 0; i < nn; i++)
99  delta[i] = tmp_delta.elem (i);
100  }
101  }
102 
103  END_INTERRUPT_WITH_EXCEPTIONS;
104 
105  return 0;
106 }
107 
108 //NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
109 //C WP, IWP, B, EPLIN, IER, RPAR, IPAR)
110 
111 static octave_idx_type
112 ddaspk_psol (const octave_idx_type&, const double&, const double *,
113  const double *, const double *, const double&,
114  const double *, double *, octave_idx_type *, double *,
115  const double&, octave_idx_type&, double *, octave_idx_type*)
116 {
117  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
118 
119  abort ();
120 
121  END_INTERRUPT_WITH_EXCEPTIONS;
122 
123  return 0;
124 }
125 
126 
127 static octave_idx_type
128 ddaspk_j (const double& time, const double *state, const double *deriv,
129  double *pd, const double& cj, double *, octave_idx_type *)
130 {
131  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
132 
133  // FIXME: would be nice to avoid copying the data.
134 
135  ColumnVector tmp_state (nn);
136  ColumnVector tmp_deriv (nn);
137 
138  for (octave_idx_type i = 0; i < nn; i++)
139  {
140  tmp_deriv.elem (i) = deriv[i];
141  tmp_state.elem (i) = state[i];
142  }
143 
144  Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
145 
146  for (octave_idx_type j = 0; j < nn; j++)
147  for (octave_idx_type i = 0; i < nn; i++)
148  pd[nn * j + i] = tmp_pd.elem (i, j);
149 
150  END_INTERRUPT_WITH_EXCEPTIONS;
151 
152  return 0;
153 }
154 
156 DASPK::do_integrate (double tout)
157 {
158  // FIXME: should handle all this option stuff just once
159  // for each new problem.
160 
161  ColumnVector retval;
162 
164  {
165  integration_error = false;
166 
167  initialized = true;
168 
169  info.resize (dim_vector (20, 1));
170 
171  for (octave_idx_type i = 0; i < 20; i++)
172  info(i) = 0;
173 
174  octave_idx_type n = size ();
175 
176  nn = n;
177 
178  info(0) = 0;
179 
180  if (stop_time_set)
181  {
182  rwork(0) = stop_time;
183  info(3) = 1;
184  }
185  else
186  info(3) = 0;
187 
188  // DAEFunc
189 
192 
193  if (user_fun)
194  {
195  octave_idx_type ires = 0;
196 
197  ColumnVector res = (*user_fun) (x, xdot, t, ires);
198 
199  if (res.length () != x.length ())
200  {
201  (*current_liboctave_error_handler)
202  ("daspk: inconsistent sizes for state and residual vectors");
203 
204  integration_error = true;
205  return retval;
206  }
207  }
208  else
209  {
210  (*current_liboctave_error_handler)
211  ("daspk: no user supplied RHS subroutine!");
212 
213  integration_error = true;
214  return retval;
215  }
216 
217  info(4) = user_jac ? 1 : 0;
218 
219  DAEFunc::reset = false;
220 
224 
225  liw = 40 + n;
226  if (eiq == 1 || eiq == 3)
227  liw += n;
228  if (ccic == 1 || eavfet == 1)
229  liw += n;
230 
231  lrw = 50 + 9*n + n*n;
232  if (eavfet == 1)
233  lrw += n;
234 
235  iwork.resize (dim_vector (liw, 1));
236  rwork.resize (dim_vector (lrw, 1));
237 
238  // DASPK_options
239 
242 
243  octave_idx_type abs_tol_len = abs_tol.length ();
244  octave_idx_type rel_tol_len = rel_tol.length ();
245 
246  if (abs_tol_len == 1 && rel_tol_len == 1)
247  {
248  info(1) = 0;
249  }
250  else if (abs_tol_len == n && rel_tol_len == n)
251  {
252  info(1) = 1;
253  }
254  else
255  {
256  (*current_liboctave_error_handler)
257  ("daspk: inconsistent sizes for tolerance arrays");
258 
259  integration_error = true;
260  return retval;
261  }
262 
263  double hmax = maximum_step_size ();
264  if (hmax >= 0.0)
265  {
266  rwork(1) = hmax;
267  info(6) = 1;
268  }
269  else
270  info(6) = 0;
271 
272  double h0 = initial_step_size ();
273  if (h0 >= 0.0)
274  {
275  rwork(2) = h0;
276  info(7) = 1;
277  }
278  else
279  info(7) = 0;
280 
281  octave_idx_type maxord = maximum_order ();
282  if (maxord >= 0)
283  {
284  if (maxord > 0 && maxord < 6)
285  {
286  info(8) = 1;
287  iwork(2) = maxord;
288  }
289  else
290  {
291  (*current_liboctave_error_handler)
292  ("daspk: invalid value for maximum order");
293  integration_error = true;
294  return retval;
295  }
296  }
297 
298  switch (eiq)
299  {
300  case 1:
301  case 3:
302  {
304 
305  if (ict.length () == n)
306  {
307  for (octave_idx_type i = 0; i < n; i++)
308  {
309  octave_idx_type val = ict(i);
310  if (val < -2 || val > 2)
311  {
312  (*current_liboctave_error_handler)
313  ("daspk: invalid value for inequality constraint type");
314  integration_error = true;
315  return retval;
316  }
317  iwork(40+i) = val;
318  }
319  }
320  else
321  {
322  (*current_liboctave_error_handler)
323  ("daspk: inequality constraint types size mismatch");
324  integration_error = true;
325  return retval;
326  }
327  }
328  // Fall through...
329 
330  case 0:
331  case 2:
332  info(9) = eiq;
333  break;
334 
335  default:
336  (*current_liboctave_error_handler)
337  ("daspk: invalid value for enforce inequality constraints option");
338  integration_error = true;
339  return retval;
340  }
341 
342  if (ccic)
343  {
344  if (ccic == 1)
345  {
346  // FIXME: this code is duplicated below.
347 
349 
350  if (av.length () == n)
351  {
352  octave_idx_type lid;
353  if (eiq == 0 || eiq == 2)
354  lid = 40;
355  else if (eiq == 1 || eiq == 3)
356  lid = 40 + n;
357  else
358  abort ();
359 
360  for (octave_idx_type i = 0; i < n; i++)
361  iwork(lid+i) = av(i) ? -1 : 1;
362  }
363  else
364  {
365  (*current_liboctave_error_handler)
366  ("daspk: algebraic variables size mismatch");
367  integration_error = true;
368  return retval;
369  }
370  }
371  else if (ccic != 2)
372  {
373  (*current_liboctave_error_handler)
374  ("daspk: invalid value for compute consistent initial condition option");
375  integration_error = true;
376  return retval;
377  }
378 
379  info(10) = ccic;
380  }
381 
382  if (eavfet)
383  {
384  info(15) = 1;
385 
386  // FIXME: this code is duplicated above.
387 
389 
390  if (av.length () == n)
391  {
392  octave_idx_type lid;
393  if (eiq == 0 || eiq == 2)
394  lid = 40;
395  else if (eiq == 1 || eiq == 3)
396  lid = 40 + n;
397  else
398  abort ();
399 
400  for (octave_idx_type i = 0; i < n; i++)
401  iwork(lid+i) = av(i) ? -1 : 1;
402  }
403  }
404 
406  {
408 
409  if (ich.length () == 6)
410  {
411  iwork(31) = NINTbig (ich(0));
412  iwork(32) = NINTbig (ich(1));
413  iwork(33) = NINTbig (ich(2));
414  iwork(34) = NINTbig (ich(3));
415 
416  rwork(13) = ich(4);
417  rwork(14) = ich(5);
418  }
419  else
420  {
421  (*current_liboctave_error_handler)
422  ("daspk: invalid initial condition heuristics option");
423  integration_error = true;
424  return retval;
425  }
426 
427  info(16) = 1;
428  }
429 
431  switch (pici)
432  {
433  case 0:
434  case 1:
435  case 2:
436  info(17) = pici;
437  break;
438 
439  default:
440  (*current_liboctave_error_handler)
441  ("daspk: invalid value for print initial condition info option");
442  integration_error = true;
443  return retval;
444  break;
445  }
446 
447  DASPK_options::reset = false;
448 
449  restart = false;
450  }
451 
452  double *px = x.fortran_vec ();
453  double *pxdot = xdot.fortran_vec ();
454 
455  octave_idx_type *pinfo = info.fortran_vec ();
456 
457  double *prel_tol = rel_tol.fortran_vec ();
458  double *pabs_tol = abs_tol.fortran_vec ();
459 
460  double *prwork = rwork.fortran_vec ();
461  octave_idx_type *piwork = iwork.fortran_vec ();
462 
463  double *dummy = 0;
464  octave_idx_type *idummy = 0;
465 
466  F77_XFCN (ddaspk, DDASPK, (ddaspk_f, nn, t, px, pxdot, tout, pinfo,
467  prel_tol, pabs_tol, istate, prwork, lrw,
468  piwork, liw, dummy, idummy, ddaspk_j,
469  ddaspk_psol));
470 
471  switch (istate)
472  {
473  case 1: // A step was successfully taken in intermediate-output
474  // mode. The code has not yet reached TOUT.
475  case 2: // The integration to TSTOP was successfully completed
476  // (T=TSTOP) by stepping exactly to TSTOP.
477  case 3: // The integration to TOUT was successfully completed
478  // (T=TOUT) by stepping past TOUT. Y(*) is obtained by
479  // interpolation. YPRIME(*) is obtained by interpolation.
480  case 4: // The initial condition calculation, with
481  // INFO(11) > 0, was successful, and INFO(14) = 1.
482  // No integration steps were taken, and the solution
483  // is not considered to have been started.
484  retval = x;
485  t = tout;
486  break;
487 
488  case -1: // A large amount of work has been expended. (~500 steps).
489  case -2: // The error tolerances are too stringent.
490  case -3: // The local error test cannot be satisfied because you
491  // specified a zero component in ATOL and the
492  // corresponding computed solution component is zero.
493  // Thus, a pure relative error test is impossible for
494  // this component.
495  case -6: // DDASPK had repeated error test failures on the last
496  // attempted step.
497  case -7: // The corrector could not converge.
498  case -8: // The matrix of partial derivatives is singular.
499  case -9: // The corrector could not converge. There were repeated
500  // error test failures in this step.
501  case -10: // The corrector could not converge because IRES was
502  // equal to minus one.
503  case -11: // IRES equal to -2 was encountered and control is being
504  // returned to the calling program.
505  case -12: // DDASPK failed to compute the initial YPRIME.
506  case -13: // Unrecoverable error encountered inside user's
507  // PSOL routine, and control is being returned to
508  // the calling program.
509  case -14: // The Krylov linear system solver could not
510  // achieve convergence.
511  case -33: // The code has encountered trouble from which it cannot
512  // recover. A message is printed explaining the trouble
513  // and control is returned to the calling program. For
514  // example, this occurs when invalid input is detected.
515  integration_error = true;
516  break;
517 
518  default:
519  integration_error = true;
520  (*current_liboctave_error_handler)
521  ("unrecognized value of istate (= %d) returned from ddaspk",
522  istate);
523  break;
524  }
525 
526  return retval;
527 }
528 
529 Matrix
531 {
532  Matrix dummy;
533  return integrate (tout, dummy);
534 }
535 
536 Matrix
537 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out)
538 {
539  Matrix retval;
540 
541  octave_idx_type n_out = tout.capacity ();
542  octave_idx_type n = size ();
543 
544  if (n_out > 0 && n > 0)
545  {
546  retval.resize (n_out, n);
547  xdot_out.resize (n_out, n);
548 
549  for (octave_idx_type i = 0; i < n; i++)
550  {
551  retval.elem (0, i) = x.elem (i);
552  xdot_out.elem (0, i) = xdot.elem (i);
553  }
554 
555  for (octave_idx_type j = 1; j < n_out; j++)
556  {
557  ColumnVector x_next = do_integrate (tout.elem (j));
558 
559  if (integration_error)
560  return retval;
561 
562  for (octave_idx_type i = 0; i < n; i++)
563  {
564  retval.elem (j, i) = x_next.elem (i);
565  xdot_out.elem (j, i) = xdot.elem (i);
566  }
567  }
568  }
569 
570  return retval;
571 }
572 
573 Matrix
574 DASPK::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
575 {
576  Matrix dummy;
577  return integrate (tout, dummy, tcrit);
578 }
579 
580 Matrix
581 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out,
582  const ColumnVector& tcrit)
583 {
584  Matrix retval;
585 
586  octave_idx_type n_out = tout.capacity ();
587  octave_idx_type n = size ();
588 
589  if (n_out > 0 && n > 0)
590  {
591  retval.resize (n_out, n);
592  xdot_out.resize (n_out, n);
593 
594  for (octave_idx_type i = 0; i < n; i++)
595  {
596  retval.elem (0, i) = x.elem (i);
597  xdot_out.elem (0, i) = xdot.elem (i);
598  }
599 
600  octave_idx_type n_crit = tcrit.capacity ();
601 
602  if (n_crit > 0)
603  {
604  octave_idx_type i_crit = 0;
605  octave_idx_type i_out = 1;
606  double next_crit = tcrit.elem (0);
607  double next_out;
608  while (i_out < n_out)
609  {
610  bool do_restart = false;
611 
612  next_out = tout.elem (i_out);
613  if (i_crit < n_crit)
614  next_crit = tcrit.elem (i_crit);
615 
616  bool save_output;
617  double t_out;
618 
619  if (next_crit == next_out)
620  {
621  set_stop_time (next_crit);
622  t_out = next_out;
623  save_output = true;
624  i_out++;
625  i_crit++;
626  do_restart = true;
627  }
628  else if (next_crit < next_out)
629  {
630  if (i_crit < n_crit)
631  {
632  set_stop_time (next_crit);
633  t_out = next_crit;
634  save_output = false;
635  i_crit++;
636  do_restart = true;
637  }
638  else
639  {
640  clear_stop_time ();
641  t_out = next_out;
642  save_output = true;
643  i_out++;
644  }
645  }
646  else
647  {
648  set_stop_time (next_crit);
649  t_out = next_out;
650  save_output = true;
651  i_out++;
652  }
653 
654  ColumnVector x_next = do_integrate (t_out);
655 
656  if (integration_error)
657  return retval;
658 
659  if (save_output)
660  {
661  for (octave_idx_type i = 0; i < n; i++)
662  {
663  retval.elem (i_out-1, i) = x_next.elem (i);
664  xdot_out.elem (i_out-1, i) = xdot.elem (i);
665  }
666  }
667 
668  if (do_restart)
669  force_restart ();
670  }
671  }
672  else
673  {
674  retval = integrate (tout, xdot_out);
675 
676  if (integration_error)
677  return retval;
678  }
679  }
680 
681  return retval;
682 }
683 
684 std::string
686 {
687  std::string retval;
688 
689  std::ostringstream buf;
690  buf << t;
691  std::string t_curr = buf.str ();
692 
693  switch (istate)
694  {
695  case 1:
696  retval = "a step was successfully taken in intermediate-output mode.";
697  break;
698 
699  case 2:
700  retval = "integration completed by stepping exactly to TOUT";
701  break;
702 
703  case 3:
704  retval = "integration to tout completed by stepping past TOUT";
705  break;
706 
707  case 4:
708  retval = "initial condition calculation completed successfully";
709  break;
710 
711  case -1:
712  retval = std::string ("a large amount of work has been expended (t =")
713  + t_curr + ")";
714  break;
715 
716  case -2:
717  retval = "the error tolerances are too stringent";
718  break;
719 
720  case -3:
721  retval = std::string ("error weight became zero during problem. (t = ")
722  + t_curr
723  + "; solution component i vanished, and atol or atol(i) == 0)";
724  break;
725 
726  case -6:
727  retval = std::string ("repeated error test failures on the last attempted step (t = ")
728  + t_curr + ")";
729  break;
730 
731  case -7:
732  retval = std::string ("the corrector could not converge (t = ")
733  + t_curr + ")";
734  break;
735 
736  case -8:
737  retval = std::string ("the matrix of partial derivatives is singular (t = ")
738  + t_curr + ")";
739  break;
740 
741  case -9:
742  retval = std::string ("the corrector could not converge (t = ")
743  + t_curr + "; repeated test failures)";
744  break;
745 
746  case -10:
747  retval = std::string ("corrector could not converge because IRES was -1 (t = ")
748  + t_curr + ")";
749  break;
750 
751  case -11:
752  retval = std::string ("return requested in user-supplied function (t = ")
753  + t_curr + ")";
754  break;
755 
756  case -12:
757  retval = "failed to compute consistent initial conditions";
758  break;
759 
760  case -13:
761  retval = std::string ("unrecoverable error encountered inside user's PSOL function (t = ")
762  + t_curr + ")";
763  break;
764 
765  case -14:
766  retval = std::string ("the Krylov linear system solver failed to converge (t = ")
767  + t_curr + ")";
768  break;
769 
770  case -33:
771  retval = "unrecoverable error (see printed message)";
772  break;
773 
774  default:
775  retval = "unknown error state";
776  break;
777  }
778 
779  return retval;
780 }