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 "DASSL.h"
00035 #include "f77-fcn.h"
00036 #include "lo-error.h"
00037 #include "lo-sstream.h"
00038 #include "quit.h"
00039
00040 typedef int (*dassl_fcn_ptr) (const double&, const double*, const double*,
00041 double*, int&, double*, int*);
00042
00043 typedef int (*dassl_jac_ptr) (const double&, const double*, const double*,
00044 double*, const double&, double*, int*);
00045
00046 extern "C"
00047 {
00048 F77_RET_T
00049 F77_FUNC (ddassl, DDASSL) (dassl_fcn_ptr, const int&, double&,
00050 double*, double*, double&, const int*,
00051 const double*, const double*, int&,
00052 double*, const int&, int*, const int&,
00053 const double*, const int*,
00054 dassl_jac_ptr);
00055 }
00056
00057 static DAEFunc::DAERHSFunc user_fun;
00058 static DAEFunc::DAEJacFunc user_jac;
00059
00060 static int nn;
00061
00062 static int
00063 ddassl_f (const double& time, const double *state, const double *deriv,
00064 double *delta, int& ires, double *, int *)
00065 {
00066 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00067
00068
00069
00070 ColumnVector tmp_deriv (nn);
00071 ColumnVector tmp_state (nn);
00072 ColumnVector tmp_delta (nn);
00073
00074 for (int i = 0; i < nn; i++)
00075 {
00076 tmp_deriv.elem (i) = deriv [i];
00077 tmp_state.elem (i) = state [i];
00078 }
00079
00080 tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires);
00081
00082 if (ires >= 0)
00083 {
00084 if (tmp_delta.length () == 0)
00085 ires = -2;
00086 else
00087 {
00088 for (int i = 0; i < nn; i++)
00089 delta [i] = tmp_delta.elem (i);
00090 }
00091 }
00092
00093 END_INTERRUPT_WITH_EXCEPTIONS;
00094
00095 return 0;
00096 }
00097
00098 static int
00099 ddassl_j (const double& time, const double *state, const double *deriv,
00100 double *pd, const double& cj, double *, int *)
00101 {
00102 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00103
00104
00105
00106 ColumnVector tmp_state (nn);
00107 ColumnVector tmp_deriv (nn);
00108
00109 for (int i = 0; i < nn; i++)
00110 {
00111 tmp_deriv.elem (i) = deriv [i];
00112 tmp_state.elem (i) = state [i];
00113 }
00114
00115 Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj);
00116
00117 for (int j = 0; j < nn; j++)
00118 for (int i = 0; i < nn; i++)
00119 pd [nn * j + i] = tmp_pd.elem (i, j);
00120
00121 END_INTERRUPT_WITH_EXCEPTIONS;
00122
00123 return 0;
00124 }
00125
00126 ColumnVector
00127 DASSL::do_integrate (double tout)
00128 {
00129 ColumnVector retval;
00130
00131 if (! initialized || restart || DAEFunc::reset|| DASSL_options::reset)
00132 {
00133 integration_error = false;
00134
00135 initialized = true;
00136
00137 info.resize (15);
00138
00139 for (int i = 0; i < 15; i++)
00140 info(i) = 0;
00141
00142 pinfo = info.fortran_vec ();
00143
00144 int n = size ();
00145
00146 liw = 21 + n;
00147 lrw = 40 + 9*n + n*n;
00148
00149 nn = n;
00150
00151 iwork.resize (liw);
00152 rwork.resize (lrw);
00153
00154 info(0) = 0;
00155
00156 if (stop_time_set)
00157 {
00158 rwork(0) = stop_time;
00159 info(3) = 1;
00160 }
00161 else
00162 info(3) = 0;
00163
00164 px = x.fortran_vec ();
00165 pxdot = xdot.fortran_vec ();
00166
00167 piwork = iwork.fortran_vec ();
00168 prwork = rwork.fortran_vec ();
00169
00170 restart = false;
00171
00172
00173
00174 user_fun = DAEFunc::function ();
00175 user_jac = DAEFunc::jacobian_function ();
00176
00177 if (user_fun)
00178 {
00179 int ires = 0;
00180
00181 ColumnVector res = (*user_fun) (x, xdot, t, ires);
00182
00183 if (res.length () != x.length ())
00184 {
00185 (*current_liboctave_error_handler)
00186 ("dassl: inconsistent sizes for state and residual vectors");
00187
00188 integration_error = true;
00189 return retval;
00190 }
00191 }
00192 else
00193 {
00194 (*current_liboctave_error_handler)
00195 ("dassl: no user supplied RHS subroutine!");
00196
00197 integration_error = true;
00198 return retval;
00199 }
00200
00201 info(4) = user_jac ? 1 : 0;
00202
00203 DAEFunc::reset = false;
00204
00205
00206
00207 double hmax = maximum_step_size ();
00208 if (hmax >= 0.0)
00209 {
00210 rwork(1) = hmax;
00211 info(6) = 1;
00212 }
00213 else
00214 info(6) = 0;
00215
00216 double h0 = initial_step_size ();
00217 if (h0 >= 0.0)
00218 {
00219 rwork(2) = h0;
00220 info(7) = 1;
00221 }
00222 else
00223 info(7) = 0;
00224
00225 if (step_limit () >= 0)
00226 {
00227 info(11) = 1;
00228 iwork(20) = step_limit ();
00229 }
00230 else
00231 info(11) = 0;
00232
00233 int maxord = maximum_order ();
00234 if (maxord >= 0)
00235 {
00236 if (maxord > 0 && maxord < 6)
00237 {
00238 info(8) = 1;
00239 iwork(2) = maxord;
00240 }
00241 else
00242 {
00243 (*current_liboctave_error_handler)
00244 ("dassl: invalid value for maximum order");
00245 integration_error = true;
00246 return retval;
00247 }
00248 }
00249
00250 int enc = enforce_nonnegativity_constraints ();
00251 info(9) = enc ? 1 : 0;
00252
00253 int ccic = compute_consistent_initial_condition ();
00254 info(10) = ccic ? 1 : 0;
00255
00256 abs_tol = absolute_tolerance ();
00257 rel_tol = relative_tolerance ();
00258
00259 int abs_tol_len = abs_tol.length ();
00260 int rel_tol_len = rel_tol.length ();
00261
00262 if (abs_tol_len == 1 && rel_tol_len == 1)
00263 {
00264 info(1) = 0;
00265 }
00266 else if (abs_tol_len == n && rel_tol_len == n)
00267 {
00268 info(1) = 1;
00269 }
00270 else
00271 {
00272 (*current_liboctave_error_handler)
00273 ("dassl: inconsistent sizes for tolerance arrays");
00274
00275 integration_error = true;
00276 return retval;
00277 }
00278
00279 pabs_tol = abs_tol.fortran_vec ();
00280 prel_tol = rel_tol.fortran_vec ();
00281
00282 DASSL_options::reset = false;
00283 }
00284
00285 static double *dummy = 0;
00286 static int *idummy = 0;
00287
00288 F77_XFCN (ddassl, DDASSL, (ddassl_f, nn, t, px, pxdot, tout, pinfo,
00289 prel_tol, pabs_tol, istate, prwork, lrw,
00290 piwork, liw, dummy, idummy, ddassl_j));
00291
00292 if (f77_exception_encountered)
00293 {
00294 integration_error = true;
00295 (*current_liboctave_error_handler) ("unrecoverable error in dassl");
00296 }
00297 else
00298 {
00299 switch (istate)
00300 {
00301 case 1:
00302
00303 case 2:
00304
00305 case 3:
00306
00307
00308 retval = x;
00309 t = tout;
00310 break;
00311
00312 case -1:
00313 case -2:
00314 case -3:
00315
00316
00317
00318
00319 case -6:
00320
00321 case -7:
00322 case -8:
00323 case -9:
00324
00325 case -10:
00326
00327 case -11:
00328
00329 case -12:
00330 case -33:
00331
00332
00333
00334 integration_error = true;
00335 break;
00336
00337 default:
00338 integration_error = true;
00339 (*current_liboctave_error_handler)
00340 ("unrecognized value of istate (= %d) returned from ddassl",
00341 istate);
00342 break;
00343 }
00344 }
00345
00346 return retval;
00347 }
00348
00349 Matrix
00350 DASSL::do_integrate (const ColumnVector& tout)
00351 {
00352 Matrix dummy;
00353 return integrate (tout, dummy);
00354 }
00355
00356 Matrix
00357 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
00358 {
00359 Matrix retval;
00360
00361 int n_out = tout.capacity ();
00362 int n = size ();
00363
00364 if (n_out > 0 && n > 0)
00365 {
00366 retval.resize (n_out, n);
00367 xdot_out.resize (n_out, n);
00368
00369 for (int i = 0; i < n; i++)
00370 {
00371 retval.elem (0, i) = x.elem (i);
00372 xdot_out.elem (0, i) = xdot.elem (i);
00373 }
00374
00375 for (int j = 1; j < n_out; j++)
00376 {
00377 ColumnVector x_next = do_integrate (tout.elem (j));
00378
00379 if (integration_error)
00380 return retval;
00381
00382 for (int i = 0; i < n; i++)
00383 {
00384 retval.elem (j, i) = x_next.elem (i);
00385 xdot_out.elem (j, i) = xdot.elem (i);
00386 }
00387 }
00388 }
00389
00390 return retval;
00391 }
00392
00393 Matrix
00394 DASSL::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
00395 {
00396 Matrix dummy;
00397 return integrate (tout, dummy, tcrit);
00398 }
00399
00400 Matrix
00401 DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out,
00402 const ColumnVector& tcrit)
00403 {
00404 Matrix retval;
00405
00406 int n_out = tout.capacity ();
00407 int n = size ();
00408
00409 if (n_out > 0 && n > 0)
00410 {
00411 retval.resize (n_out, n);
00412 xdot_out.resize (n_out, n);
00413
00414 for (int i = 0; i < n; i++)
00415 {
00416 retval.elem (0, i) = x.elem (i);
00417 xdot_out.elem (0, i) = xdot.elem (i);
00418 }
00419
00420 int n_crit = tcrit.capacity ();
00421
00422 if (n_crit > 0)
00423 {
00424 int i_crit = 0;
00425 int i_out = 1;
00426 double next_crit = tcrit.elem (0);
00427 double next_out;
00428 while (i_out < n_out)
00429 {
00430 bool do_restart = false;
00431
00432 next_out = tout.elem (i_out);
00433 if (i_crit < n_crit)
00434 next_crit = tcrit.elem (i_crit);
00435
00436 bool save_output;
00437 double t_out;
00438
00439 if (next_crit == next_out)
00440 {
00441 set_stop_time (next_crit);
00442 t_out = next_out;
00443 save_output = true;
00444 i_out++;
00445 i_crit++;
00446 do_restart = true;
00447 }
00448 else if (next_crit < next_out)
00449 {
00450 if (i_crit < n_crit)
00451 {
00452 set_stop_time (next_crit);
00453 t_out = next_crit;
00454 save_output = false;
00455 i_crit++;
00456 do_restart = true;
00457 }
00458 else
00459 {
00460 clear_stop_time ();
00461 t_out = next_out;
00462 save_output = true;
00463 i_out++;
00464 }
00465 }
00466 else
00467 {
00468 set_stop_time (next_crit);
00469 t_out = next_out;
00470 save_output = true;
00471 i_out++;
00472 }
00473
00474 ColumnVector x_next = do_integrate (t_out);
00475
00476 if (integration_error)
00477 return retval;
00478
00479 if (save_output)
00480 {
00481 for (int i = 0; i < n; i++)
00482 {
00483 retval.elem (i_out-1, i) = x_next.elem (i);
00484 xdot_out.elem (i_out-1, i) = xdot.elem (i);
00485 }
00486 }
00487
00488 if (do_restart)
00489 force_restart ();
00490 }
00491 }
00492 else
00493 {
00494 retval = integrate (tout, xdot_out);
00495
00496 if (integration_error)
00497 return retval;
00498 }
00499 }
00500
00501 return retval;
00502 }
00503
00504 std::string
00505 DASSL::error_message (void) const
00506 {
00507 std::string retval;
00508
00509 OSSTREAM buf;
00510 buf << t << OSSTREAM_ENDS;
00511 std::string t_curr = OSSTREAM_STR (buf);
00512 OSSTREAM_FREEZE (buf);
00513
00514 switch (istate)
00515 {
00516 case 1:
00517 retval = "a step was successfully taken in intermediate-output mode.";
00518 break;
00519
00520 case 2:
00521 retval = "integration completed by stepping exactly to TOUT";
00522 break;
00523
00524 case 3:
00525 retval = "integration to tout completed by stepping past TOUT";
00526 break;
00527
00528 case -1:
00529 retval = std::string ("a large amount of work has been expended (t =")
00530 + t_curr + ")";
00531 break;
00532
00533 case -2:
00534 retval = "the error tolerances are too stringent";
00535 break;
00536
00537 case -3:
00538 retval = std::string ("error weight became zero during problem. (t = ")
00539 + t_curr
00540 + "; solution component i vanished, and atol or atol(i) == 0)";
00541 break;
00542
00543 case -6:
00544 retval = std::string ("repeated error test failures on the last attempted step (t = ")
00545 + t_curr + ")";
00546 break;
00547
00548 case -7:
00549 retval = std::string ("the corrector could not converge (t = ")
00550 + t_curr + ")";
00551 break;
00552
00553 case -8:
00554 retval = std::string ("the matrix of partial derivatives is singular (t = ")
00555 + t_curr + ")";
00556 break;
00557
00558 case -9:
00559 retval = std::string ("the corrector could not converge (t = ")
00560 + t_curr + "; repeated test failures)";
00561 break;
00562
00563 case -10:
00564 retval = std::string ("corrector could not converge because IRES was -1 (t = ")
00565 + t_curr + ")";
00566 break;
00567
00568 case -11:
00569 retval = std::string ("return requested in user-supplied function (t = ")
00570 + t_curr + ")";
00571 break;
00572
00573 case -12:
00574 retval = "failed to compute consistent initial conditions";
00575 break;
00576
00577 case -33:
00578 retval = "unrecoverable error (see printed message)";
00579 break;
00580
00581 default:
00582 retval = "unknown error state";
00583 break;
00584 }
00585
00586 return retval;
00587 }
00588
00589
00590
00591
00592
00593