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 "LSODE.h"
00035 #include "f77-fcn.h"
00036 #include "lo-error.h"
00037 #include "lo-sstream.h"
00038 #include "quit.h"
00039
00040 typedef int (*lsode_fcn_ptr) (const int&, const double&, double*,
00041 double*, int&);
00042
00043 typedef int (*lsode_jac_ptr) (const int&, const double&, double*,
00044 const int&, const int&, double*, const
00045 int&);
00046
00047 extern "C"
00048 {
00049 F77_RET_T
00050 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, int&, double*, double&,
00051 double&, int&, double&, const double*, int&,
00052 int&, int&, double*, int&, int*, int&,
00053 lsode_jac_ptr, int&);
00054 }
00055
00056 static ODEFunc::ODERHSFunc user_fun;
00057 static ODEFunc::ODEJacFunc user_jac;
00058 static ColumnVector *tmp_x;
00059
00060 static int
00061 lsode_f (const int& neq, const double& time, double *,
00062 double *deriv, int& ierr)
00063 {
00064 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00065
00066 ColumnVector tmp_deriv;
00067
00068
00069
00070
00071
00072 tmp_deriv = (*user_fun) (*tmp_x, time);
00073
00074 if (tmp_deriv.length () == 0)
00075 ierr = -1;
00076 else
00077 {
00078 for (int i = 0; i < neq; i++)
00079 deriv [i] = tmp_deriv.elem (i);
00080 }
00081
00082 END_INTERRUPT_WITH_EXCEPTIONS;
00083
00084 return 0;
00085 }
00086
00087 static int
00088 lsode_j (const int& neq, const double& time, double *,
00089 const int&, const int&, double *pd, const int& nrowpd)
00090 {
00091 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00092
00093 Matrix tmp_jac (neq, neq);
00094
00095
00096
00097
00098
00099 tmp_jac = (*user_jac) (*tmp_x, time);
00100
00101 for (int j = 0; j < neq; j++)
00102 for (int i = 0; i < neq; i++)
00103 pd [nrowpd * j + i] = tmp_jac (i, j);
00104
00105 END_INTERRUPT_WITH_EXCEPTIONS;
00106
00107 return 0;
00108 }
00109
00110 ColumnVector
00111 LSODE::do_integrate (double tout)
00112 {
00113 ColumnVector retval;
00114
00115 static int nn = 0;
00116
00117 if (! initialized || restart || ODEFunc::reset || LSODE_options::reset)
00118 {
00119 integration_error = false;
00120
00121 initialized = true;
00122
00123 istate = 1;
00124
00125 int n = size ();
00126
00127 nn = n;
00128
00129 int max_maxord = 0;
00130
00131 if (integration_method () == "stiff")
00132 {
00133 max_maxord = 5;
00134
00135 if (jac)
00136 method_flag = 21;
00137 else
00138 method_flag = 22;
00139
00140 liw = 20 + n;
00141 lrw = 22 + n * (9 + n);
00142 }
00143 else
00144 {
00145 max_maxord = 12;
00146
00147 method_flag = 10;
00148
00149 liw = 20;
00150 lrw = 22 + 16 * n;
00151 }
00152
00153 maxord = maximum_order ();
00154
00155 if (maxord >= 0)
00156 {
00157 if (maxord > 0 && maxord <= max_maxord)
00158 {
00159 iwork(4) = maxord;
00160 iopt = 1;
00161 }
00162 else
00163 {
00164 (*current_liboctave_error_handler)
00165 ("lsode: invalid value for maximum order");
00166 integration_error = true;
00167 return retval;
00168 }
00169 }
00170
00171 iwork.resize (liw);
00172
00173 for (int i = 4; i < 9; i++)
00174 iwork(i) = 0;
00175
00176 rwork.resize (lrw);
00177
00178 for (int i = 4; i < 9; i++)
00179 rwork(i) = 0;
00180
00181 if (stop_time_set)
00182 {
00183 itask = 4;
00184 rwork(0) = stop_time;
00185 iopt = 1;
00186 }
00187 else
00188 {
00189 itask = 1;
00190 }
00191
00192 px = x.fortran_vec ();
00193
00194 piwork = iwork.fortran_vec ();
00195 prwork = rwork.fortran_vec ();
00196
00197 restart = false;
00198
00199
00200
00201
00202
00203
00204
00205 tmp_x = &x;
00206
00207 user_fun = function ();
00208 user_jac = jacobian_function ();
00209
00210 ColumnVector xdot = (*user_fun) (x, t);
00211
00212 if (x.length () != xdot.length ())
00213 {
00214 (*current_liboctave_error_handler)
00215 ("lsode: inconsistent sizes for state and derivative vectors");
00216
00217 integration_error = true;
00218 return retval;
00219 }
00220
00221 ODEFunc::reset = false;
00222
00223
00224
00225 rel_tol = relative_tolerance ();
00226 abs_tol = absolute_tolerance ();
00227
00228 int abs_tol_len = abs_tol.length ();
00229
00230 if (abs_tol_len == 1)
00231 itol = 1;
00232 else if (abs_tol_len == n)
00233 itol = 2;
00234 else
00235 {
00236 (*current_liboctave_error_handler)
00237 ("lsode: inconsistent sizes for state and absolute tolerance vectors");
00238
00239 integration_error = true;
00240 return retval;
00241 }
00242
00243 double iss = initial_step_size ();
00244 if (iss >= 0.0)
00245 {
00246 rwork(4) = iss;
00247 iopt = 1;
00248 }
00249
00250 double maxss = maximum_step_size ();
00251 if (maxss >= 0.0)
00252 {
00253 rwork(5) = maxss;
00254 iopt = 1;
00255 }
00256
00257 double minss = minimum_step_size ();
00258 if (minss >= 0.0)
00259 {
00260 rwork(6) = minss;
00261 iopt = 1;
00262 }
00263
00264 int sl = step_limit ();
00265 if (sl > 0)
00266 {
00267 iwork(5) = sl;
00268 iopt = 1;
00269 }
00270
00271 pabs_tol = abs_tol.fortran_vec ();
00272
00273 LSODE_options::reset = false;
00274 }
00275
00276 F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol,
00277 pabs_tol, itask, istate, iopt, prwork, lrw,
00278 piwork, liw, lsode_j, method_flag));
00279
00280 if (f77_exception_encountered)
00281 {
00282 integration_error = true;
00283 (*current_liboctave_error_handler) ("unrecoverable error in lsode");
00284 }
00285 else
00286 {
00287 switch (istate)
00288 {
00289 case 1:
00290 case 2:
00291 retval = x;
00292 t = tout;
00293 break;
00294
00295 case -1:
00296 case -2:
00297 case -3:
00298 case -4:
00299 case -5:
00300
00301 case -6:
00302
00303 case -13:
00304 integration_error = true;
00305 break;
00306
00307 default:
00308 integration_error = true;
00309 (*current_liboctave_error_handler)
00310 ("unrecognized value of istate (= %d) returned from lsode",
00311 istate);
00312 break;
00313 }
00314 }
00315
00316 return retval;
00317 }
00318
00319 std::string
00320 LSODE::error_message (void) const
00321 {
00322 std::string retval;
00323
00324 OSSTREAM buf;
00325 buf << t << OSSTREAM_ENDS;
00326 std::string t_curr = OSSTREAM_STR (buf);
00327 OSSTREAM_FREEZE (buf);
00328
00329 switch (istate)
00330 {
00331 case 1:
00332 retval = "prior to initial integration step";
00333 break;
00334
00335 case 2:
00336 retval = "successful exit";
00337 break;
00338
00339 case 3:
00340 retval = "prior to continuation call with modified parameters";
00341 break;
00342
00343 case -1:
00344 retval = std::string ("excess work on this call (t = ")
00345 + t_curr + "; perhaps wrong integration method)";
00346 break;
00347
00348 case -2:
00349 retval = "excess accuracy requested (tolerances too small)";
00350 break;
00351
00352 case -3:
00353 retval = "invalid input detected (see printed message)";
00354 break;
00355
00356 case -4:
00357 retval = std::string ("repeated error test failures (t = ")
00358 + t_curr + "check all inputs)";
00359 break;
00360
00361 case -5:
00362 retval = std::string ("repeated convergence failures (t = ")
00363 + t_curr
00364 + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)";
00365 break;
00366
00367 case -6:
00368 retval = std::string ("error weight became zero during problem. (t = ")
00369 + t_curr
00370 + "; solution component i vanished, and atol or atol(i) == 0)";
00371 break;
00372
00373 case -13:
00374 retval = "return requested in user-supplied function (t = "
00375 + t_curr + ")";
00376 break;
00377
00378 default:
00379 retval = "unknown error state";
00380 break;
00381 }
00382
00383 return retval;
00384 }
00385
00386 Matrix
00387 LSODE::do_integrate (const ColumnVector& tout)
00388 {
00389 Matrix retval;
00390
00391 int n_out = tout.capacity ();
00392 int n = size ();
00393
00394 if (n_out > 0 && n > 0)
00395 {
00396 retval.resize (n_out, n);
00397
00398 for (int i = 0; i < n; i++)
00399 retval.elem (0, i) = x.elem (i);
00400
00401 for (int j = 1; j < n_out; j++)
00402 {
00403 ColumnVector x_next = do_integrate (tout.elem (j));
00404
00405 if (integration_error)
00406 return retval;
00407
00408 for (int i = 0; i < n; i++)
00409 retval.elem (j, i) = x_next.elem (i);
00410 }
00411 }
00412
00413 return retval;
00414 }
00415
00416 Matrix
00417 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
00418 {
00419 Matrix retval;
00420
00421 int n_out = tout.capacity ();
00422 int n = size ();
00423
00424 if (n_out > 0 && n > 0)
00425 {
00426 retval.resize (n_out, n);
00427
00428 for (int i = 0; i < n; i++)
00429 retval.elem (0, i) = x.elem (i);
00430
00431 int n_crit = tcrit.capacity ();
00432
00433 if (n_crit > 0)
00434 {
00435 int i_crit = 0;
00436 int i_out = 1;
00437 double next_crit = tcrit.elem (0);
00438 double next_out;
00439 while (i_out < n_out)
00440 {
00441 bool do_restart = false;
00442
00443 next_out = tout.elem (i_out);
00444 if (i_crit < n_crit)
00445 next_crit = tcrit.elem (i_crit);
00446
00447 int save_output;
00448 double t_out;
00449
00450 if (next_crit == next_out)
00451 {
00452 set_stop_time (next_crit);
00453 t_out = next_out;
00454 save_output = 1;
00455 i_out++;
00456 i_crit++;
00457 do_restart = true;
00458 }
00459 else if (next_crit < next_out)
00460 {
00461 if (i_crit < n_crit)
00462 {
00463 set_stop_time (next_crit);
00464 t_out = next_crit;
00465 save_output = 0;
00466 i_crit++;
00467 do_restart = true;
00468 }
00469 else
00470 {
00471 clear_stop_time ();
00472 t_out = next_out;
00473 save_output = 1;
00474 i_out++;
00475 }
00476 }
00477 else
00478 {
00479 set_stop_time (next_crit);
00480 t_out = next_out;
00481 save_output = 1;
00482 i_out++;
00483 }
00484
00485 ColumnVector x_next = do_integrate (t_out);
00486
00487 if (integration_error)
00488 return retval;
00489
00490 if (save_output)
00491 {
00492 for (int i = 0; i < n; i++)
00493 retval.elem (i_out-1, i) = x_next.elem (i);
00494 }
00495
00496 if (do_restart)
00497 force_restart ();
00498 }
00499 }
00500 else
00501 {
00502 retval = do_integrate (tout);
00503
00504 if (integration_error)
00505 return retval;
00506 }
00507 }
00508
00509 return retval;
00510 }
00511
00512
00513
00514
00515
00516