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