Mercurial > octave-nkf
annotate liboctave/DASPK.cc @ 7948:af10baa63915 ss-3-1-50
3.1.50 snapshot
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Fri, 18 Jul 2008 17:42:48 -0400 |
| parents | 29980c6b8604 |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 3912 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1996, 1997, 2002, 2003, 2004, 2005, 2006, 2007 |
| 4 John W. Eaton | |
| 3912 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 3912 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 3912 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <cfloat> | |
| 29 | |
| 5765 | 30 #include <sstream> |
| 31 | |
| 3912 | 32 #include "DASPK.h" |
| 33 #include "f77-fcn.h" | |
| 34 #include "lo-error.h" | |
| 7231 | 35 #include "lo-math.h" |
| 4155 | 36 #include "quit.h" |
| 3912 | 37 |
| 5275 | 38 typedef octave_idx_type (*daspk_fcn_ptr) (const double&, const double*, |
| 3912 | 39 const double*, const double&, |
| 5275 | 40 double*, octave_idx_type&, double*, octave_idx_type*); |
| 3912 | 41 |
| 5275 | 42 typedef octave_idx_type (*daspk_jac_ptr) (const double&, const double*, |
| 3912 | 43 const double*, double*, |
| 5275 | 44 const double&, double*, octave_idx_type*); |
| 3912 | 45 |
| 5275 | 46 typedef octave_idx_type (*daspk_psol_ptr) (const octave_idx_type&, const double&, |
| 3912 | 47 const double*, const double*, |
| 48 const double*, const double&, | |
| 5275 | 49 const double*, double*, octave_idx_type*, |
| 50 double*, const double&, octave_idx_type&, | |
| 51 double*, octave_idx_type*); | |
| 3912 | 52 |
| 53 extern "C" | |
| 4552 | 54 { |
| 55 F77_RET_T | |
| 5275 | 56 F77_FUNC (ddaspk, DDASPK) (daspk_fcn_ptr, const octave_idx_type&, double&, |
| 57 double*, double*, double&, const octave_idx_type*, | |
| 58 const double*, const double*, octave_idx_type&, | |
| 59 double*, const octave_idx_type&, octave_idx_type*, const octave_idx_type&, | |
| 60 const double*, const octave_idx_type*, | |
| 4552 | 61 daspk_jac_ptr, daspk_psol_ptr); |
| 62 } | |
| 3912 | 63 |
| 64 static DAEFunc::DAERHSFunc user_fun; | |
| 65 static DAEFunc::DAEJacFunc user_jac; | |
| 5275 | 66 static octave_idx_type nn; |
| 3912 | 67 |
| 5275 | 68 static octave_idx_type |
| 3912 | 69 ddaspk_f (const double& time, const double *state, const double *deriv, |
| 5275 | 70 const double&, double *delta, octave_idx_type& ires, double *, octave_idx_type *) |
| 3912 | 71 { |
| 4180 | 72 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 73 | |
| 3912 | 74 ColumnVector tmp_deriv (nn); |
| 75 ColumnVector tmp_state (nn); | |
| 76 ColumnVector tmp_delta (nn); | |
| 77 | |
| 5275 | 78 for (octave_idx_type i = 0; i < nn; i++) |
| 3912 | 79 { |
| 80 tmp_deriv.elem (i) = deriv [i]; | |
| 81 tmp_state.elem (i) = state [i]; | |
| 82 } | |
| 83 | |
| 84 tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires); | |
| 85 | |
| 86 if (ires >= 0) | |
| 87 { | |
| 88 if (tmp_delta.length () == 0) | |
| 89 ires = -2; | |
| 90 else | |
| 91 { | |
| 5275 | 92 for (octave_idx_type i = 0; i < nn; i++) |
| 3912 | 93 delta [i] = tmp_delta.elem (i); |
| 94 } | |
| 95 } | |
| 96 | |
| 4180 | 97 END_INTERRUPT_WITH_EXCEPTIONS; |
| 98 | |
| 3912 | 99 return 0; |
| 100 } | |
| 101 | |
| 102 //NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT, | |
| 103 //C WP, IWP, B, EPLIN, IER, RPAR, IPAR) | |
| 104 | |
| 5275 | 105 static octave_idx_type |
| 106 ddaspk_psol (const octave_idx_type&, const double&, const double *, | |
| 4662 | 107 const double *, const double *, const double&, |
| 5275 | 108 const double *, double *, octave_idx_type *, double *, |
| 109 const double&, octave_idx_type&, double *, octave_idx_type*) | |
| 3912 | 110 { |
| 4180 | 111 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 112 | |
| 3912 | 113 abort (); |
| 4180 | 114 |
| 115 END_INTERRUPT_WITH_EXCEPTIONS; | |
| 116 | |
| 3946 | 117 return 0; |
| 3912 | 118 } |
| 119 | |
| 3991 | 120 |
| 5275 | 121 static octave_idx_type |
| 3991 | 122 ddaspk_j (const double& time, const double *state, const double *deriv, |
| 5275 | 123 double *pd, const double& cj, double *, octave_idx_type *) |
| 3912 | 124 { |
| 4180 | 125 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 126 | |
| 5775 | 127 // FIXME -- would be nice to avoid copying the data. |
| 3991 | 128 |
| 3912 | 129 ColumnVector tmp_state (nn); |
| 130 ColumnVector tmp_deriv (nn); | |
| 131 | |
| 5275 | 132 for (octave_idx_type i = 0; i < nn; i++) |
| 3991 | 133 { |
| 134 tmp_deriv.elem (i) = deriv [i]; | |
| 135 tmp_state.elem (i) = state [i]; | |
| 136 } | |
| 3912 | 137 |
| 3991 | 138 Matrix tmp_pd = user_jac (tmp_state, tmp_deriv, time, cj); |
| 3912 | 139 |
| 5275 | 140 for (octave_idx_type j = 0; j < nn; j++) |
| 141 for (octave_idx_type i = 0; i < nn; i++) | |
| 3991 | 142 pd [nn * j + i] = tmp_pd.elem (i, j); |
| 3912 | 143 |
| 4180 | 144 END_INTERRUPT_WITH_EXCEPTIONS; |
| 145 | |
| 3912 | 146 return 0; |
| 147 } | |
| 148 | |
| 149 ColumnVector | |
| 150 DASPK::do_integrate (double tout) | |
| 151 { | |
| 5775 | 152 // FIXME -- should handle all this option stuff just once |
| 4047 | 153 // for each new problem. |
| 154 | |
| 3912 | 155 ColumnVector retval; |
| 156 | |
| 4049 | 157 if (! initialized || restart || DAEFunc::reset|| DASPK_options::reset) |
| 3912 | 158 { |
| 4049 | 159 integration_error = false; |
| 160 | |
| 161 initialized = true; | |
| 162 | |
| 163 info.resize (20); | |
| 3912 | 164 |
| 5275 | 165 for (octave_idx_type i = 0; i < 20; i++) |
| 4049 | 166 info(i) = 0; |
| 167 | |
| 168 pinfo = info.fortran_vec (); | |
| 3912 | 169 |
| 5275 | 170 octave_idx_type n = size (); |
| 3912 | 171 |
| 4049 | 172 nn = n; |
| 173 | |
| 174 info(0) = 0; | |
| 3912 | 175 |
| 4049 | 176 if (stop_time_set) |
| 177 { | |
| 178 rwork(0) = stop_time; | |
| 179 info(3) = 1; | |
| 180 } | |
| 181 else | |
| 182 info(3) = 0; | |
| 3912 | 183 |
| 4049 | 184 px = x.fortran_vec (); |
| 185 pxdot = xdot.fortran_vec (); | |
| 186 | |
| 187 // DAEFunc | |
| 188 | |
| 189 user_fun = DAEFunc::function (); | |
| 190 user_jac = DAEFunc::jacobian_function (); | |
| 3912 | 191 |
| 4049 | 192 if (user_fun) |
| 193 { | |
| 5275 | 194 octave_idx_type ires = 0; |
| 3912 | 195 |
| 4049 | 196 ColumnVector res = (*user_fun) (x, xdot, t, ires); |
| 3912 | 197 |
| 4049 | 198 if (res.length () != x.length ()) |
| 199 { | |
| 200 (*current_liboctave_error_handler) | |
| 201 ("daspk: inconsistent sizes for state and residual vectors"); | |
| 3912 | 202 |
| 4049 | 203 integration_error = true; |
| 204 return retval; | |
| 205 } | |
| 206 } | |
| 207 else | |
| 3912 | 208 { |
| 209 (*current_liboctave_error_handler) | |
| 4049 | 210 ("daspk: no user supplied RHS subroutine!"); |
| 211 | |
| 212 integration_error = true; | |
| 213 return retval; | |
| 214 } | |
| 215 | |
| 216 info(4) = user_jac ? 1 : 0; | |
| 217 | |
| 218 DAEFunc::reset = false; | |
| 219 | |
| 5275 | 220 octave_idx_type eiq = enforce_inequality_constraints (); |
| 221 octave_idx_type ccic = compute_consistent_initial_condition (); | |
| 222 octave_idx_type eavfet = exclude_algebraic_variables_from_error_test (); | |
| 4144 | 223 |
| 224 liw = 40 + n; | |
| 225 if (eiq == 1 || eiq == 3) | |
| 226 liw += n; | |
| 227 if (ccic == 1 || eavfet == 1) | |
| 228 liw += n; | |
| 229 | |
| 4842 | 230 lrw = 50 + 9*n + n*n; |
| 4144 | 231 if (eavfet == 1) |
| 232 lrw += n; | |
| 233 | |
| 234 iwork.resize (liw); | |
| 235 rwork.resize (lrw); | |
| 236 | |
| 237 piwork = iwork.fortran_vec (); | |
| 238 prwork = rwork.fortran_vec (); | |
| 239 | |
| 4049 | 240 // DASPK_options |
| 241 | |
| 4144 | 242 abs_tol = absolute_tolerance (); |
| 243 rel_tol = relative_tolerance (); | |
| 4049 | 244 |
| 5275 | 245 octave_idx_type abs_tol_len = abs_tol.length (); |
| 246 octave_idx_type rel_tol_len = rel_tol.length (); | |
| 4049 | 247 |
| 248 if (abs_tol_len == 1 && rel_tol_len == 1) | |
| 249 { | |
| 250 info(1) = 0; | |
| 251 } | |
| 252 else if (abs_tol_len == n && rel_tol_len == n) | |
| 253 { | |
| 254 info(1) = 1; | |
| 255 } | |
| 256 else | |
| 257 { | |
| 258 (*current_liboctave_error_handler) | |
| 259 ("daspk: inconsistent sizes for tolerance arrays"); | |
| 3912 | 260 |
| 3995 | 261 integration_error = true; |
| 3912 | 262 return retval; |
| 263 } | |
| 264 | |
| 4144 | 265 pabs_tol = abs_tol.fortran_vec (); |
| 266 prel_tol = rel_tol.fortran_vec (); | |
| 267 | |
| 4049 | 268 double hmax = maximum_step_size (); |
| 269 if (hmax >= 0.0) | |
| 270 { | |
| 271 rwork(1) = hmax; | |
| 272 info(6) = 1; | |
| 273 } | |
| 274 else | |
| 275 info(6) = 0; | |
| 3998 | 276 |
| 4049 | 277 double h0 = initial_step_size (); |
| 278 if (h0 >= 0.0) | |
| 4047 | 279 { |
| 4049 | 280 rwork(2) = h0; |
| 281 info(7) = 1; | |
| 4047 | 282 } |
| 283 else | |
| 4049 | 284 info(7) = 0; |
| 285 | |
| 5275 | 286 octave_idx_type maxord = maximum_order (); |
| 4049 | 287 if (maxord >= 0) |
| 4047 | 288 { |
| 4049 | 289 if (maxord > 0 && maxord < 6) |
| 290 { | |
| 291 info(8) = 1; | |
| 292 iwork(2) = maxord; | |
| 293 } | |
| 294 else | |
| 295 { | |
| 296 (*current_liboctave_error_handler) | |
| 297 ("daspk: invalid value for maximum order"); | |
| 298 integration_error = true; | |
| 299 return retval; | |
| 300 } | |
| 301 } | |
| 302 | |
| 303 switch (eiq) | |
| 304 { | |
| 305 case 1: | |
| 306 case 3: | |
| 307 { | |
| 5275 | 308 Array<octave_idx_type> ict = inequality_constraint_types (); |
| 4049 | 309 |
| 310 if (ict.length () == n) | |
| 311 { | |
| 5275 | 312 for (octave_idx_type i = 0; i < n; i++) |
| 4049 | 313 { |
| 5275 | 314 octave_idx_type val = ict(i); |
| 4049 | 315 if (val < -2 || val > 2) |
| 316 { | |
| 317 (*current_liboctave_error_handler) | |
| 318 ("daspk: invalid value for inequality constraint type"); | |
| 319 integration_error = true; | |
| 320 return retval; | |
| 321 } | |
| 322 iwork(40+i) = val; | |
| 323 } | |
| 324 } | |
| 325 else | |
| 326 { | |
| 327 (*current_liboctave_error_handler) | |
| 328 ("daspk: inequality constraint types size mismatch"); | |
| 329 integration_error = true; | |
| 330 return retval; | |
| 331 } | |
| 332 } | |
| 333 // Fall through... | |
| 334 | |
| 4144 | 335 case 0: |
| 4049 | 336 case 2: |
| 337 info(9) = eiq; | |
| 338 break; | |
| 339 | |
| 340 default: | |
| 4047 | 341 (*current_liboctave_error_handler) |
| 4049 | 342 ("daspk: invalid value for enforce inequality constraints option"); |
| 4047 | 343 integration_error = true; |
| 344 return retval; | |
| 345 } | |
| 4049 | 346 |
| 347 if (ccic) | |
| 348 { | |
| 349 if (ccic == 1) | |
| 350 { | |
| 5775 | 351 // FIXME -- this code is duplicated below. |
| 4049 | 352 |
| 5275 | 353 Array<octave_idx_type> av = algebraic_variables (); |
| 4047 | 354 |
| 4049 | 355 if (av.length () == n) |
| 356 { | |
| 5275 | 357 octave_idx_type lid; |
| 4049 | 358 if (eiq == 0 || eiq == 2) |
| 359 lid = 40; | |
| 360 else if (eiq == 1 || eiq == 3) | |
| 361 lid = 40 + n; | |
| 362 else | |
| 363 abort (); | |
| 4047 | 364 |
| 5275 | 365 for (octave_idx_type i = 0; i < n; i++) |
| 4049 | 366 iwork(lid+i) = av(i) ? -1 : 1; |
| 367 } | |
| 368 else | |
| 369 { | |
| 370 (*current_liboctave_error_handler) | |
| 371 ("daspk: algebraic variables size mismatch"); | |
| 372 integration_error = true; | |
| 373 return retval; | |
| 374 } | |
| 375 } | |
| 376 else if (ccic != 2) | |
| 377 { | |
| 378 (*current_liboctave_error_handler) | |
| 379 ("daspk: invalid value for compute consistent initial condition option"); | |
| 380 integration_error = true; | |
| 381 return retval; | |
| 382 } | |
| 4047 | 383 |
| 4049 | 384 info(10) = ccic; |
| 385 } | |
| 4047 | 386 |
| 4049 | 387 if (eavfet) |
| 388 { | |
| 389 info(15) = 1; | |
| 4047 | 390 |
| 5775 | 391 // FIXME -- this code is duplicated above. |
| 4047 | 392 |
| 5275 | 393 Array<octave_idx_type> av = algebraic_variables (); |
| 4047 | 394 |
| 395 if (av.length () == n) | |
| 396 { | |
| 5275 | 397 octave_idx_type lid; |
| 4047 | 398 if (eiq == 0 || eiq == 2) |
| 399 lid = 40; | |
| 400 else if (eiq == 1 || eiq == 3) | |
| 401 lid = 40 + n; | |
| 402 else | |
| 403 abort (); | |
| 404 | |
| 5275 | 405 for (octave_idx_type i = 0; i < n; i++) |
| 4047 | 406 iwork(lid+i) = av(i) ? -1 : 1; |
| 407 } | |
| 4049 | 408 } |
| 409 | |
| 410 if (use_initial_condition_heuristics ()) | |
| 411 { | |
| 412 Array<double> ich = initial_condition_heuristics (); | |
| 413 | |
| 414 if (ich.length () == 6) | |
| 415 { | |
| 5275 | 416 iwork(31) = NINTbig (ich(0)); |
| 417 iwork(32) = NINTbig (ich(1)); | |
| 418 iwork(33) = NINTbig (ich(2)); | |
| 419 iwork(34) = NINTbig (ich(3)); | |
| 4049 | 420 |
| 421 rwork(13) = ich(4); | |
| 422 rwork(14) = ich(5); | |
| 423 } | |
| 4047 | 424 else |
| 425 { | |
| 426 (*current_liboctave_error_handler) | |
| 4049 | 427 ("daspk: invalid initial condition heuristics option"); |
| 4047 | 428 integration_error = true; |
| 429 return retval; | |
| 430 } | |
| 4049 | 431 |
| 432 info(16) = 1; | |
| 4047 | 433 } |
| 4049 | 434 |
| 5275 | 435 octave_idx_type pici = print_initial_condition_info (); |
| 4049 | 436 switch (pici) |
| 4047 | 437 { |
| 4049 | 438 case 0: |
| 439 case 1: | |
| 440 case 2: | |
| 441 info(17) = pici; | |
| 442 break; | |
| 443 | |
| 444 default: | |
| 4047 | 445 (*current_liboctave_error_handler) |
| 4049 | 446 ("daspk: invalid value for print initial condition info option"); |
| 4047 | 447 integration_error = true; |
| 448 return retval; | |
| 4049 | 449 break; |
| 4047 | 450 } |
| 451 | |
| 4049 | 452 DASPK_options::reset = false; |
| 4047 | 453 |
| 4049 | 454 restart = false; |
| 4047 | 455 } |
| 456 | |
| 4049 | 457 static double *dummy = 0; |
| 5275 | 458 static octave_idx_type *idummy = 0; |
| 4047 | 459 |
| 4049 | 460 F77_XFCN (ddaspk, DDASPK, (ddaspk_f, nn, t, px, pxdot, tout, pinfo, |
| 3998 | 461 prel_tol, pabs_tol, istate, prwork, lrw, |
| 3912 | 462 piwork, liw, dummy, idummy, ddaspk_j, |
| 463 ddaspk_psol)); | |
| 464 | |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
465 switch (istate) |
| 3912 | 466 { |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
467 case 1: // A step was successfully taken in intermediate-output |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
468 // mode. The code has not yet reached TOUT. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
469 case 2: // The integration to TSTOP was successfully completed |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
470 // (T=TSTOP) by stepping exactly to TSTOP. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
471 case 3: // The integration to TOUT was successfully completed |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
472 // (T=TOUT) by stepping past TOUT. Y(*) is obtained by |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
473 // interpolation. YPRIME(*) is obtained by interpolation. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
474 case 4: // The initial condition calculation, with |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
475 // INFO(11) > 0, was successful, and INFO(14) = 1. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
476 // No integration steps were taken, and the solution |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
477 // is not considered to have been started. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
478 retval = x; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
479 t = tout; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
480 break; |
| 3912 | 481 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
482 case -1: // A large amount of work has been expended. (~500 steps). |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
483 case -2: // The error tolerances are too stringent. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
484 case -3: // The local error test cannot be satisfied because you |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
485 // specified a zero component in ATOL and the |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
486 // corresponding computed solution component is zero. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
487 // Thus, a pure relative error test is impossible for |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
488 // this component. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
489 case -6: // DDASPK had repeated error test failures on the last |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
490 // attempted step. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
491 case -7: // The corrector could not converge. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
492 case -8: // The matrix of partial derivatives is singular. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
493 case -9: // The corrector could not converge. There were repeated |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
494 // error test failures in this step. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
495 case -10: // The corrector could not converge because IRES was |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
496 // equal to minus one. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
497 case -11: // IRES equal to -2 was encountered and control is being |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
498 // returned to the calling program. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
499 case -12: // DDASPK failed to compute the initial YPRIME. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
500 case -13: // Unrecoverable error encountered inside user's |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
501 // PSOL routine, and control is being returned to |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
502 // the calling program. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
503 case -14: // The Krylov linear system solver could not |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
504 // achieve convergence. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
505 case -33: // The code has encountered trouble from which it cannot |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
506 // recover. A message is printed explaining the trouble |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
507 // and control is returned to the calling program. For |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
508 // example, this occurs when invalid input is detected. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
509 integration_error = true; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
510 break; |
| 3996 | 511 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
512 default: |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
513 integration_error = true; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
514 (*current_liboctave_error_handler) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
515 ("unrecognized value of istate (= %d) returned from ddaspk", |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
516 istate); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
517 break; |
| 3912 | 518 } |
| 519 | |
| 520 return retval; | |
| 521 } | |
| 522 | |
| 523 Matrix | |
| 524 DASPK::do_integrate (const ColumnVector& tout) | |
| 525 { | |
| 526 Matrix dummy; | |
| 527 return integrate (tout, dummy); | |
| 528 } | |
| 529 | |
| 530 Matrix | |
| 531 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out) | |
| 532 { | |
| 533 Matrix retval; | |
| 4049 | 534 |
| 5275 | 535 octave_idx_type n_out = tout.capacity (); |
| 536 octave_idx_type n = size (); | |
| 3912 | 537 |
| 538 if (n_out > 0 && n > 0) | |
| 539 { | |
| 540 retval.resize (n_out, n); | |
| 541 xdot_out.resize (n_out, n); | |
| 542 | |
| 5275 | 543 for (octave_idx_type i = 0; i < n; i++) |
| 3912 | 544 { |
| 545 retval.elem (0, i) = x.elem (i); | |
| 546 xdot_out.elem (0, i) = xdot.elem (i); | |
| 547 } | |
| 548 | |
| 5275 | 549 for (octave_idx_type j = 1; j < n_out; j++) |
| 3912 | 550 { |
| 551 ColumnVector x_next = do_integrate (tout.elem (j)); | |
| 552 | |
| 553 if (integration_error) | |
| 554 return retval; | |
| 555 | |
| 5275 | 556 for (octave_idx_type i = 0; i < n; i++) |
| 3912 | 557 { |
| 558 retval.elem (j, i) = x_next.elem (i); | |
| 559 xdot_out.elem (j, i) = xdot.elem (i); | |
| 560 } | |
| 561 } | |
| 562 } | |
| 563 | |
| 564 return retval; | |
| 565 } | |
| 566 | |
| 567 Matrix | |
| 568 DASPK::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit) | |
| 569 { | |
| 570 Matrix dummy; | |
| 571 return integrate (tout, dummy, tcrit); | |
| 572 } | |
| 573 | |
| 574 Matrix | |
| 575 DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out, | |
| 576 const ColumnVector& tcrit) | |
| 577 { | |
| 578 Matrix retval; | |
| 4049 | 579 |
| 5275 | 580 octave_idx_type n_out = tout.capacity (); |
| 581 octave_idx_type n = size (); | |
| 3912 | 582 |
| 583 if (n_out > 0 && n > 0) | |
| 584 { | |
| 585 retval.resize (n_out, n); | |
| 586 xdot_out.resize (n_out, n); | |
| 587 | |
| 5275 | 588 for (octave_idx_type i = 0; i < n; i++) |
| 3912 | 589 { |
| 590 retval.elem (0, i) = x.elem (i); | |
| 591 xdot_out.elem (0, i) = xdot.elem (i); | |
| 592 } | |
| 593 | |
| 5275 | 594 octave_idx_type n_crit = tcrit.capacity (); |
| 3912 | 595 |
| 596 if (n_crit > 0) | |
| 597 { | |
| 5275 | 598 octave_idx_type i_crit = 0; |
| 599 octave_idx_type i_out = 1; | |
| 3912 | 600 double next_crit = tcrit.elem (0); |
| 601 double next_out; | |
| 602 while (i_out < n_out) | |
| 603 { | |
| 604 bool do_restart = false; | |
| 605 | |
| 606 next_out = tout.elem (i_out); | |
| 607 if (i_crit < n_crit) | |
| 608 next_crit = tcrit.elem (i_crit); | |
| 609 | |
| 610 bool save_output; | |
| 611 double t_out; | |
| 612 | |
| 613 if (next_crit == next_out) | |
| 614 { | |
| 615 set_stop_time (next_crit); | |
| 616 t_out = next_out; | |
| 617 save_output = true; | |
| 618 i_out++; | |
| 619 i_crit++; | |
| 620 do_restart = true; | |
| 621 } | |
| 622 else if (next_crit < next_out) | |
| 623 { | |
| 624 if (i_crit < n_crit) | |
| 625 { | |
| 626 set_stop_time (next_crit); | |
| 627 t_out = next_crit; | |
| 628 save_output = false; | |
| 629 i_crit++; | |
| 630 do_restart = true; | |
| 631 } | |
| 632 else | |
| 633 { | |
| 634 clear_stop_time (); | |
| 635 t_out = next_out; | |
| 636 save_output = true; | |
| 637 i_out++; | |
| 638 } | |
| 639 } | |
| 640 else | |
| 641 { | |
| 642 set_stop_time (next_crit); | |
| 643 t_out = next_out; | |
| 644 save_output = true; | |
| 645 i_out++; | |
| 646 } | |
| 647 | |
| 648 ColumnVector x_next = do_integrate (t_out); | |
| 649 | |
| 650 if (integration_error) | |
| 651 return retval; | |
| 652 | |
| 653 if (save_output) | |
| 654 { | |
| 5275 | 655 for (octave_idx_type i = 0; i < n; i++) |
| 3912 | 656 { |
| 657 retval.elem (i_out-1, i) = x_next.elem (i); | |
| 658 xdot_out.elem (i_out-1, i) = xdot.elem (i); | |
| 659 } | |
| 660 } | |
| 661 | |
| 662 if (do_restart) | |
| 663 force_restart (); | |
| 664 } | |
| 665 } | |
| 666 else | |
| 667 { | |
| 668 retval = integrate (tout, xdot_out); | |
| 669 | |
| 670 if (integration_error) | |
| 671 return retval; | |
| 672 } | |
| 673 } | |
| 674 | |
| 675 return retval; | |
| 676 } | |
| 677 | |
| 3995 | 678 std::string |
| 679 DASPK::error_message (void) const | |
| 680 { | |
| 681 std::string retval; | |
| 682 | |
| 5765 | 683 std::ostringstream buf; |
| 684 buf << t; | |
| 685 std::string t_curr = buf.str (); | |
| 4043 | 686 |
| 3997 | 687 switch (istate) |
| 3995 | 688 { |
| 3996 | 689 case 1: |
| 690 retval = "a step was successfully taken in intermediate-output mode."; | |
| 691 break; | |
| 692 | |
| 693 case 2: | |
| 694 retval = "integration completed by stepping exactly to TOUT"; | |
| 695 break; | |
| 696 | |
| 697 case 3: | |
| 698 retval = "integration to tout completed by stepping past TOUT"; | |
| 699 break; | |
| 700 | |
| 701 case 4: | |
| 702 retval = "initial condition calculation completed successfully"; | |
| 703 break; | |
| 704 | |
| 705 case -1: | |
| 4043 | 706 retval = std::string ("a large amount of work has been expended (t =") |
| 707 + t_curr + ")"; | |
| 3996 | 708 break; |
| 709 | |
| 710 case -2: | |
| 711 retval = "the error tolerances are too stringent"; | |
| 712 break; | |
| 713 | |
| 714 case -3: | |
| 4043 | 715 retval = std::string ("error weight became zero during problem. (t = ") |
| 716 + t_curr | |
| 717 + "; solution component i vanished, and atol or atol(i) == 0)"; | |
| 3996 | 718 break; |
| 719 | |
| 720 case -6: | |
| 4043 | 721 retval = std::string ("repeated error test failures on the last attempted step (t = ") |
| 722 + t_curr + ")"; | |
| 3996 | 723 break; |
| 724 | |
| 725 case -7: | |
| 4043 | 726 retval = std::string ("the corrector could not converge (t = ") |
| 727 + t_curr + ")"; | |
| 3996 | 728 break; |
| 729 | |
| 730 case -8: | |
| 4043 | 731 retval = std::string ("the matrix of partial derivatives is singular (t = ") |
| 732 + t_curr + ")"; | |
| 3996 | 733 break; |
| 734 | |
| 735 case -9: | |
| 4043 | 736 retval = std::string ("the corrector could not converge (t = ") |
| 737 + t_curr + "; repeated test failures)"; | |
| 3996 | 738 break; |
| 739 | |
| 740 case -10: | |
| 4043 | 741 retval = std::string ("corrector could not converge because IRES was -1 (t = ") |
| 742 + t_curr + ")"; | |
| 3996 | 743 break; |
| 744 | |
| 745 case -11: | |
| 4043 | 746 retval = std::string ("return requested in user-supplied function (t = ") |
| 747 + t_curr + ")"; | |
| 3996 | 748 break; |
| 749 | |
| 750 case -12: | |
| 751 retval = "failed to compute consistent initial conditions"; | |
| 752 break; | |
| 753 | |
| 754 case -13: | |
| 4043 | 755 retval = std::string ("unrecoverable error encountered inside user's PSOL function (t = ") |
| 756 + t_curr + ")"; | |
| 3996 | 757 break; |
| 758 | |
| 759 case -14: | |
| 4043 | 760 retval = std::string ("the Krylov linear system solver failed to converge (t = ") |
| 761 + t_curr + ")"; | |
| 3996 | 762 break; |
| 763 | |
| 764 case -33: | |
| 765 retval = "unrecoverable error (see printed message)"; | |
| 766 break; | |
| 767 | |
| 3995 | 768 default: |
| 769 retval = "unknown error state"; | |
| 770 break; | |
| 771 } | |
| 772 | |
| 773 return retval; | |
| 774 } | |
| 775 | |
| 3912 | 776 /* |
| 777 ;;; Local Variables: *** | |
| 778 ;;; mode: C++ *** | |
| 779 ;;; End: *** | |
| 780 */ |