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
00035 #include "Array.h"
00036 #include "Array.cc"
00037
00038 #include "ODESSA.h"
00039 #include "f77-fcn.h"
00040 #include "lo-error.h"
00041 #include "lo-sstream.h"
00042 #include "quit.h"
00043
00044 typedef int (*odessa_fcn_ptr) (int*, const double&, double*,
00045 double*, double*);
00046
00047 typedef int (*odessa_jac_ptr) (int*, const double&, double*,
00048 double*, const int&, const int&,
00049 double*, const int&);
00050
00051 typedef int (*odessa_dfdp_ptr) (int*, const double&, double*,
00052 double*, double*, const int&);
00053
00054
00055 extern "C"
00056 {
00057 F77_RET_T
00058 F77_FUNC (dodessa, DODESSA) (odessa_fcn_ptr, odessa_dfdp_ptr, int*,
00059 double*, double*, double&, double&,
00060 int&, double&, const double*, int&,
00061 int&, int*, double*, int&, int*, int&,
00062 odessa_jac_ptr, int&);
00063 }
00064
00065 INSTANTIATE_ARRAY (Matrix);
00066
00067 static ODESFunc::ODES_fsub user_fsub;
00068 static ODESFunc::ODES_bsub user_bsub;
00069 static ODESFunc::ODES_jsub user_jsub;
00070
00071
00072 static int
00073 odessa_f (int* neq, const double& t, double *state,
00074 double *par, double *fval)
00075 {
00076 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00077
00078 int n = neq[0];
00079 int n_par = neq[1];
00080
00081
00082
00083 ColumnVector tmp_state (n);
00084 for (int i = 0; i < n; i++)
00085 tmp_state(i) = state[i];
00086
00087 ColumnVector tmp_param (n_par);
00088 for (int i = 0; i < n_par; i++)
00089 tmp_param(i) = par[i];
00090
00091 ColumnVector tmp_fval = user_fsub (tmp_state, t, tmp_param);
00092
00093 if (tmp_fval.length () == 0)
00094 octave_jump_to_enclosing_context ();
00095 else
00096 {
00097 for (int i = 0; i < n; i++)
00098 fval[i] = tmp_fval(i);
00099 }
00100
00101 END_INTERRUPT_WITH_EXCEPTIONS;
00102
00103 return 0;
00104 }
00105
00106 static int
00107 odessa_j (int* neq, const double& t, double *state,
00108 double *par, const int& , const int& ,
00109 double *pd, const int& nrowpd)
00110 {
00111 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00112
00113 int n = neq[0];
00114 int n_par = neq[1];
00115
00116
00117 ColumnVector tmp_state (n);
00118 for (int i = 0; i < n; i++)
00119 tmp_state(i) = state[i];
00120
00121 ColumnVector tmp_param (n_par);
00122 for (int i = 0; i < n_par; i++)
00123 tmp_param(i) = par[i];
00124
00125 Matrix tmp_fval = user_jsub (tmp_state, t, tmp_param);
00126
00127 if (tmp_fval.length () == 0)
00128 octave_jump_to_enclosing_context ();
00129 else
00130 {
00131 for (int j = 0; j < n; j++)
00132 for (int i = 0; i < nrowpd; i++)
00133 pd[nrowpd*j+i] = tmp_fval(i,j);
00134 }
00135
00136 END_INTERRUPT_WITH_EXCEPTIONS;
00137
00138 return 0;
00139 }
00140
00141 static int
00142 odessa_b (int* neq, const double& t, double *state,
00143 double *par, double *dfdp, const int& jpar)
00144
00145 {
00146 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00147
00148 int n = neq[0];
00149 int n_par = neq[1];
00150
00151
00152 ColumnVector tmp_state (n);
00153 for (int i = 0; i < n; i++)
00154 tmp_state(i) = state[i];
00155
00156 ColumnVector tmp_param (n_par);
00157 for (int i = 0; i < n_par; i++)
00158 tmp_param(i) = par[i];
00159
00160 ColumnVector tmp_fval = user_bsub (tmp_state, t, tmp_param, jpar);
00161
00162 if (tmp_fval.length () == 0)
00163 octave_jump_to_enclosing_context ();
00164 else
00165 {
00166 for (int i = 0; i < n; i++)
00167 dfdp[i] = tmp_fval(i);
00168 }
00169
00170 END_INTERRUPT_WITH_EXCEPTIONS;
00171
00172 return 0;
00173 }
00174
00175 ODESSA::ODESSA (void) : ODES (), ODESSA_options ()
00176 {
00177 initialized = false;
00178
00179 neq.resize(2);
00180 n = size ();
00181
00182 iopt.resize(4);
00183
00184 itask = 1;
00185 iopt(0) = 0;
00186 isopt = 0;
00187 iopt(1) = isopt;
00188 npar = 0;
00189 neq(0) = n;
00190 neq(1) = npar;
00191
00192 sanity_checked = false;
00193 }
00194
00195 ODESSA::ODESSA (const ColumnVector& s, double tm, ODESFunc& f)
00196 : ODES (s, tm, f), ODESSA_options ()
00197 {
00198 initialized = false;
00199
00200 neq.resize(2);
00201 n = size ();
00202
00203 iopt.resize(4);
00204 itask = 1;
00205 iopt(0) = 0;
00206 isopt = 0;
00207 iopt(1) = isopt;
00208
00209 sanity_checked = false;
00210
00211 npar = 0;
00212 neq(0) = n;
00213 neq(1) = npar;
00214
00215 y.resize (n, 1, 0.0);
00216 }
00217
00218 ODESSA::ODESSA (const ColumnVector& s, const ColumnVector& xtheta,
00219 const Matrix& sensitivity_guess, double tm, ODESFunc& f)
00220 : ODES (s, xtheta, tm, f)
00221 {
00222 initialized = false;
00223
00224 neq.resize(2);
00225 n = s.length();
00226 npar = xtheta.length();
00227
00228 neq(0) = n;
00229 neq(1) = npar;
00230
00231 sx0 = sensitivity_guess;
00232 par.resize (npar);
00233
00234 for (int i = 0; i < npar; i++)
00235 {
00236 par(i) = xtheta(i);
00237 }
00238
00239 sanity_checked = false;
00240
00241 npar = xtheta.length ();
00242
00243 iopt.resize(4);
00244 itask = 1;
00245 iopt(0) = 0;
00246 isopt = 1;
00247 iopt(1) = isopt;
00248
00249 y.resize (n, npar+1, 0.0);
00250 }
00251
00252 void
00253 ODESSA::integrate (double tout)
00254 {
00255 ODESSA_result retval;
00256
00257 if (! initialized)
00258 {
00259
00260 for (int i = 0; i < n; i++)
00261 y(i,0) = x(i);
00262
00263 if (npar > 0)
00264 {
00265 for (int j = 0; j < npar; j++)
00266 for (int i = 0; i < n; i++)
00267 y(i,j+1) = sx0(i,j);
00268 }
00269
00270 integration_error = false;
00271
00272 user_fsub = ODESFunc::fsub_function ();
00273 user_bsub = ODESFunc::bsub_function ();
00274 user_jsub = ODESFunc::jsub_function ();
00275
00276 int idf;
00277
00278 if (user_bsub)
00279 idf = 1;
00280 else
00281 idf = 0;
00282
00283 iopt(2) = idf;
00284
00285 if (restart)
00286 {
00287 restart = false;
00288 istate = 1;
00289 }
00290
00291 int max_maxord = 0;
00292
00293 if (integration_method () == "stiff")
00294 {
00295 if (user_jsub)
00296 method_flag = 21;
00297 else
00298 method_flag = 22;
00299
00300 max_maxord = 5;
00301
00302 if (isopt)
00303 {
00304 liw = 21 + n + npar;
00305 lrw = 22 + 8*(npar+1)*n + n*n + n;
00306 }
00307 else
00308 {
00309 liw = 20 + n;
00310 lrw = 22 + 9*n + n*n;
00311 }
00312 }
00313 else
00314 {
00315 max_maxord = 12;
00316
00317 if (isopt)
00318 {
00319 if (user_jsub)
00320 method_flag = 11;
00321 else
00322 method_flag = 12;
00323 liw = 21 + n + npar;
00324 lrw = 22 + 15*(npar+1)*n + n*n + n;
00325 }
00326 else
00327 {
00328 method_flag = 10;
00329 liw = 20 + n;
00330 lrw = 22 + 16 * n;
00331 }
00332 }
00333
00334 if (iwork.length () != liw)
00335 {
00336 iwork.resize (liw);
00337
00338 for (int i = 4; i < 10; i++)
00339 iwork.elem (i) = 0;
00340 }
00341
00342 if (rwork.length () != lrw)
00343 {
00344 rwork.resize (lrw);
00345
00346 for (int i = 4; i < 10; i++)
00347 rwork.elem (i) = 0.0;
00348 }
00349
00350 maxord = maximum_order ();
00351
00352 if (maxord >= 0)
00353 {
00354 if (maxord > 0 && maxord <= max_maxord)
00355 {
00356 iwork(4) = maxord;
00357 iopt(0) = 1;
00358 }
00359 else
00360 {
00361 (*current_liboctave_error_handler)
00362 ("odessa: invalid value for maximum order");
00363 integration_error = true;
00364 return;
00365 }
00366 }
00367
00368 initialized = true;
00369 }
00370
00371 integration_error = false;
00372
00373
00374
00375
00376
00377
00378 if (! sanity_checked)
00379 {
00380 ColumnVector fval = user_fsub (x, t, theta);
00381
00382 if (fval.length () != x.length ())
00383 {
00384 (*current_liboctave_error_handler)
00385 ("odessa: inconsistent sizes for state and residual vectors");
00386
00387 integration_error = true;
00388 return;
00389 }
00390
00391 sanity_checked = true;
00392 }
00393
00394 if (stop_time_set)
00395 {
00396 itask = 4;
00397 rwork.elem (0) = stop_time;
00398 iopt(0) = 1;
00399 }
00400 else
00401 {
00402 itask = 1;
00403 }
00404
00405 double rel_tol = relative_tolerance ();
00406
00407 Array<double> abs_tol = absolute_tolerance ();
00408
00409
00410 int abs_tol_len = abs_tol.length ();
00411
00412 int itol;
00413
00414 if (abs_tol_len == 1)
00415 itol = 1;
00416 else if (abs_tol_len == n)
00417 itol = 2;
00418 else
00419 {
00420 (*current_liboctave_error_handler)
00421 ("odessa: inconsistent sizes for state and absolute tolerance vectors");
00422
00423 integration_error = 1;
00424 return;
00425 }
00426
00427 if (initial_step_size () >= 0.0)
00428 {
00429 rwork.elem (4) = initial_step_size ();
00430 iopt(0) = 1;
00431 }
00432
00433 if (maximum_step_size () >= 0.0)
00434 {
00435 rwork.elem (5) = maximum_step_size ();
00436 iopt(0) = 1;
00437 }
00438
00439 if (minimum_step_size () >= 0.0)
00440 {
00441 rwork.elem (6) = minimum_step_size ();
00442 iopt(0) = 1;
00443 }
00444
00445 if (step_limit () > 0)
00446 {
00447 iwork.elem (5) = step_limit ();
00448 iopt(0) = 1;
00449 }
00450
00451 py = y.fortran_vec ();
00452 piwork = iwork.fortran_vec ();
00453 prwork = rwork.fortran_vec ();
00454 ppar = par.fortran_vec ();
00455 piopt = iopt.fortran_vec ();
00456 pneq = neq.fortran_vec ();
00457
00458 const double *pabs_tol = abs_tol.fortran_vec ();
00459
00460 F77_XFCN (dodessa, DODESSA, (odessa_f, odessa_b, pneq, py, ppar, t,
00461 tout, itol, rel_tol, pabs_tol, itask,
00462 istate, piopt, prwork, lrw, piwork, liw,
00463 odessa_j, method_flag));
00464
00465 if (f77_exception_encountered)
00466 {
00467 integration_error = true;
00468 (*current_liboctave_error_handler) ("unrecoverable error in odessa");
00469 }
00470 else
00471 {
00472 switch (istate)
00473 {
00474 case 1:
00475 case 2:
00476 t = tout;
00477 break;
00478
00479 case -1:
00480 case -2:
00481 case -3:
00482 case -4:
00483 case -5:
00484
00485 case -6:
00486
00487 case -13:
00488 integration_error = 1;
00489 break;
00490
00491 default:
00492 integration_error = 1;
00493 (*current_liboctave_error_handler)
00494 ("unrecognized value of istate (= %d) returned from odessa",
00495 istate);
00496 break;
00497 }
00498 }
00499 }
00500
00501 std::string
00502 ODESSA::error_message (void) const
00503 {
00504 std::string retval;
00505
00506 OSSTREAM buf;
00507 buf << t << OSSTREAM_ENDS;
00508 std::string t_curr = OSSTREAM_STR (buf);
00509 OSSTREAM_FREEZE (buf);
00510
00511 switch (istate)
00512 {
00513 case 1:
00514 retval = "prior to initial integration step";
00515 break;
00516
00517 case 2:
00518 retval = "successful exit";
00519 break;
00520
00521 case 3:
00522 retval = "prior to continuation call with modified parameters";
00523 break;
00524
00525 case -1:
00526 retval = std::string ("excess work on this call (t = ")
00527 + t_curr + "; perhaps wrong integration method)";
00528 break;
00529
00530 case -2:
00531 retval = "excess accuracy requested (tolerances too small)";
00532 break;
00533
00534 case -3:
00535 retval = "invalid input detected (see printed message)";
00536 break;
00537
00538 case -4:
00539 retval = std::string ("repeated error test failures (t = ")
00540 + t_curr + "check all inputs)";
00541 break;
00542
00543 case -5:
00544 retval = std::string ("repeated convergence failures (t = ")
00545 + t_curr
00546 + "perhaps bad jacobian supplied or wrong choice of integration method or tolerances)";
00547 break;
00548
00549 case -6:
00550 retval = std::string ("error weight became zero during problem. (t = ")
00551 + t_curr
00552 + "; solution component i vanished, and atol or atol(i) == 0)";
00553 break;
00554
00555 case -13:
00556 retval = "return requested in user-supplied function (t = "
00557 + t_curr + ")";
00558 break;
00559
00560 default:
00561 retval = "unknown error state";
00562 break;
00563 }
00564
00565 return retval;
00566 }
00567
00568
00569 ODESSA_result
00570 ODESSA::integrate (const ColumnVector& tout)
00571 {
00572 ODESSA_result retval;
00573
00574 Matrix x_out;
00575
00576 Array<Matrix> x_s_out;
00577
00578 int n_out = tout.capacity ();
00579
00580 if (n_out > 0 && n > 0)
00581 {
00582 x_out.resize (n_out, n);
00583
00584 x_s_out.resize (npar, Matrix (n_out, n, 0.0));
00585
00586 for (int j = 0; j < n_out; j++)
00587 {
00588 integrate (tout(j));
00589
00590 if (integration_error)
00591 {
00592 retval = ODESSA_result (x_out, x_s_out);
00593 return retval;
00594 }
00595
00596 for (int i = 0; i < n; i++)
00597 {
00598 x_out(j,i) = y(i,0);
00599
00600 for (int k = 0; k < npar; k++)
00601 {
00602 x_s_out(k)(j,i) = y(i,k+1);
00603 }
00604 }
00605 }
00606 }
00607
00608 retval = ODESSA_result (x_out, x_s_out);
00609
00610 return retval;
00611 }
00612
00613 ODESSA_result
00614 ODESSA::integrate (const ColumnVector& tout, const ColumnVector& tcrit)
00615 {
00616 ODESSA_result retval;
00617
00618 Matrix x_out;
00619
00620 Array<Matrix> x_s_out;
00621
00622 int n_out = tout.capacity ();
00623
00624 if (n_out > 0 && n > 0)
00625 {
00626 x_out.resize (n_out, n);
00627
00628 x_s_out.resize (npar, Matrix (n_out, n, 0.0));
00629
00630 int n_crit = tcrit.capacity ();
00631
00632 if (n_crit > 0)
00633 {
00634 int i_crit = 0;
00635 int i_out = 0;
00636 double next_crit = tcrit(0);
00637 double next_out;
00638 while (i_out < n_out)
00639 {
00640 bool do_restart = false;
00641
00642 next_out = tout(i_out);
00643 if (i_crit < n_crit)
00644 next_crit = tcrit(i_crit);
00645
00646 int save_output;
00647 double t_out;
00648
00649 if (next_crit == next_out)
00650 {
00651 set_stop_time (next_crit);
00652 t_out = next_out;
00653 save_output = true;
00654 i_out++;
00655 i_crit++;
00656 do_restart = true;
00657 }
00658 else if (next_crit < next_out)
00659 {
00660 if (i_crit < n_crit)
00661 {
00662 set_stop_time (next_crit);
00663 t_out = next_crit;
00664 save_output = false;
00665 i_crit++;
00666 do_restart = true;
00667 }
00668 else
00669 {
00670 clear_stop_time ();
00671 t_out = next_out;
00672 save_output = true;
00673 i_out++;
00674 }
00675 }
00676 else
00677 {
00678 set_stop_time (next_crit);
00679 t_out = next_out;
00680 save_output = true;
00681 i_out++;
00682 }
00683 integrate (t_out);
00684 if (integration_error)
00685 {
00686 retval = ODESSA_result (x_out, x_s_out);
00687 return retval;
00688 }
00689 if (save_output)
00690 {
00691 for (int i = 0; i < n; i++)
00692 {
00693 x_out(i_out-1,i) = y(i,0);
00694
00695 for (int k = 0; k < npar; k++)
00696 {
00697 x_s_out(k)(i_out-1,i) = y(i,k+1);
00698 }
00699 }
00700 }
00701
00702 if (do_restart)
00703 force_restart ();
00704 }
00705
00706 retval = ODESSA_result (x_out, x_s_out);
00707 }
00708 else
00709 {
00710 retval = integrate (tout);
00711
00712 if (integration_error)
00713 return retval;
00714 }
00715 }
00716
00717 return retval;
00718 }
00719
00720
00721
00722
00723
00724