00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include <cfloat>
00032 #include <cmath>
00033
00034 #include "DASPK.h"
00035 #include "f77-fcn.h"
00036 #include "lo-error.h"
00037 #include "lo-sstream.h"
00038 #include "quit.h"
00039
00040 typedef int (*daspk_fcn_ptr) (const double&, const double*,
00041 const double*, const double&,
00042 double*, int&, double*, int*);
00043
00044 typedef int (*daspk_jac_ptr) (const double&, const double*,
00045 const double*, double*,
00046 const double&, double*, int*);
00047
00048 typedef int (*daspk_psol_ptr) (const int&, const double&,
00049 const double*, const double*,
00050 const double*, const double&,
00051 const double*, double*, int*,
00052 double*, const double&, int&,
00053 double*, int*);
00054
00055 extern "C"
00056 {
00057 F77_RET_T
00058 F77_FUNC (ddaspk, DDASPK) (daspk_fcn_ptr, const int&, double&,
00059 double*, double*, double&, const int*,
00060 const double*, const double*, int&,
00061 double*, const int&, int*, const int&,
00062 const double*, const int*,
00063 daspk_jac_ptr, daspk_psol_ptr);
00064 }
00065
00066 static DAEFunc::DAERHSFunc user_fun;
00067 static DAEFunc::DAEJacFunc user_jac;
00068 static int nn;
00069
00070 static int
00071 ddaspk_f (const double& time, const double *state, const double *deriv,
00072 const double&, double *delta, int& ires, double *, int *)
00073 {
00074 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00075
00076 ColumnVector tmp_deriv (nn);
00077 ColumnVector tmp_state (nn);
00078 ColumnVector tmp_delta (nn);
00079
00080 for (int i = 0; i < nn; i++)
00081 {
00082 tmp_deriv.elem (i) = deriv [i];
00083 tmp_state.elem (i) = state [i];
00084 }
00085
00086 tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires);
00087
00088 if (ires >= 0)
00089 {
00090 if (tmp_delta.length () == 0)
00091 ires = -2;
00092 else
00093 {
00094 for (int i = 0; i < nn; i++)
00095 delta [i] = tmp_delta.elem (i);
00096 }
00097 }
00098
00099 END_INTERRUPT_WITH_EXCEPTIONS;
00100
00101 return 0;
00102 }
00103
00104
00105
00106
00107 static int
00108 ddaspk_psol (const int&, const double&, const double *,
00109 const double *, const double *, const double&,
00110 const double *, double *, int *, double *,
00111 const double&, int&, double *, int*)
00112 {
00113 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00114
00115 abort ();
00116
00117 END_INTERRUPT_WITH_EXCEPTIONS;
00118
00119 return 0;
00120 }
00121
00122
00123 static int
00124 ddaspk_j (const double& time, const double *state, const double *deriv,
00125 double *pd, const double& cj, double *, int *)
00126 {
00127 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00128
00129
00130
00131 ColumnVector tmp_state (nn);
00132 ColumnVector tmp_deriv (nn);
00133
00134 for (int i = 0; i < nn; i++)
00135 {
00136 tmp_deriv.elem (i) = deriv [i];
00137 tmp_state.elem (i) = state [i];
00138 }
00139
00140 Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
00141
00142 for (int j = 0; j < nn; j++)
00143 for (int i = 0; i < nn; i++)
00144 pd [nn * j + i] = tmp_pd.elem (i, j);
00145
00146 END_INTERRUPT_WITH_EXCEPTIONS;
00147
00148 return 0;
00149 }
00150
00151 ColumnVector
00152 DASPK::do_integrate (double tout)
00153 {
00154
00155
00156
00157 ColumnVector retval;
00158
00159 if (! initialized || restart || DAEFunc::reset|| DASPK_options::reset)
00160 {
00161 integration_error = false;
00162
00163 initialized = true;
00164
00165 info.resize (20);
00166
00167 for (int i = 0; i < 20; i++)
00168 info(i) = 0;
00169
00170 pinfo = info.fortran_vec ();
00171
00172 int n = size ();
00173
00174 nn = n;
00175
00176 info(0) = 0;
00177
00178 if (stop_time_set)
00179 {
00180 rwork(0) = stop_time;
00181 info(3) = 1;
00182 }
00183 else
00184 info(3) = 0;
00185
00186 px = x.fortran_vec ();
00187 pxdot = xdot.fortran_vec ();
00188
00189
00190
00191 user_fun = DAEFunc::function ();
00192 user_jac = DAEFunc::jacobian_function ();
00193
00194 if (user_fun)
00195 {
00196 int ires = 0;
00197
00198 ColumnVector res = (*user_fun) (x, xdot, t, ires);
00199
00200 if (res.length () != x.length ())
00201 {
00202 (*current_liboctave_error_handler)
00203 ("daspk: inconsistent sizes for state and residual vectors");
00204
00205 integration_error = true;
00206 return retval;
00207 }
00208 }
00209 else
00210 {
00211 (*current_liboctave_error_handler)
00212 ("daspk: no user supplied RHS subroutine!");
00213
00214 integration_error = true;
00215 return retval;
00216 }
00217
00218 info(4) = user_jac ? 1 : 0;
00219
00220 DAEFunc::reset = false;
00221
00222 int eiq = enforce_inequality_constraints ();
00223 int ccic = compute_consistent_initial_condition ();
00224 int eavfet = exclude_algebraic_variables_from_error_test ();
00225
00226 liw = 40 + n;
00227 if (eiq == 1 || eiq == 3)
00228 liw += n;
00229 if (ccic == 1 || eavfet == 1)
00230 liw += n;
00231
00232 lrw = 50 + 9*n + n*n;
00233 if (eavfet == 1)
00234 lrw += n;
00235
00236 iwork.resize (liw);
00237 rwork.resize (lrw);
00238
00239 piwork = iwork.fortran_vec ();
00240 prwork = rwork.fortran_vec ();
00241
00242
00243
00244 abs_tol = absolute_tolerance ();
00245 rel_tol = relative_tolerance ();
00246
00247 int abs_tol_len = abs_tol.length ();
00248 int rel_tol_len = rel_tol.length ();
00249
00250 if (abs_tol_len == 1 && rel_tol_len == 1)
00251 {
00252 info(1) = 0;
00253 }
00254 else if (abs_tol_len == n && rel_tol_len == n)
00255 {
00256 info(1) = 1;
00257 }
00258 else
00259 {
00260 (*current_liboctave_error_handler)
00261 ("daspk: inconsistent sizes for tolerance arrays");
00262
00263 integration_error = true;
00264 return retval;
00265 }
00266
00267 pabs_tol = abs_tol.fortran_vec ();
00268 prel_tol = rel_tol.fortran_vec ();
00269
00270 double hmax = maximum_step_size ();
00271 if (hmax >= 0.0)
00272 {
00273 rwork(1) = hmax;
00274 info(6) = 1;
00275 }
00276 else
00277 info(6) = 0;
00278
00279 double h0 = initial_step_size ();
00280 if (h0 >= 0.0)
00281 {
00282 rwork(2) = h0;
00283 info(7) = 1;
00284 }
00285 else
00286 info(7) = 0;
00287
00288 int maxord = maximum_order ();
00289 if (maxord >= 0)
00290 {
00291 if (maxord > 0 && maxord < 6)
00292 {
00293 info(8) = 1;
00294 iwork(2) = maxord;
00295 }
00296 else
00297 {
00298 (*current_liboctave_error_handler)
00299 ("daspk: invalid value for maximum order");
00300 integration_error = true;
00301 return retval;
00302 }
00303 }
00304
00305 switch (eiq)
00306 {
00307 case 1:
00308 case 3:
00309 {
00310 Array<int> ict = inequality_constraint_types ();
00311
00312 if (ict.length () == n)
00313 {
00314 for (int i = 0; i < n; i++)
00315 {
00316 int val = ict(i);
00317 if (val < -2 || val > 2)
00318 {
00319 (*current_liboctave_error_handler)
00320 ("daspk: invalid value for inequality constraint type");
00321 integration_error = true;
00322 return retval;
00323 }
00324 iwork(40+i) = val;
00325 }
00326 }
00327 else
00328 {
00329 (*current_liboctave_error_handler)
00330 ("daspk: inequality constraint types size mismatch");
00331 integration_error = true;
00332 return retval;
00333 }
00334 }
00335
00336
00337 case 0:
00338 case 2:
00339 info(9) = eiq;
00340 break;
00341
00342 default:
00343 (*current_liboctave_error_handler)
00344 ("daspk: invalid value for enforce inequality constraints option");
00345 integration_error = true;
00346 return retval;
00347 }
00348
00349 if (ccic)
00350 {
00351 if (ccic == 1)
00352 {
00353
00354
00355 Array<int> av = algebraic_variables ();
00356
00357 if (av.length () == n)
00358 {
00359 int lid;
00360 if (eiq == 0 || eiq == 2)
00361 lid = 40;
00362 else if (eiq == 1 || eiq == 3)
00363 lid = 40 + n;
00364 else
00365 abort ();
00366
00367 for (int i = 0; i < n; i++)
00368 iwork(lid+i) = av(i) ? -1 : 1;
00369 }
00370 else
00371 {
00372 (*current_liboctave_error_handler)
00373 ("daspk: algebraic variables size mismatch");
00374 integration_error = true;
00375 return retval;
00376 }
00377 }
00378 else if (ccic != 2)
00379 {
00380 (*current_liboctave_error_handler)
00381 ("daspk: invalid value for compute consistent initial condition option");
00382 integration_error = true;
00383 return retval;
00384 }
00385
00386 info(10) = ccic;
00387 }
00388
00389 if (eavfet)
00390 {
00391 info(15) = 1;
00392
00393
00394
00395 Array<int> av = algebraic_variables ();
00396
00397 if (av.length () == n)
00398 {
00399 int lid;
00400 if (eiq == 0 || eiq == 2)
00401 lid = 40;
00402 else if (eiq == 1 || eiq == 3)
00403 lid = 40 + n;
00404 else
00405 abort ();
00406
00407 for (int i = 0; i < n; i++)
00408 iwork(lid+i) = av(i) ? -1 : 1;
00409 }
00410 }
00411
00412 if (use_initial_condition_heuristics ())
00413 {
00414 Array<double> ich = initial_condition_heuristics ();
00415
00416 if (ich.length () == 6)
00417 {
00418 iwork(31) = NINT (ich(0));
00419 iwork(32) = NINT (ich(1));
00420 iwork(33) = NINT (ich(2));
00421 iwork(34) = NINT (ich(3));
00422
00423 rwork(13) = ich(4);
00424 rwork(14) = ich(5);
00425 }
00426 else
00427 {
00428 (*current_liboctave_error_handler)
00429 ("daspk: invalid initial condition heuristics option");
00430 integration_error = true;
00431 return retval;
00432 }
00433
00434 info(16) = 1;
00435 }
00436
00437 int pici = print_initial_condition_info ();
00438 switch (pici)
00439 {
00440 case 0:
00441 case 1:
00442 case 2:
00443 info(17) = pici;
00444 break;
00445
00446 default:
00447 (*current_liboctave_error_handler)
00448 ("daspk: invalid value for print initial condition info option");
00449 integration_error = true;
00450 return retval;
00451 break;
00452 }
00453
00454 DASPK_options::reset = false;
00455
00456 restart = false;
00457 }
00458
00459 static double *dummy = 0;
00460 static int *idummy = 0;
00461
00462 F77_XFCN (ddaspk, DDASPK, (ddaspk_f, nn, t, px, pxdot, tout, pinfo,
00463 prel_tol, pabs_tol, istate, prwork, lrw,
00464 piwork, liw, dummy, idummy, ddaspk_j,
00465 ddaspk_psol));
00466
00467 if (f77_exception_encountered)
00468 {
00469 integration_error = true;
00470 (*current_liboctave_error_handler) ("unrecoverable error in daspk");
00471 }
00472 else
00473 {
00474 switch (istate)
00475 {
00476 case 1:
00477
00478 case 2:
00479
00480 case 3:
00481
00482
00483 case 4:
00484
00485
00486
00487 retval = x;
00488 t = tout;
00489 break;
00490
00491 case -1:
00492 case -2:
00493 case -3:
00494
00495
00496
00497
00498 case -6:
00499
00500 case -7:
00501 case -8:
00502 case -9:
00503
00504 case -10:
00505
00506 case -11:
00507
00508 case -12:
00509 case -13:
00510
00511
00512 case -14:
00513
00514 case -33:
00515
00516
00517
00518 integration_error = true;
00519 break;
00520
00521 default:
00522 integration_error = true;
00523 (*current_liboctave_error_handler)
00524 ("unrecognized value of istate (= %d) returned from ddaspk",
00525 istate);
00526 break;
00527 }
00528 }
00529
00530 return retval;
00531 }
00532
00533 Matrix
00534 DASPK::do_integrate (const ColumnVector& tout)
00535 {
00536 Matrix dummy;
00537 return integrate (tout, dummy);
00538 }
00539
00540 Matrix
00541 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out)
00542 {
00543 Matrix retval;
00544
00545 int n_out = tout.capacity ();
00546 int n = size ();
00547
00548 if (n_out > 0 && n > 0)
00549 {
00550 retval.resize (n_out, n);
00551 xdot_out.resize (n_out, n);
00552
00553 for (int i = 0; i < n; i++)
00554 {
00555 retval.elem (0, i) = x.elem (i);
00556 xdot_out.elem (0, i) = xdot.elem (i);
00557 }
00558
00559 for (int j = 1; j < n_out; j++)
00560 {
00561 ColumnVector x_next = do_integrate (tout.elem (j));
00562
00563 if (integration_error)
00564 return retval;
00565
00566 for (int i = 0; i < n; i++)
00567 {
00568 retval.elem (j, i) = x_next.elem (i);
00569 xdot_out.elem (j, i) = xdot.elem (i);
00570 }
00571 }
00572 }
00573
00574 return retval;
00575 }
00576
00577 Matrix
00578 DASPK::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
00579 {
00580 Matrix dummy;
00581 return integrate (tout, dummy, tcrit);
00582 }
00583
00584 Matrix
00585 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out,
00586 const ColumnVector& tcrit)
00587 {
00588 Matrix retval;
00589
00590 int n_out = tout.capacity ();
00591 int n = size ();
00592
00593 if (n_out > 0 && n > 0)
00594 {
00595 retval.resize (n_out, n);
00596 xdot_out.resize (n_out, n);
00597
00598 for (int i = 0; i < n; i++)
00599 {
00600 retval.elem (0, i) = x.elem (i);
00601 xdot_out.elem (0, i) = xdot.elem (i);
00602 }
00603
00604 int n_crit = tcrit.capacity ();
00605
00606 if (n_crit > 0)
00607 {
00608 int i_crit = 0;
00609 int i_out = 1;
00610 double next_crit = tcrit.elem (0);
00611 double next_out;
00612 while (i_out < n_out)
00613 {
00614 bool do_restart = false;
00615
00616 next_out = tout.elem (i_out);
00617 if (i_crit < n_crit)
00618 next_crit = tcrit.elem (i_crit);
00619
00620 bool save_output;
00621 double t_out;
00622
00623 if (next_crit == next_out)
00624 {
00625 set_stop_time (next_crit);
00626 t_out = next_out;
00627 save_output = true;
00628 i_out++;
00629 i_crit++;
00630 do_restart = true;
00631 }
00632 else if (next_crit < next_out)
00633 {
00634 if (i_crit < n_crit)
00635 {
00636 set_stop_time (next_crit);
00637 t_out = next_crit;
00638 save_output = false;
00639 i_crit++;
00640 do_restart = true;
00641 }
00642 else
00643 {
00644 clear_stop_time ();
00645 t_out = next_out;
00646 save_output = true;
00647 i_out++;
00648 }
00649 }
00650 else
00651 {
00652 set_stop_time (next_crit);
00653 t_out = next_out;
00654 save_output = true;
00655 i_out++;
00656 }
00657
00658 ColumnVector x_next = do_integrate (t_out);
00659
00660 if (integration_error)
00661 return retval;
00662
00663 if (save_output)
00664 {
00665 for (int i = 0; i < n; i++)
00666 {
00667 retval.elem (i_out-1, i) = x_next.elem (i);
00668 xdot_out.elem (i_out-1, i) = xdot.elem (i);
00669 }
00670 }
00671
00672 if (do_restart)
00673 force_restart ();
00674 }
00675 }
00676 else
00677 {
00678 retval = integrate (tout, xdot_out);
00679
00680 if (integration_error)
00681 return retval;
00682 }
00683 }
00684
00685 return retval;
00686 }
00687
00688 std::string
00689 DASPK::error_message (void) const
00690 {
00691 std::string retval;
00692
00693 OSSTREAM buf;
00694 buf << t << OSSTREAM_ENDS;
00695 std::string t_curr = OSSTREAM_STR (buf);
00696 OSSTREAM_FREEZE (buf);
00697
00698 switch (istate)
00699 {
00700 case 1:
00701 retval = "a step was successfully taken in intermediate-output mode.";
00702 break;
00703
00704 case 2:
00705 retval = "integration completed by stepping exactly to TOUT";
00706 break;
00707
00708 case 3:
00709 retval = "integration to tout completed by stepping past TOUT";
00710 break;
00711
00712 case 4:
00713 retval = "initial condition calculation completed successfully";
00714 break;
00715
00716 case -1:
00717 retval = std::string ("a large amount of work has been expended (t =")
00718 + t_curr + ")";
00719 break;
00720
00721 case -2:
00722 retval = "the error tolerances are too stringent";
00723 break;
00724
00725 case -3:
00726 retval = std::string ("error weight became zero during problem. (t = ")
00727 + t_curr
00728 + "; solution component i vanished, and atol or atol(i) == 0)";
00729 break;
00730
00731 case -6:
00732 retval = std::string ("repeated error test failures on the last attempted step (t = ")
00733 + t_curr + ")";
00734 break;
00735
00736 case -7:
00737 retval = std::string ("the corrector could not converge (t = ")
00738 + t_curr + ")";
00739 break;
00740
00741 case -8:
00742 retval = std::string ("the matrix of partial derivatives is singular (t = ")
00743 + t_curr + ")";
00744 break;
00745
00746 case -9:
00747 retval = std::string ("the corrector could not converge (t = ")
00748 + t_curr + "; repeated test failures)";
00749 break;
00750
00751 case -10:
00752 retval = std::string ("corrector could not converge because IRES was -1 (t = ")
00753 + t_curr + ")";
00754 break;
00755
00756 case -11:
00757 retval = std::string ("return requested in user-supplied function (t = ")
00758 + t_curr + ")";
00759 break;
00760
00761 case -12:
00762 retval = "failed to compute consistent initial conditions";
00763 break;
00764
00765 case -13:
00766 retval = std::string ("unrecoverable error encountered inside user's PSOL function (t = ")
00767 + t_curr + ")";
00768 break;
00769
00770 case -14:
00771 retval = std::string ("the Krylov linear system solver failed to converge (t = ")
00772 + t_curr + ")";
00773 break;
00774
00775 case -33:
00776 retval = "unrecoverable error (see printed message)";
00777 break;
00778
00779 default:
00780 retval = "unknown error state";
00781 break;
00782 }
00783
00784 return retval;
00785 }
00786
00787
00788
00789
00790
00791