Mercurial > hg > octave-nkf
annotate liboctave/numeric/LSODE.cc @ 19139:d5db1019ad4a
load-save.cc: Return string rather printing to octave_stdout
author | Andreas Weber <andy.weber.aw@gmail.com> |
---|---|
date | Mon, 21 Jul 2014 18:27:10 +0200 |
parents | 49a5a4be04a1 |
children | 4197fc428c7d |
rev | line source |
---|---|
3 | 1 /* |
2 | |
17744
d63878346099
maint: Update copyright notices for release.
John W. Eaton <jwe@octave.org>
parents:
15271
diff
changeset
|
3 Copyright (C) 1993-2013 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++) |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
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 *, | |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
92 const octave_idx_type&, const octave_idx_type&, double *pd, |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
93 const octave_idx_type& nrowpd) |
3 | 94 { |
4180 | 95 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
96 | |
1251 | 97 Matrix tmp_jac (neq, neq); |
3 | 98 |
1360 | 99 // NOTE: this won't work if LSODE passes copies of the state vector. |
100 // In that case we have to create a temporary vector object | |
101 // and copy. | |
102 | |
1251 | 103 tmp_jac = (*user_jac) (*tmp_x, time); |
3 | 104 |
5275 | 105 for (octave_idx_type j = 0; j < neq; j++) |
106 for (octave_idx_type i = 0; i < neq; i++) | |
15020
560317fd5977
maint: Cuddle open bracket used for indexing C++ arrays in source code.
Rik <rik@octave.org>
parents:
15018
diff
changeset
|
107 pd[nrowpd * j + i] = tmp_jac (i, j); |
3 | 108 |
4180 | 109 END_INTERRUPT_WITH_EXCEPTIONS; |
110 | |
3 | 111 return 0; |
112 } | |
113 | |
114 ColumnVector | |
1842 | 115 LSODE::do_integrate (double tout) |
3 | 116 { |
1945 | 117 ColumnVector retval; |
118 | |
5275 | 119 static octave_idx_type nn = 0; |
4049 | 120 |
121 if (! initialized || restart || ODEFunc::reset || LSODE_options::reset) | |
1945 | 122 { |
4049 | 123 integration_error = false; |
1945 | 124 |
4049 | 125 initialized = true; |
126 | |
127 istate = 1; | |
128 | |
5275 | 129 octave_idx_type n = size (); |
4049 | 130 |
131 nn = n; | |
3955 | 132 |
5275 | 133 octave_idx_type max_maxord = 0; |
4231 | 134 |
4049 | 135 if (integration_method () == "stiff") |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
136 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
137 max_maxord = 5; |
4231 | 138 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
139 if (jac) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
140 method_flag = 21; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
141 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
142 method_flag = 22; |
3955 | 143 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
144 liw = 20 + n; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
145 lrw = 22 + n * (9 + n); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
146 } |
4049 | 147 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
148 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
149 max_maxord = 12; |
4231 | 150 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
151 method_flag = 10; |
3955 | 152 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
153 liw = 20; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
154 lrw = 22 + 16 * n; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
155 } |
4049 | 156 |
4231 | 157 maxord = maximum_order (); |
158 | |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
159 iwork.resize (dim_vector (liw, 1)); |
5552 | 160 |
161 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
|
162 iwork(i) = 0; |
5552 | 163 |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
164 rwork.resize (dim_vector (lrw, 1)); |
5552 | 165 |
166 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
|
167 rwork(i) = 0; |
5552 | 168 |
4231 | 169 if (maxord >= 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
170 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
171 if (maxord > 0 && maxord <= max_maxord) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
172 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
173 iwork(4) = maxord; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
174 iopt = 1; |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
175 } |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
176 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
177 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
178 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
179 ("lsode: invalid value for maximum order"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
180 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
181 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
182 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
183 } |
4231 | 184 |
4049 | 185 if (stop_time_set) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
186 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
187 itask = 4; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
188 rwork(0) = stop_time; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
189 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
190 } |
4049 | 191 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
192 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
193 itask = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
194 } |
258 | 195 |
4049 | 196 restart = false; |
197 | |
198 // ODEFunc | |
3 | 199 |
4049 | 200 // NOTE: this won't work if LSODE passes copies of the state vector. |
201 // In that case we have to create a temporary vector object | |
202 // and copy. | |
3 | 203 |
4049 | 204 tmp_x = &x; |
205 | |
206 user_fun = function (); | |
207 user_jac = jacobian_function (); | |
208 | |
2343 | 209 ColumnVector xdot = (*user_fun) (x, t); |
210 | |
211 if (x.length () != xdot.length ()) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
212 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
213 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
214 ("lsode: inconsistent sizes for state and derivative vectors"); |
2343 | 215 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
216 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
217 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
218 } |
2343 | 219 |
4049 | 220 ODEFunc::reset = false; |
221 | |
222 // LSODE_options | |
223 | |
224 rel_tol = relative_tolerance (); | |
225 abs_tol = absolute_tolerance (); | |
226 | |
5275 | 227 octave_idx_type abs_tol_len = abs_tol.length (); |
4049 | 228 |
229 if (abs_tol_len == 1) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
230 itol = 1; |
4049 | 231 else if (abs_tol_len == n) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
232 itol = 2; |
4049 | 233 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
234 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
235 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
236 ("lsode: inconsistent sizes for state and absolute tolerance vectors"); |
4049 | 237 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
238 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
239 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
240 } |
2343 | 241 |
4049 | 242 double iss = initial_step_size (); |
243 if (iss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
244 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 rwork(4) = iss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
247 } |
4049 | 248 |
249 double maxss = maximum_step_size (); | |
250 if (maxss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
251 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
252 rwork(5) = maxss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
253 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
254 } |
4049 | 255 |
256 double minss = minimum_step_size (); | |
257 if (minss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
258 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
259 rwork(6) = minss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
260 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
261 } |
4049 | 262 |
5275 | 263 octave_idx_type sl = step_limit (); |
4049 | 264 if (sl > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
265 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
266 iwork(5) = sl; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
267 iopt = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
268 } |
4049 | 269 |
270 LSODE_options::reset = false; | |
3 | 271 } |
272 | |
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
|
273 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
|
274 |
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 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
|
276 |
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 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
|
278 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
|
279 |
4583 | 280 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
|
281 pabs_tol, itask, istate, iopt, prwork, lrw, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
282 piwork, liw, lsode_j, method_flag)); |
3 | 283 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
284 switch (istate) |
3 | 285 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
286 case 1: // prior to initial integration step. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
287 case 2: // lsode was successful. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
288 retval = x; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
289 t = tout; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
290 break; |
3996 | 291 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
292 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
|
293 case -2: // excess accuracy requested (tolerances too small). |
8807
401d54a83690
use 'invalid', not 'illegal'
John W. Eaton <jwe@octave.org>
parents:
7482
diff
changeset
|
294 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
|
295 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
|
296 case -5: // repeated convergence failures (perhaps bad Jacobian |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
297 // 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
|
298 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
|
299 // 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
|
300 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
|
301 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
302 break; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
303 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
304 default: |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
305 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
306 (*current_liboctave_error_handler) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
307 ("unrecognized value of istate (= %d) returned from lsode", |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
308 istate); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
309 break; |
3 | 310 } |
311 | |
1945 | 312 return retval; |
3 | 313 } |
314 | |
3959 | 315 std::string |
316 LSODE::error_message (void) const | |
317 { | |
318 std::string retval; | |
319 | |
5765 | 320 std::ostringstream buf; |
321 buf << t; | |
322 std::string t_curr = buf.str (); | |
4042 | 323 |
3959 | 324 switch (istate) |
325 { | |
326 case 1: | |
3996 | 327 retval = "prior to initial integration step"; |
3959 | 328 break; |
329 | |
330 case 2: | |
331 retval = "successful exit"; | |
332 break; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
333 |
3959 | 334 case 3: |
335 retval = "prior to continuation call with modified parameters"; | |
336 break; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
337 |
3996 | 338 case -1: |
4042 | 339 retval = std::string ("excess work on this call (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
340 + t_curr + "; perhaps wrong integration method)"; |
3996 | 341 break; |
342 | |
343 case -2: | |
344 retval = "excess accuracy requested (tolerances too small)"; | |
345 break; | |
346 | |
347 case -3: | |
348 retval = "invalid input detected (see printed message)"; | |
349 break; | |
350 | |
351 case -4: | |
4042 | 352 retval = std::string ("repeated error test failures (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
353 + t_curr + "; check all inputs)"; |
3996 | 354 break; |
355 | |
356 case -5: | |
4042 | 357 retval = std::string ("repeated convergence failures (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
358 + t_curr |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
359 + "; perhaps bad Jacobian supplied or wrong choice of integration method or tolerances)"; |
3996 | 360 break; |
361 | |
362 case -6: | |
4042 | 363 retval = std::string ("error weight became zero during problem. (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
364 + t_curr |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
365 + "; solution component i vanished, and atol or atol(i) == 0)"; |
3996 | 366 break; |
367 | |
368 case -13: | |
4042 | 369 retval = "return requested in user-supplied function (t = " |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
370 + t_curr + ")"; |
3996 | 371 break; |
372 | |
3959 | 373 default: |
374 retval = "unknown error state"; | |
375 break; | |
376 } | |
377 | |
378 return retval; | |
379 } | |
380 | |
3 | 381 Matrix |
1842 | 382 LSODE::do_integrate (const ColumnVector& tout) |
3 | 383 { |
384 Matrix retval; | |
4049 | 385 |
5275 | 386 octave_idx_type n_out = tout.capacity (); |
387 octave_idx_type n = size (); | |
3 | 388 |
389 if (n_out > 0 && n > 0) | |
390 { | |
391 retval.resize (n_out, n); | |
392 | |
5275 | 393 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
|
394 retval.elem (0, i) = x.elem (i); |
3 | 395 |
5275 | 396 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
|
397 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
398 ColumnVector x_next = do_integrate (tout.elem (j)); |
258 | 399 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
400 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
401 return retval; |
258 | 402 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
403 for (octave_idx_type i = 0; i < n; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
404 retval.elem (j, i) = x_next.elem (i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
405 } |
3 | 406 } |
407 | |
408 return retval; | |
409 } | |
410 | |
411 Matrix | |
3511 | 412 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit) |
3 | 413 { |
414 Matrix retval; | |
4049 | 415 |
5275 | 416 octave_idx_type n_out = tout.capacity (); |
417 octave_idx_type n = size (); | |
3 | 418 |
419 if (n_out > 0 && n > 0) | |
420 { | |
421 retval.resize (n_out, n); | |
422 | |
5275 | 423 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
|
424 retval.elem (0, i) = x.elem (i); |
3 | 425 |
5275 | 426 octave_idx_type n_crit = tcrit.capacity (); |
3 | 427 |
428 if (n_crit > 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
429 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
430 octave_idx_type i_crit = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
431 octave_idx_type i_out = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
432 double next_crit = tcrit.elem (0); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
433 double next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
434 while (i_out < n_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
435 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
436 bool do_restart = false; |
3 | 437 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
438 next_out = tout.elem (i_out); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
439 if (i_crit < n_crit) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
440 next_crit = tcrit.elem (i_crit); |
3 | 441 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
442 octave_idx_type save_output; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
443 double t_out; |
3 | 444 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
445 if (next_crit == next_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
446 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
447 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
448 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
449 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
450 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
451 i_crit++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
452 do_restart = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
453 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
454 else if (next_crit < next_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
455 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
456 if (i_crit < n_crit) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
457 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
458 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
459 t_out = next_crit; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
460 save_output = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
461 i_crit++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
462 do_restart = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
463 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
464 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
465 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
466 clear_stop_time (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
467 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
468 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
469 i_out++; |
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 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
472 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
473 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
474 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
475 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
476 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
477 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
478 } |
3 | 479 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
480 ColumnVector x_next = do_integrate (t_out); |
3 | 481 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
482 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
483 return retval; |
258 | 484 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
485 if (save_output) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
486 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
487 for (octave_idx_type i = 0; i < n; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
488 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
|
489 } |
3 | 490 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
491 if (do_restart) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
492 force_restart (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
493 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
494 } |
3 | 495 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
496 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
497 retval = do_integrate (tout); |
258 | 498 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
499 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
500 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
501 } |
3 | 502 } |
503 | |
504 return retval; | |
505 } |