Mercurial > hg > octave-lyh
annotate liboctave/LSODE.cc @ 14444:245963d3d628 stable
pkg: bug fix - accessing non-existent variable for error message
author | Miguel Bazdresch <lmb@2pif.info> |
---|---|
date | Thu, 01 Mar 2012 14:58:59 +0000 |
parents | 72c96de7a403 |
children | 3d8ace26c5b4 |
rev | line source |
---|---|
3 | 1 /* |
2 | |
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
12600
diff
changeset
|
3 Copyright (C) 1993-2012 John W. Eaton |
3 | 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. | |
3 | 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/>. | |
3 | 20 |
21 */ | |
22 | |
238 | 23 #ifdef HAVE_CONFIG_H |
1192 | 24 #include <config.h> |
3 | 25 #endif |
26 | |
1367 | 27 #include <cfloat> |
28 | |
5765 | 29 #include <sstream> |
30 | |
1842 | 31 #include "LSODE.h" |
1847 | 32 #include "f77-fcn.h" |
227 | 33 #include "lo-error.h" |
7231 | 34 #include "lo-math.h" |
4180 | 35 #include "quit.h" |
3 | 36 |
11495 | 37 typedef octave_idx_type (*lsode_fcn_ptr) (const octave_idx_type&, |
38 const double&, double*, | |
39 double*, octave_idx_type&); | |
3507 | 40 |
11495 | 41 typedef octave_idx_type (*lsode_jac_ptr) (const octave_idx_type&, |
42 const double&, double*, | |
43 const octave_idx_type&, | |
44 const octave_idx_type&, | |
45 double*, const octave_idx_type&); | |
3507 | 46 |
3 | 47 extern "C" |
4552 | 48 { |
49 F77_RET_T | |
11495 | 50 F77_FUNC (dlsode, DLSODE) (lsode_fcn_ptr, octave_idx_type&, double*, |
51 double&, double&, octave_idx_type&, double&, | |
52 const double*, octave_idx_type&, | |
53 octave_idx_type&, octave_idx_type&, | |
54 double*, octave_idx_type&, octave_idx_type*, | |
55 octave_idx_type&, lsode_jac_ptr, | |
56 octave_idx_type&); | |
4552 | 57 } |
3 | 58 |
532 | 59 static ODEFunc::ODERHSFunc user_fun; |
60 static ODEFunc::ODEJacFunc user_jac; | |
3 | 61 static ColumnVector *tmp_x; |
62 | |
5275 | 63 static octave_idx_type |
64 lsode_f (const octave_idx_type& neq, const double& time, double *, | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
65 double *deriv, octave_idx_type& ierr) |
3 | 66 { |
4180 | 67 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
68 | |
2343 | 69 ColumnVector tmp_deriv; |
3 | 70 |
1360 | 71 // NOTE: this won't work if LSODE passes copies of the state vector. |
72 // In that case we have to create a temporary vector object | |
73 // and copy. | |
74 | |
1251 | 75 tmp_deriv = (*user_fun) (*tmp_x, time); |
3 | 76 |
258 | 77 if (tmp_deriv.length () == 0) |
1251 | 78 ierr = -1; |
258 | 79 else |
80 { | |
5275 | 81 for (octave_idx_type i = 0; i < neq; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
82 deriv [i] = tmp_deriv.elem (i); |
258 | 83 } |
3 | 84 |
4180 | 85 END_INTERRUPT_WITH_EXCEPTIONS; |
86 | |
3 | 87 return 0; |
88 } | |
89 | |
5275 | 90 static octave_idx_type |
91 lsode_j (const octave_idx_type& neq, const double& time, double *, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
92 const octave_idx_type&, const octave_idx_type&, double *pd, const octave_idx_type& nrowpd) |
3 | 93 { |
4180 | 94 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
95 | |
1251 | 96 Matrix tmp_jac (neq, neq); |
3 | 97 |
1360 | 98 // NOTE: this won't work if LSODE passes copies of the state vector. |
99 // In that case we have to create a temporary vector object | |
100 // and copy. | |
101 | |
1251 | 102 tmp_jac = (*user_jac) (*tmp_x, time); |
3 | 103 |
5275 | 104 for (octave_idx_type j = 0; j < neq; j++) |
105 for (octave_idx_type i = 0; i < neq; i++) | |
1251 | 106 pd [nrowpd * j + i] = tmp_jac (i, j); |
3 | 107 |
4180 | 108 END_INTERRUPT_WITH_EXCEPTIONS; |
109 | |
3 | 110 return 0; |
111 } | |
112 | |
113 ColumnVector | |
1842 | 114 LSODE::do_integrate (double tout) |
3 | 115 { |
1945 | 116 ColumnVector retval; |
117 | |
5275 | 118 static octave_idx_type nn = 0; |
4049 | 119 |
120 if (! initialized || restart || ODEFunc::reset || LSODE_options::reset) | |
1945 | 121 { |
4049 | 122 integration_error = false; |
1945 | 123 |
4049 | 124 initialized = true; |
125 | |
126 istate = 1; | |
127 | |
5275 | 128 octave_idx_type n = size (); |
4049 | 129 |
130 nn = n; | |
3955 | 131 |
5275 | 132 octave_idx_type max_maxord = 0; |
4231 | 133 |
4049 | 134 if (integration_method () == "stiff") |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
135 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
136 max_maxord = 5; |
4231 | 137 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
138 if (jac) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
139 method_flag = 21; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
140 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
141 method_flag = 22; |
3955 | 142 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
143 liw = 20 + n; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
144 lrw = 22 + n * (9 + n); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
145 } |
4049 | 146 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
147 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
148 max_maxord = 12; |
4231 | 149 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
150 method_flag = 10; |
3955 | 151 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
152 liw = 20; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
153 lrw = 22 + 16 * n; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
154 } |
4049 | 155 |
4231 | 156 maxord = maximum_order (); |
157 | |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
158 iwork.resize (dim_vector (liw, 1)); |
5552 | 159 |
160 for (octave_idx_type i = 4; i < 9; i++) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
161 iwork(i) = 0; |
5552 | 162 |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
163 rwork.resize (dim_vector (lrw, 1)); |
5552 | 164 |
165 for (octave_idx_type i = 4; i < 9; i++) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
166 rwork(i) = 0; |
5552 | 167 |
4231 | 168 if (maxord >= 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
169 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
170 if (maxord > 0 && maxord <= max_maxord) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
171 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
172 iwork(4) = maxord; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
173 iopt = 1; |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
174 } |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
175 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
176 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
177 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
178 ("lsode: invalid value for maximum order"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
179 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
180 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
181 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
182 } |
4231 | 183 |
4049 | 184 if (stop_time_set) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
185 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
186 itask = 4; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
187 rwork(0) = stop_time; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
188 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
189 } |
4049 | 190 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
191 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
192 itask = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
193 } |
258 | 194 |
4049 | 195 restart = false; |
196 | |
197 // ODEFunc | |
3 | 198 |
4049 | 199 // NOTE: this won't work if LSODE passes copies of the state vector. |
200 // In that case we have to create a temporary vector object | |
201 // and copy. | |
3 | 202 |
4049 | 203 tmp_x = &x; |
204 | |
205 user_fun = function (); | |
206 user_jac = jacobian_function (); | |
207 | |
2343 | 208 ColumnVector xdot = (*user_fun) (x, t); |
209 | |
210 if (x.length () != xdot.length ()) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
211 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
212 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
213 ("lsode: inconsistent sizes for state and derivative vectors"); |
2343 | 214 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
215 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
216 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
217 } |
2343 | 218 |
4049 | 219 ODEFunc::reset = false; |
220 | |
221 // LSODE_options | |
222 | |
223 rel_tol = relative_tolerance (); | |
224 abs_tol = absolute_tolerance (); | |
225 | |
5275 | 226 octave_idx_type abs_tol_len = abs_tol.length (); |
4049 | 227 |
228 if (abs_tol_len == 1) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
229 itol = 1; |
4049 | 230 else if (abs_tol_len == n) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
231 itol = 2; |
4049 | 232 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
233 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
234 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
235 ("lsode: inconsistent sizes for state and absolute tolerance vectors"); |
4049 | 236 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
237 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
238 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
239 } |
2343 | 240 |
4049 | 241 double iss = initial_step_size (); |
242 if (iss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
243 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
244 rwork(4) = iss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 } |
4049 | 247 |
248 double maxss = maximum_step_size (); | |
249 if (maxss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
250 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
251 rwork(5) = maxss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
252 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
253 } |
4049 | 254 |
255 double minss = minimum_step_size (); | |
256 if (minss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
257 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
258 rwork(6) = minss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
259 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
260 } |
4049 | 261 |
5275 | 262 octave_idx_type sl = step_limit (); |
4049 | 263 if (sl > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
264 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
265 iwork(5) = sl; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
266 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
267 } |
4049 | 268 |
269 LSODE_options::reset = false; | |
3 | 270 } |
271 | |
11502
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
272 double *px = x.fortran_vec (); |
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
273 |
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
274 double *pabs_tol = abs_tol.fortran_vec (); |
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
275 |
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
276 octave_idx_type *piwork = iwork.fortran_vec (); |
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
277 double *prwork = rwork.fortran_vec (); |
4638800cd660
delete data pointer members from liboctave ODE/DAE classes; make destuctors virtual in ODE/DAE base classes
John W. Eaton <jwe@octave.org>
parents:
11495
diff
changeset
|
278 |
4583 | 279 F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
280 pabs_tol, itask, istate, iopt, prwork, lrw, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
281 piwork, liw, lsode_j, method_flag)); |
3 | 282 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
283 switch (istate) |
3 | 284 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
285 case 1: // prior to initial integration step. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
286 case 2: // lsode was successful. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
287 retval = x; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
288 t = tout; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
289 break; |
3996 | 290 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
291 case -1: // excess work done on this call (perhaps wrong mf). |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
292 case -2: // excess accuracy requested (tolerances too small). |
8807
401d54a83690
use 'invalid', not 'illegal'
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
293 case -3: // invalid input detected (see printed message). |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
294 case -4: // repeated error test failures (check all inputs). |
12600
7c000c70f873
LSODE.cc: Add semicolon to error messages to prevent run-together text.
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
295 case -5: // repeated convergence failures (perhaps bad Jacobian |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
296 // supplied or wrong choice of mf or tolerances). |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
297 case -6: // error weight became zero during problem. (solution |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
298 // component i vanished, and atol or atol(i) = 0.) |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
299 case -13: // return requested in user-supplied function. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
300 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
301 break; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
302 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
303 default: |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
304 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
305 (*current_liboctave_error_handler) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
306 ("unrecognized value of istate (= %d) returned from lsode", |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
307 istate); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
308 break; |
3 | 309 } |
310 | |
1945 | 311 return retval; |
3 | 312 } |
313 | |
3959 | 314 std::string |
315 LSODE::error_message (void) const | |
316 { | |
317 std::string retval; | |
318 | |
5765 | 319 std::ostringstream buf; |
320 buf << t; | |
321 std::string t_curr = buf.str (); | |
4042 | 322 |
3959 | 323 switch (istate) |
324 { | |
325 case 1: | |
3996 | 326 retval = "prior to initial integration step"; |
3959 | 327 break; |
328 | |
329 case 2: | |
330 retval = "successful exit"; | |
331 break; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
332 |
3959 | 333 case 3: |
334 retval = "prior to continuation call with modified parameters"; | |
335 break; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
336 |
3996 | 337 case -1: |
4042 | 338 retval = std::string ("excess work on this call (t = ") |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
339 + t_curr + "; perhaps wrong integration method)"; |
3996 | 340 break; |
341 | |
342 case -2: | |
343 retval = "excess accuracy requested (tolerances too small)"; | |
344 break; | |
345 | |
346 case -3: | |
347 retval = "invalid input detected (see printed message)"; | |
348 break; | |
349 | |
350 case -4: | |
4042 | 351 retval = std::string ("repeated error test failures (t = ") |
12600
7c000c70f873
LSODE.cc: Add semicolon to error messages to prevent run-together text.
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
352 + t_curr + "; check all inputs)"; |
3996 | 353 break; |
354 | |
355 case -5: | |
4042 | 356 retval = std::string ("repeated convergence failures (t = ") |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
357 + t_curr |
12600
7c000c70f873
LSODE.cc: Add semicolon to error messages to prevent run-together text.
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
358 + "; perhaps bad Jacobian supplied or wrong choice of integration method or tolerances)"; |
3996 | 359 break; |
360 | |
361 case -6: | |
4042 | 362 retval = std::string ("error weight became zero during problem. (t = ") |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
363 + t_curr |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
364 + "; solution component i vanished, and atol or atol(i) == 0)"; |
3996 | 365 break; |
366 | |
367 case -13: | |
4042 | 368 retval = "return requested in user-supplied function (t = " |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
369 + t_curr + ")"; |
3996 | 370 break; |
371 | |
3959 | 372 default: |
373 retval = "unknown error state"; | |
374 break; | |
375 } | |
376 | |
377 return retval; | |
378 } | |
379 | |
3 | 380 Matrix |
1842 | 381 LSODE::do_integrate (const ColumnVector& tout) |
3 | 382 { |
383 Matrix retval; | |
4049 | 384 |
5275 | 385 octave_idx_type n_out = tout.capacity (); |
386 octave_idx_type n = size (); | |
3 | 387 |
388 if (n_out > 0 && n > 0) | |
389 { | |
390 retval.resize (n_out, n); | |
391 | |
5275 | 392 for (octave_idx_type i = 0; i < n; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
393 retval.elem (0, i) = x.elem (i); |
3 | 394 |
5275 | 395 for (octave_idx_type j = 1; j < n_out; j++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
396 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
397 ColumnVector x_next = do_integrate (tout.elem (j)); |
258 | 398 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
399 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
400 return retval; |
258 | 401 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
402 for (octave_idx_type i = 0; i < n; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
403 retval.elem (j, i) = x_next.elem (i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
404 } |
3 | 405 } |
406 | |
407 return retval; | |
408 } | |
409 | |
410 Matrix | |
3511 | 411 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit) |
3 | 412 { |
413 Matrix retval; | |
4049 | 414 |
5275 | 415 octave_idx_type n_out = tout.capacity (); |
416 octave_idx_type n = size (); | |
3 | 417 |
418 if (n_out > 0 && n > 0) | |
419 { | |
420 retval.resize (n_out, n); | |
421 | |
5275 | 422 for (octave_idx_type i = 0; i < n; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
423 retval.elem (0, i) = x.elem (i); |
3 | 424 |
5275 | 425 octave_idx_type n_crit = tcrit.capacity (); |
3 | 426 |
427 if (n_crit > 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
428 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
429 octave_idx_type i_crit = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
430 octave_idx_type i_out = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
431 double next_crit = tcrit.elem (0); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
432 double next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
433 while (i_out < n_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
434 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
435 bool do_restart = false; |
3 | 436 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
437 next_out = tout.elem (i_out); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
438 if (i_crit < n_crit) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
439 next_crit = tcrit.elem (i_crit); |
3 | 440 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
441 octave_idx_type save_output; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
442 double t_out; |
3 | 443 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
444 if (next_crit == next_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
445 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
446 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
447 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
448 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
449 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
450 i_crit++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
451 do_restart = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
452 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
453 else if (next_crit < next_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
454 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
455 if (i_crit < n_crit) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
456 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
457 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
458 t_out = next_crit; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
459 save_output = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
460 i_crit++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
461 do_restart = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
462 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
463 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
464 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
465 clear_stop_time (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
466 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
467 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
468 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
469 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
470 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
471 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
472 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
473 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
474 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
475 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
476 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
477 } |
3 | 478 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
479 ColumnVector x_next = do_integrate (t_out); |
3 | 480 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
481 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
482 return retval; |
258 | 483 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
484 if (save_output) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
485 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
486 for (octave_idx_type i = 0; i < n; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
487 retval.elem (i_out-1, i) = x_next.elem (i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
488 } |
3 | 489 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
490 if (do_restart) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
491 force_restart (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
492 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
493 } |
3 | 494 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
495 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
496 retval = do_integrate (tout); |
258 | 497 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
498 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
499 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
500 } |
3 | 501 } |
502 | |
503 return retval; | |
504 } |