Mercurial > octave-libgccjit
annotate liboctave/DASRT.cc @ 8920:eb63fbe60fab
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 10:41:27 -0500 |
| parents | 29980c6b8604 |
| children | 4c0cdbe0acca |
| rev | line source |
|---|---|
| 3990 | 1 /* |
| 2 | |
| 8920 | 3 Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008 John W. Eaton |
| 3990 | 4 |
| 5 This file is part of Octave. | |
| 6 | |
| 7 Octave is free software; you can redistribute it and/or modify it | |
| 8 under the terms of the GNU General Public License as published by the | |
| 7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
| 10 option) any later version. | |
| 3990 | 11 |
| 12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 15 for more details. | |
| 16 | |
| 17 You should have received a copy of the GNU General Public License | |
| 7016 | 18 along with Octave; see the file COPYING. If not, see |
| 19 <http://www.gnu.org/licenses/>. | |
| 3990 | 20 |
| 21 */ | |
| 22 | |
| 23 #ifdef HAVE_CONFIG_H | |
| 24 #include <config.h> | |
| 25 #endif | |
| 26 | |
| 27 #include <cfloat> | |
| 28 | |
| 5765 | 29 #include <sstream> |
| 30 | |
| 3990 | 31 #include "DASRT.h" |
| 32 #include "f77-fcn.h" | |
| 33 #include "lo-error.h" | |
| 7231 | 34 #include "lo-math.h" |
| 4180 | 35 #include "quit.h" |
| 3990 | 36 |
| 5275 | 37 typedef octave_idx_type (*dasrt_fcn_ptr) (const double&, const double*, const double*, |
| 38 double*, octave_idx_type&, double*, octave_idx_type*); | |
| 3991 | 39 |
| 5275 | 40 typedef octave_idx_type (*dasrt_jac_ptr) (const double&, const double*, const double*, |
| 41 double*, const double&, double*, octave_idx_type*); | |
| 3991 | 42 |
| 5275 | 43 typedef octave_idx_type (*dasrt_constr_ptr) (const octave_idx_type&, const double&, const double*, |
| 44 const octave_idx_type&, double*, double*, octave_idx_type*); | |
| 3991 | 45 |
| 3990 | 46 extern "C" |
| 4552 | 47 { |
| 48 F77_RET_T | |
| 5275 | 49 F77_FUNC (ddasrt, DDASRT) (dasrt_fcn_ptr, const octave_idx_type&, double&, |
| 50 double*, double*, const double&, octave_idx_type*, | |
| 51 const double*, const double*, octave_idx_type&, double*, | |
| 52 const octave_idx_type&, octave_idx_type*, const octave_idx_type&, double*, | |
| 53 octave_idx_type*, dasrt_jac_ptr, dasrt_constr_ptr, | |
| 54 const octave_idx_type&, octave_idx_type*); | |
| 4552 | 55 } |
| 3990 | 56 |
| 57 static DAEFunc::DAERHSFunc user_fsub; | |
| 58 static DAEFunc::DAEJacFunc user_jsub; | |
| 59 static DAERTFunc::DAERTConstrFunc user_csub; | |
| 3993 | 60 |
| 5275 | 61 static octave_idx_type nn; |
| 3990 | 62 |
| 5275 | 63 static octave_idx_type |
| 3991 | 64 ddasrt_f (const double& t, const double *state, const double *deriv, |
| 5275 | 65 double *delta, octave_idx_type& ires, double *, octave_idx_type *) |
| 3990 | 66 { |
| 4180 | 67 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 68 | |
| 3990 | 69 ColumnVector tmp_state (nn); |
| 4133 | 70 ColumnVector tmp_deriv (nn); |
| 3990 | 71 |
| 5275 | 72 for (octave_idx_type i = 0; i < nn; i++) |
| 4133 | 73 { |
| 74 tmp_state(i) = state[i]; | |
| 75 tmp_deriv(i) = deriv[i]; | |
| 76 } | |
| 3990 | 77 |
| 3994 | 78 ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires); |
| 3990 | 79 |
| 80 if (tmp_fval.length () == 0) | |
| 81 ires = -2; | |
| 82 else | |
| 83 { | |
| 5275 | 84 for (octave_idx_type i = 0; i < nn; i++) |
| 3990 | 85 delta[i] = tmp_fval(i); |
| 86 } | |
| 87 | |
| 4180 | 88 END_INTERRUPT_WITH_EXCEPTIONS; |
| 89 | |
| 3990 | 90 return 0; |
| 91 } | |
| 92 | |
| 5275 | 93 octave_idx_type |
| 3991 | 94 ddasrt_j (const double& time, const double *state, const double *deriv, |
| 5275 | 95 double *pd, const double& cj, double *, octave_idx_type *) |
| 3990 | 96 { |
| 4180 | 97 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 98 | |
| 5775 | 99 // FIXME -- would be nice to avoid copying the data. |
| 3990 | 100 |
| 3991 | 101 ColumnVector tmp_state (nn); |
| 102 ColumnVector tmp_deriv (nn); | |
| 3990 | 103 |
| 5275 | 104 for (octave_idx_type i = 0; i < nn; i++) |
| 3991 | 105 { |
| 106 tmp_deriv.elem (i) = deriv [i]; | |
| 107 tmp_state.elem (i) = state [i]; | |
| 108 } | |
| 3990 | 109 |
| 3994 | 110 Matrix tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj); |
| 3990 | 111 |
| 5275 | 112 for (octave_idx_type j = 0; j < nn; j++) |
| 113 for (octave_idx_type i = 0; i < nn; i++) | |
| 3991 | 114 pd [nn * j + i] = tmp_pd.elem (i, j); |
| 3990 | 115 |
| 4180 | 116 END_INTERRUPT_WITH_EXCEPTIONS; |
| 117 | |
| 3990 | 118 return 0; |
| 119 } | |
| 120 | |
| 5275 | 121 static octave_idx_type |
| 122 ddasrt_g (const octave_idx_type& neq, const double& t, const double *state, | |
| 123 const octave_idx_type& ng, double *gout, double *, octave_idx_type *) | |
| 3990 | 124 { |
| 4180 | 125 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 126 | |
| 5275 | 127 octave_idx_type n = neq; |
| 3990 | 128 |
| 129 ColumnVector tmp_state (n); | |
| 5275 | 130 for (octave_idx_type i = 0; i < n; i++) |
| 3990 | 131 tmp_state(i) = state[i]; |
| 132 | |
| 3994 | 133 ColumnVector tmp_fval = (*user_csub) (tmp_state, t); |
| 3990 | 134 |
| 5275 | 135 for (octave_idx_type i = 0; i < ng; i++) |
| 3990 | 136 gout[i] = tmp_fval(i); |
| 137 | |
| 4180 | 138 END_INTERRUPT_WITH_EXCEPTIONS; |
| 139 | |
| 3990 | 140 return 0; |
| 141 } | |
| 142 | |
| 143 void | |
| 144 DASRT::integrate (double tout) | |
| 145 { | |
| 146 DASRT_result retval; | |
| 147 | |
| 4049 | 148 // I suppose this is the safe thing to do. If this is the first |
| 149 // call, or if anything about the problem has changed, we should | |
| 150 // start completely fresh. | |
| 151 | |
| 152 if (! initialized || restart | |
| 153 || DAEFunc::reset || DAERTFunc::reset || DASRT_options::reset) | |
| 3990 | 154 { |
| 155 integration_error = false; | |
| 156 | |
| 4049 | 157 initialized = true; |
| 158 | |
| 159 info.resize (15); | |
| 160 | |
| 5275 | 161 for (octave_idx_type i = 0; i < 15; i++) |
| 4049 | 162 info(i) = 0; |
| 163 | |
| 164 pinfo = info.fortran_vec (); | |
| 165 | |
| 5275 | 166 octave_idx_type n = size (); |
| 4049 | 167 |
| 168 nn = n; | |
| 3990 | 169 |
| 4133 | 170 // DAERTFunc |
| 171 | |
| 172 user_csub = DAERTFunc::constraint_function (); | |
| 173 | |
| 174 if (user_csub) | |
| 175 { | |
| 176 ColumnVector tmp = (*user_csub) (x, t); | |
| 177 ng = tmp.length (); | |
| 178 } | |
| 179 else | |
| 180 ng = 0; | |
| 181 | |
| 5275 | 182 octave_idx_type maxord = maximum_order (); |
| 4133 | 183 if (maxord >= 0) |
| 184 { | |
| 185 if (maxord > 0 && maxord < 6) | |
| 186 { | |
| 187 info(8) = 1; | |
| 188 iwork(2) = maxord; | |
| 189 } | |
| 190 else | |
| 191 { | |
| 192 (*current_liboctave_error_handler) | |
| 193 ("dassl: invalid value for maximum order"); | |
| 194 integration_error = true; | |
| 195 return; | |
| 196 } | |
| 197 } | |
| 198 | |
| 4428 | 199 liw = 21 + n; |
| 4133 | 200 lrw = 50 + 9*n + n*n + 3*ng; |
| 4049 | 201 |
| 202 iwork.resize (liw); | |
| 203 rwork.resize (lrw); | |
| 204 | |
| 205 info(0) = 0; | |
| 206 | |
| 207 if (stop_time_set) | |
| 208 { | |
| 209 info(3) = 1; | |
| 210 rwork(0) = stop_time; | |
| 211 } | |
| 3990 | 212 else |
| 4049 | 213 info(3) = 0; |
| 3990 | 214 |
| 3994 | 215 px = x.fortran_vec (); |
| 216 pxdot = xdot.fortran_vec (); | |
| 217 | |
| 4049 | 218 piwork = iwork.fortran_vec (); |
| 219 prwork = rwork.fortran_vec (); | |
| 220 | |
| 221 restart = false; | |
| 3994 | 222 |
| 4049 | 223 // DAEFunc |
| 224 | |
| 225 user_fsub = DAEFunc::function (); | |
| 226 user_jsub = DAEFunc::jacobian_function (); | |
| 227 | |
| 228 if (user_fsub) | |
| 3990 | 229 { |
| 5275 | 230 octave_idx_type ires = 0; |
| 3990 | 231 |
| 3994 | 232 ColumnVector fval = (*user_fsub) (x, xdot, t, ires); |
| 3990 | 233 |
| 234 if (fval.length () != x.length ()) | |
| 235 { | |
| 236 (*current_liboctave_error_handler) | |
| 4049 | 237 ("dasrt: inconsistent sizes for state and residual vectors"); |
| 3990 | 238 |
| 239 integration_error = true; | |
| 240 return; | |
| 241 } | |
| 4049 | 242 } |
| 243 else | |
| 244 { | |
| 245 (*current_liboctave_error_handler) | |
| 246 ("dasrt: no user supplied RHS subroutine!"); | |
| 3990 | 247 |
| 4049 | 248 integration_error = true; |
| 249 return; | |
| 250 } | |
| 251 | |
| 252 info(4) = user_jsub ? 1 : 0; | |
| 253 | |
| 254 DAEFunc::reset = false; | |
| 255 | |
| 4625 | 256 jroot.resize (ng, 1); |
| 4049 | 257 |
| 258 pjroot = jroot.fortran_vec (); | |
| 259 | |
| 260 DAERTFunc::reset = false; | |
| 261 | |
| 262 // DASRT_options | |
| 263 | |
| 4050 | 264 double mss = maximum_step_size (); |
| 265 if (mss >= 0.0) | |
| 4049 | 266 { |
| 4050 | 267 rwork(1) = mss; |
| 4049 | 268 info(6) = 1; |
| 269 } | |
| 270 else | |
| 271 info(6) = 0; | |
| 3990 | 272 |
| 4050 | 273 double iss = initial_step_size (); |
| 274 if (iss >= 0.0) | |
| 4049 | 275 { |
| 4050 | 276 rwork(2) = iss; |
| 4049 | 277 info(7) = 1; |
| 278 } | |
| 279 else | |
| 280 info(7) = 0; | |
| 281 | |
| 282 if (step_limit () >= 0) | |
| 283 { | |
| 284 info(11) = 1; | |
| 4428 | 285 iwork(20) = step_limit (); |
| 4049 | 286 } |
| 287 else | |
| 288 info(11) = 0; | |
| 3990 | 289 |
| 290 abs_tol = absolute_tolerance (); | |
| 291 rel_tol = relative_tolerance (); | |
| 292 | |
| 5275 | 293 octave_idx_type abs_tol_len = abs_tol.length (); |
| 294 octave_idx_type rel_tol_len = rel_tol.length (); | |
| 3998 | 295 |
| 296 if (abs_tol_len == 1 && rel_tol_len == 1) | |
| 297 { | |
| 298 info.elem (1) = 0; | |
| 299 } | |
| 300 else if (abs_tol_len == n && rel_tol_len == n) | |
| 301 { | |
| 302 info.elem (1) = 1; | |
| 303 } | |
| 304 else | |
| 305 { | |
| 306 (*current_liboctave_error_handler) | |
| 4049 | 307 ("dasrt: inconsistent sizes for tolerance arrays"); |
| 3998 | 308 |
| 309 integration_error = true; | |
| 310 return; | |
| 311 } | |
| 312 | |
| 313 pabs_tol = abs_tol.fortran_vec (); | |
| 314 prel_tol = rel_tol.fortran_vec (); | |
| 3990 | 315 |
| 4049 | 316 DASRT_options::reset = false; |
| 3990 | 317 } |
| 318 | |
| 4049 | 319 static double *dummy = 0; |
| 5275 | 320 static octave_idx_type *idummy = 0; |
| 3990 | 321 |
| 4575 | 322 F77_XFCN (ddasrt, DDASRT, (ddasrt_f, nn, t, px, pxdot, tout, pinfo, |
| 323 prel_tol, pabs_tol, istate, prwork, lrw, | |
| 324 piwork, liw, dummy, idummy, ddasrt_j, | |
| 325 ddasrt_g, ng, pjroot)); | |
| 3990 | 326 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
327 switch (istate) |
| 3990 | 328 { |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
329 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
|
330 // 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
|
331 case 2: // The integration to TOUT was successfully completed |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
332 // (T=TOUT) by stepping exactly to TOUT. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
333 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
|
334 // (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
|
335 // interpolation. YPRIME(*) is obtained by interpolation. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
336 t = tout; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
337 break; |
| 3990 | 338 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
339 case 4: // The integration was successfully completed |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
340 // by finding one or more roots of G at T. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
341 break; |
| 3990 | 342 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
343 case -1: // A large amount of work has been expended. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
344 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
|
345 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
|
346 // 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
|
347 // corresponding computed solution component is zero. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
348 // 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
|
349 // this component. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
350 case -6: // DDASRT 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
|
351 // attempted step. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
352 case -7: // The corrector could not converge. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
353 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
|
354 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
|
355 // error test failures in this step. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
356 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
|
357 // equal to minus one. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
358 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
|
359 // returned to the calling program. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
360 case -12: // DASSL failed to compute the initial YPRIME. |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
361 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
|
362 // 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
|
363 // 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
|
364 // 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
|
365 integration_error = true; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
366 break; |
| 3996 | 367 |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
368 default: |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
369 integration_error = true; |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
370 (*current_liboctave_error_handler) |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
371 ("unrecognized value of istate (= %d) returned from ddasrt", |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
372 istate); |
|
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
373 break; |
| 3990 | 374 } |
| 375 } | |
| 376 | |
| 377 DASRT_result | |
| 378 DASRT::integrate (const ColumnVector& tout) | |
| 379 { | |
| 380 DASRT_result retval; | |
| 381 | |
| 382 Matrix x_out; | |
| 383 Matrix xdot_out; | |
| 3994 | 384 ColumnVector t_out = tout; |
| 3990 | 385 |
| 5275 | 386 octave_idx_type n_out = tout.capacity (); |
| 387 octave_idx_type n = size (); | |
| 3990 | 388 |
| 389 if (n_out > 0 && n > 0) | |
| 390 { | |
| 391 x_out.resize (n_out, n); | |
| 392 xdot_out.resize (n_out, n); | |
| 393 | |
| 5275 | 394 for (octave_idx_type i = 0; i < n; i++) |
| 3994 | 395 { |
| 396 x_out(0,i) = x(i); | |
| 397 xdot_out(0,i) = xdot(i); | |
| 398 } | |
| 399 | |
| 5275 | 400 for (octave_idx_type j = 1; j < n_out; j++) |
| 3990 | 401 { |
| 402 integrate (tout(j)); | |
| 3994 | 403 |
| 3990 | 404 if (integration_error) |
| 405 { | |
| 406 retval = DASRT_result (x_out, xdot_out, t_out); | |
| 407 return retval; | |
| 408 } | |
| 3992 | 409 |
| 3997 | 410 if (istate == 4) |
| 3990 | 411 t_out(j) = t; |
| 412 else | |
| 413 t_out(j) = tout(j); | |
| 414 | |
| 5275 | 415 for (octave_idx_type i = 0; i < n; i++) |
| 3990 | 416 { |
| 3994 | 417 x_out(j,i) = x(i); |
| 418 xdot_out(j,i) = xdot(i); | |
| 3990 | 419 } |
| 3992 | 420 |
| 3997 | 421 if (istate == 4) |
| 3992 | 422 { |
| 3994 | 423 x_out.resize (j+1, n); |
| 424 xdot_out.resize (j+1, n); | |
| 425 t_out.resize (j+1); | |
| 426 break; | |
| 3992 | 427 } |
| 3990 | 428 } |
| 429 } | |
| 430 | |
| 431 retval = DASRT_result (x_out, xdot_out, t_out); | |
| 432 | |
| 433 return retval; | |
| 434 } | |
| 435 | |
| 436 DASRT_result | |
| 437 DASRT::integrate (const ColumnVector& tout, const ColumnVector& tcrit) | |
| 438 { | |
| 439 DASRT_result retval; | |
| 440 | |
| 441 Matrix x_out; | |
| 442 Matrix xdot_out; | |
| 3994 | 443 ColumnVector t_outs = tout; |
| 3990 | 444 |
| 5275 | 445 octave_idx_type n_out = tout.capacity (); |
| 446 octave_idx_type n = size (); | |
| 3990 | 447 |
| 448 if (n_out > 0 && n > 0) | |
| 449 { | |
| 450 x_out.resize (n_out, n); | |
| 451 xdot_out.resize (n_out, n); | |
| 452 | |
| 5275 | 453 octave_idx_type n_crit = tcrit.capacity (); |
| 3990 | 454 |
| 455 if (n_crit > 0) | |
| 456 { | |
| 5275 | 457 octave_idx_type i_crit = 0; |
| 458 octave_idx_type i_out = 1; | |
| 3990 | 459 double next_crit = tcrit(0); |
| 460 double next_out; | |
| 461 while (i_out < n_out) | |
| 462 { | |
| 463 bool do_restart = false; | |
| 464 | |
| 465 next_out = tout(i_out); | |
| 466 if (i_crit < n_crit) | |
| 467 next_crit = tcrit(i_crit); | |
| 468 | |
| 5275 | 469 octave_idx_type save_output; |
| 3990 | 470 double t_out; |
| 471 | |
| 472 if (next_crit == next_out) | |
| 473 { | |
| 474 set_stop_time (next_crit); | |
| 475 t_out = next_out; | |
| 476 save_output = 1; | |
| 477 i_out++; | |
| 478 i_crit++; | |
| 479 do_restart = true; | |
| 480 } | |
| 481 else if (next_crit < next_out) | |
| 482 { | |
| 483 if (i_crit < n_crit) | |
| 484 { | |
| 485 set_stop_time (next_crit); | |
| 486 t_out = next_crit; | |
| 487 save_output = 0; | |
| 488 i_crit++; | |
| 489 do_restart = true; | |
| 490 } | |
| 491 else | |
| 492 { | |
| 493 clear_stop_time (); | |
| 494 t_out = next_out; | |
| 495 save_output = 1; | |
| 496 i_out++; | |
| 497 } | |
| 498 } | |
| 499 else | |
| 500 { | |
| 501 set_stop_time (next_crit); | |
| 502 t_out = next_out; | |
| 503 save_output = 1; | |
| 504 i_out++; | |
| 505 } | |
| 506 | |
| 507 integrate (t_out); | |
| 508 | |
| 509 if (integration_error) | |
| 510 { | |
| 511 retval = DASRT_result (x_out, xdot_out, t_outs); | |
| 512 return retval; | |
| 513 } | |
| 3992 | 514 |
| 3997 | 515 if (istate == 4) |
| 3990 | 516 t_out = t; |
| 517 | |
| 518 if (save_output) | |
| 519 { | |
| 5275 | 520 for (octave_idx_type i = 0; i < n; i++) |
| 3990 | 521 { |
| 3994 | 522 x_out(i_out-1,i) = x(i); |
| 523 xdot_out(i_out-1,i) = xdot(i); | |
| 3990 | 524 } |
| 3994 | 525 |
| 3990 | 526 t_outs(i_out-1) = t_out; |
| 3994 | 527 |
| 3997 | 528 if (istate == 4) |
| 3990 | 529 { |
| 530 x_out.resize (i_out, n); | |
| 531 xdot_out.resize (i_out, n); | |
| 532 t_outs.resize (i_out); | |
| 533 i_out = n_out; | |
| 534 } | |
| 535 } | |
| 536 | |
| 537 if (do_restart) | |
| 538 force_restart (); | |
| 539 } | |
| 540 | |
| 541 retval = DASRT_result (x_out, xdot_out, t_outs); | |
| 542 } | |
| 543 else | |
| 544 { | |
| 545 retval = integrate (tout); | |
| 546 | |
| 547 if (integration_error) | |
| 548 return retval; | |
| 549 } | |
| 550 } | |
| 551 | |
| 552 return retval; | |
| 553 } | |
| 554 | |
| 3995 | 555 std::string |
| 556 DASRT::error_message (void) const | |
| 557 { | |
| 558 std::string retval; | |
| 559 | |
| 5765 | 560 std::ostringstream buf; |
| 561 buf << t; | |
| 562 std::string t_curr = buf.str (); | |
| 4043 | 563 |
| 3997 | 564 switch (istate) |
| 3995 | 565 { |
| 3996 | 566 case 1: |
| 567 retval = "a step was successfully taken in intermediate-output mode."; | |
| 568 break; | |
| 569 | |
| 570 case 2: | |
| 571 retval = "integration completed by stepping exactly to TOUT"; | |
| 572 break; | |
| 573 | |
| 574 case 3: | |
| 575 retval = "integration to tout completed by stepping past TOUT"; | |
| 576 break; | |
| 577 | |
| 578 case 4: | |
| 579 retval = "integration completed by finding one or more roots of G at T"; | |
| 580 break; | |
| 581 | |
| 582 case -1: | |
| 4043 | 583 retval = std::string ("a large amount of work has been expended (t =") |
| 584 + t_curr + ")"; | |
| 3996 | 585 break; |
| 586 | |
| 587 case -2: | |
| 588 retval = "the error tolerances are too stringent"; | |
| 589 break; | |
| 590 | |
| 591 case -3: | |
| 4043 | 592 retval = std::string ("error weight became zero during problem. (t = ") |
| 593 + t_curr | |
| 594 + "; solution component i vanished, and atol or atol(i) == 0)"; | |
| 3996 | 595 break; |
| 596 | |
| 597 case -6: | |
| 4043 | 598 retval = std::string ("repeated error test failures on the last attempted step (t = ") |
| 599 + t_curr + ")"; | |
| 3996 | 600 break; |
| 601 | |
| 602 case -7: | |
| 4043 | 603 retval = std::string ("the corrector could not converge (t = ") |
| 604 + t_curr + ")"; | |
| 3996 | 605 break; |
| 606 | |
| 607 case -8: | |
| 4043 | 608 retval = std::string ("the matrix of partial derivatives is singular (t = ") |
| 609 + t_curr + ")"; | |
| 3996 | 610 break; |
| 611 | |
| 612 case -9: | |
| 4043 | 613 retval = std::string ("the corrector could not converge (t = ") |
| 614 + t_curr + "; repeated test failures)"; | |
| 3996 | 615 break; |
| 616 | |
| 617 case -10: | |
| 4043 | 618 retval = std::string ("corrector could not converge because IRES was -1 (t = ") |
| 619 + t_curr + ")"; | |
| 3996 | 620 break; |
| 621 | |
| 622 case -11: | |
| 4043 | 623 retval = std::string ("return requested in user-supplied function (t = ") |
| 624 + t_curr + ")"; | |
| 3996 | 625 break; |
| 626 | |
| 627 case -12: | |
| 628 retval = "failed to compute consistent initial conditions"; | |
| 629 break; | |
| 630 | |
| 631 case -33: | |
| 632 retval = "unrecoverable error (see printed message)"; | |
| 633 break; | |
| 634 | |
| 3995 | 635 default: |
| 636 retval = "unknown error state"; | |
| 637 break; | |
| 638 } | |
| 639 | |
| 640 return retval; | |
| 641 } | |
| 642 | |
| 3990 | 643 /* |
| 644 ;;; Local Variables: *** | |
| 645 ;;; mode: C++ *** | |
| 646 ;;; End: *** | |
| 647 */ |