Mercurial > hg > octave-nkf
annotate liboctave/numeric/DASRT.cc @ 20685:7fa1970a655d
pkg.m: drop check of nargout value, the interpreter already does that.
* scripts/pkg/pkg.m: the interpreter already checks if there was any variable
that got no value assigned, there's no need to make the code more
complicated to cover that. Also, there's no point in calling describe()
with different nargout since it doesn't check nargout.
author | Carnë Draug <carandraug@octave.org> |
---|---|
date | Thu, 03 Sep 2015 16:21:08 +0100 |
parents | a9574e3c6e9e |
children |
rev | line source |
---|---|
3990 | 1 /* |
2 | |
19898
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17769
diff
changeset
|
3 Copyright (C) 2002-2015 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 |
11495 | 37 typedef octave_idx_type (*dasrt_fcn_ptr) (const double&, const double*, |
38 const double*, double*, | |
39 octave_idx_type&, double*, | |
40 octave_idx_type*); | |
3991 | 41 |
11495 | 42 typedef octave_idx_type (*dasrt_jac_ptr) (const double&, const double*, |
43 const double*, double*, | |
44 const double&, double*, | |
45 octave_idx_type*); | |
3991 | 46 |
11495 | 47 typedef octave_idx_type (*dasrt_constr_ptr) (const octave_idx_type&, |
48 const double&, const double*, | |
49 const octave_idx_type&, | |
50 double*, double*, | |
51 octave_idx_type*); | |
3991 | 52 |
3990 | 53 extern "C" |
4552 | 54 { |
55 F77_RET_T | |
11495 | 56 F77_FUNC (ddasrt, DDASRT) (dasrt_fcn_ptr, const octave_idx_type&, |
57 double&, double*, double*, const double&, | |
58 octave_idx_type*, const double*, | |
59 const double*, octave_idx_type&, double*, | |
60 const octave_idx_type&, octave_idx_type*, | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
61 const octave_idx_type&, double*, |
11495 | 62 octave_idx_type*, dasrt_jac_ptr, |
63 dasrt_constr_ptr, const octave_idx_type&, | |
64 octave_idx_type*); | |
4552 | 65 } |
3990 | 66 |
67 static DAEFunc::DAERHSFunc user_fsub; | |
68 static DAEFunc::DAEJacFunc user_jsub; | |
69 static DAERTFunc::DAERTConstrFunc user_csub; | |
3993 | 70 |
5275 | 71 static octave_idx_type nn; |
3990 | 72 |
5275 | 73 static octave_idx_type |
3991 | 74 ddasrt_f (const double& t, const double *state, const double *deriv, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
75 double *delta, octave_idx_type& ires, double *, octave_idx_type *) |
3990 | 76 { |
4180 | 77 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
78 | |
3990 | 79 ColumnVector tmp_state (nn); |
4133 | 80 ColumnVector tmp_deriv (nn); |
3990 | 81 |
5275 | 82 for (octave_idx_type i = 0; i < nn; i++) |
4133 | 83 { |
84 tmp_state(i) = state[i]; | |
85 tmp_deriv(i) = deriv[i]; | |
86 } | |
3990 | 87 |
3994 | 88 ColumnVector tmp_fval = (*user_fsub) (tmp_state, tmp_deriv, t, ires); |
3990 | 89 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20439
diff
changeset
|
90 if (tmp_fval.numel () == 0) |
3990 | 91 ires = -2; |
92 else | |
93 { | |
5275 | 94 for (octave_idx_type i = 0; i < nn; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
95 delta[i] = tmp_fval(i); |
3990 | 96 } |
97 | |
4180 | 98 END_INTERRUPT_WITH_EXCEPTIONS; |
99 | |
3990 | 100 return 0; |
101 } | |
102 | |
5275 | 103 octave_idx_type |
3991 | 104 ddasrt_j (const double& time, const double *state, const double *deriv, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
105 double *pd, const double& cj, double *, octave_idx_type *) |
3990 | 106 { |
4180 | 107 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
108 | |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
109 // FIXME: would be nice to avoid copying the data. |
3990 | 110 |
3991 | 111 ColumnVector tmp_state (nn); |
112 ColumnVector tmp_deriv (nn); | |
3990 | 113 |
5275 | 114 for (octave_idx_type i = 0; i < nn; i++) |
3991 | 115 { |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
116 tmp_deriv.elem (i) = deriv[i]; |
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
117 tmp_state.elem (i) = state[i]; |
3991 | 118 } |
3990 | 119 |
3994 | 120 Matrix tmp_pd = (*user_jsub) (tmp_state, tmp_deriv, time, cj); |
3990 | 121 |
5275 | 122 for (octave_idx_type j = 0; j < nn; j++) |
123 for (octave_idx_type i = 0; i < nn; i++) | |
15020
560317fd5977
maint: Cuddle open bracket used for indexing C++ arrays in source code.
Rik <rik@octave.org>
parents:
15018
diff
changeset
|
124 pd[nn * j + i] = tmp_pd.elem (i, j); |
3990 | 125 |
4180 | 126 END_INTERRUPT_WITH_EXCEPTIONS; |
127 | |
3990 | 128 return 0; |
129 } | |
130 | |
5275 | 131 static octave_idx_type |
132 ddasrt_g (const octave_idx_type& neq, const double& t, const double *state, | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
133 const octave_idx_type& ng, double *gout, double *, octave_idx_type *) |
3990 | 134 { |
4180 | 135 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
136 | |
5275 | 137 octave_idx_type n = neq; |
3990 | 138 |
139 ColumnVector tmp_state (n); | |
5275 | 140 for (octave_idx_type i = 0; i < n; i++) |
3990 | 141 tmp_state(i) = state[i]; |
142 | |
3994 | 143 ColumnVector tmp_fval = (*user_csub) (tmp_state, t); |
3990 | 144 |
5275 | 145 for (octave_idx_type i = 0; i < ng; i++) |
3990 | 146 gout[i] = tmp_fval(i); |
147 | |
4180 | 148 END_INTERRUPT_WITH_EXCEPTIONS; |
149 | |
3990 | 150 return 0; |
151 } | |
152 | |
153 void | |
154 DASRT::integrate (double tout) | |
155 { | |
156 DASRT_result retval; | |
157 | |
4049 | 158 // I suppose this is the safe thing to do. If this is the first |
159 // call, or if anything about the problem has changed, we should | |
160 // start completely fresh. | |
161 | |
162 if (! initialized || restart | |
163 || DAEFunc::reset || DAERTFunc::reset || DASRT_options::reset) | |
3990 | 164 { |
165 integration_error = false; | |
166 | |
4049 | 167 initialized = true; |
168 | |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
169 info.resize (dim_vector (15, 1)); |
4049 | 170 |
5275 | 171 for (octave_idx_type i = 0; i < 15; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
172 info(i) = 0; |
4049 | 173 |
5275 | 174 octave_idx_type n = size (); |
4049 | 175 |
176 nn = n; | |
3990 | 177 |
4133 | 178 // DAERTFunc |
179 | |
180 user_csub = DAERTFunc::constraint_function (); | |
181 | |
182 if (user_csub) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
183 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
184 ColumnVector tmp = (*user_csub) (x, t); |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20439
diff
changeset
|
185 ng = tmp.numel (); |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
186 } |
4133 | 187 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
188 ng = 0; |
4133 | 189 |
5275 | 190 octave_idx_type maxord = maximum_order (); |
4133 | 191 if (maxord >= 0) |
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 if (maxord > 0 && maxord < 6) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
194 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
195 info(8) = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
196 iwork(2) = maxord; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
197 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
198 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
199 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
200 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
201 ("dassl: invalid value for maximum order"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
202 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
203 return; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
204 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
205 } |
4133 | 206 |
4428 | 207 liw = 21 + n; |
4133 | 208 lrw = 50 + 9*n + n*n + 3*ng; |
4049 | 209 |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
210 iwork.resize (dim_vector (liw, 1)); |
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
211 rwork.resize (dim_vector (lrw, 1)); |
4049 | 212 |
213 info(0) = 0; | |
214 | |
215 if (stop_time_set) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
216 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
217 info(3) = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
218 rwork(0) = stop_time; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
219 } |
3990 | 220 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
221 info(3) = 0; |
3990 | 222 |
4049 | 223 restart = false; |
3994 | 224 |
4049 | 225 // DAEFunc |
226 | |
227 user_fsub = DAEFunc::function (); | |
228 user_jsub = DAEFunc::jacobian_function (); | |
229 | |
230 if (user_fsub) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
231 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
232 octave_idx_type ires = 0; |
3990 | 233 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
234 ColumnVector fval = (*user_fsub) (x, xdot, t, ires); |
3990 | 235 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20439
diff
changeset
|
236 if (fval.numel () != x.numel ()) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
237 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
238 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
239 ("dasrt: inconsistent sizes for state and residual vectors"); |
3990 | 240 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
241 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
242 return; |
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 } |
4049 | 245 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
247 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
248 ("dasrt: no user supplied RHS subroutine!"); |
3990 | 249 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
250 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
251 return; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
252 } |
4049 | 253 |
254 info(4) = user_jsub ? 1 : 0; | |
255 | |
256 DAEFunc::reset = false; | |
257 | |
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
258 jroot.resize (dim_vector (ng, 1), 1); |
4049 | 259 |
260 DAERTFunc::reset = false; | |
261 | |
262 // DASRT_options | |
263 | |
4050 | 264 double mss = maximum_step_size (); |
265 if (mss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
266 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
267 rwork(1) = mss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
268 info(6) = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
269 } |
4049 | 270 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
271 info(6) = 0; |
3990 | 272 |
4050 | 273 double iss = initial_step_size (); |
274 if (iss >= 0.0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
275 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
276 rwork(2) = iss; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
277 info(7) = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
278 } |
4049 | 279 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
280 info(7) = 0; |
4049 | 281 |
282 if (step_limit () >= 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
283 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
284 info(11) = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
285 iwork(20) = step_limit (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
286 } |
4049 | 287 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
288 info(11) = 0; |
3990 | 289 |
290 abs_tol = absolute_tolerance (); | |
291 rel_tol = relative_tolerance (); | |
292 | |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20439
diff
changeset
|
293 octave_idx_type abs_tol_len = abs_tol.numel (); |
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20439
diff
changeset
|
294 octave_idx_type rel_tol_len = rel_tol.numel (); |
3998 | 295 |
296 if (abs_tol_len == 1 && rel_tol_len == 1) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
297 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
298 info.elem (1) = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
299 } |
3998 | 300 else if (abs_tol_len == n && rel_tol_len == n) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
301 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
302 info.elem (1) = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
303 } |
3998 | 304 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
305 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
306 (*current_liboctave_error_handler) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
307 ("dasrt: inconsistent sizes for tolerance arrays"); |
3998 | 308 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
309 integration_error = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
310 return; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
311 } |
3998 | 312 |
4049 | 313 DASRT_options::reset = false; |
3990 | 314 } |
315 | |
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
|
316 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
|
317 double *pxdot = xdot.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
|
318 |
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
|
319 octave_idx_type *pinfo = info.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
|
320 |
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
|
321 double *prel_tol = rel_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
|
322 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
|
323 |
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
|
324 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
|
325 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
|
326 |
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
|
327 octave_idx_type *pjroot = jroot.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
|
328 |
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
|
329 double *dummy = 0; |
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
|
330 octave_idx_type *idummy = 0; |
3990 | 331 |
4575 | 332 F77_XFCN (ddasrt, DDASRT, (ddasrt_f, nn, t, px, pxdot, tout, pinfo, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
333 prel_tol, pabs_tol, istate, prwork, lrw, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
334 piwork, liw, dummy, idummy, ddasrt_j, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
335 ddasrt_g, ng, pjroot)); |
3990 | 336 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
337 switch (istate) |
3990 | 338 { |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
339 case 1: // A step was successfully taken in intermediate-output |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
340 // mode. The code has not yet reached TOUT. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
341 case 2: // The integration to TOUT was successfully completed |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
342 // (T=TOUT) by stepping exactly to TOUT. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
343 case 3: // The integration to TOUT was successfully completed |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
344 // (T=TOUT) by stepping past TOUT. Y(*) is obtained by |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
345 // interpolation. YPRIME(*) is obtained by interpolation. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
346 t = tout; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
347 break; |
3990 | 348 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
349 case 4: // The integration was successfully completed |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
350 // by finding one or more roots of G at T. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
351 break; |
3990 | 352 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
353 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
|
354 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
|
355 case -3: // The local error test cannot be satisfied because you |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
356 // specified a zero component in ATOL and the |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
357 // corresponding computed solution component is zero. |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
358 // Thus, a pure relative error test is impossible for |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
359 // this component. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
360 case -6: // DDASRT had repeated error test failures on the last |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
361 // attempted step. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
362 case -7: // The corrector could not converge. |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
363 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
|
364 case -9: // The corrector could not converge. There were repeated |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
365 // error test failures in this step. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
366 case -10: // The corrector could not converge because IRES was |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
367 // equal to minus one. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
368 case -11: // IRES equal to -2 was encountered and control is being |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
369 // returned to the calling program. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
370 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
|
371 case -33: // The code has encountered trouble from which it cannot |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
372 // recover. A message is printed explaining the trouble |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
373 // and control is returned to the calling program. For |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
374 // example, this occurs when invalid input is detected. |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
375 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
376 break; |
3996 | 377 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
378 default: |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
379 integration_error = true; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
380 (*current_liboctave_error_handler) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
381 ("unrecognized value of istate (= %d) returned from ddasrt", |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
382 istate); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
383 break; |
3990 | 384 } |
385 } | |
386 | |
387 DASRT_result | |
388 DASRT::integrate (const ColumnVector& tout) | |
389 { | |
390 DASRT_result retval; | |
391 | |
392 Matrix x_out; | |
393 Matrix xdot_out; | |
3994 | 394 ColumnVector t_out = tout; |
3990 | 395 |
20439
5dfaaaae784f
Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
396 octave_idx_type n_out = tout.numel (); |
5275 | 397 octave_idx_type n = size (); |
3990 | 398 |
399 if (n_out > 0 && n > 0) | |
400 { | |
401 x_out.resize (n_out, n); | |
402 xdot_out.resize (n_out, n); | |
403 | |
5275 | 404 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
|
405 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
406 x_out(0,i) = x(i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
407 xdot_out(0,i) = xdot(i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
408 } |
3994 | 409 |
5275 | 410 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
|
411 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
412 integrate (tout(j)); |
3994 | 413 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
414 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
415 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
416 retval = DASRT_result (x_out, xdot_out, t_out); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
417 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
418 } |
3992 | 419 |
3997 | 420 if (istate == 4) |
3990 | 421 t_out(j) = t; |
422 else | |
423 t_out(j) = tout(j); | |
424 | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
425 for (octave_idx_type i = 0; i < n; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
426 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
427 x_out(j,i) = x(i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
428 xdot_out(j,i) = xdot(i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
429 } |
3992 | 430 |
3997 | 431 if (istate == 4) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
432 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
433 x_out.resize (j+1, n); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
434 xdot_out.resize (j+1, n); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
435 t_out.resize (j+1); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
436 break; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
437 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
438 } |
3990 | 439 } |
440 | |
441 retval = DASRT_result (x_out, xdot_out, t_out); | |
442 | |
443 return retval; | |
444 } | |
445 | |
446 DASRT_result | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11574
diff
changeset
|
447 DASRT::integrate (const ColumnVector& tout, const ColumnVector& tcrit) |
3990 | 448 { |
449 DASRT_result retval; | |
450 | |
451 Matrix x_out; | |
452 Matrix xdot_out; | |
3994 | 453 ColumnVector t_outs = tout; |
3990 | 454 |
20439
5dfaaaae784f
Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
455 octave_idx_type n_out = tout.numel (); |
5275 | 456 octave_idx_type n = size (); |
3990 | 457 |
458 if (n_out > 0 && n > 0) | |
459 { | |
460 x_out.resize (n_out, n); | |
461 xdot_out.resize (n_out, n); | |
462 | |
20439
5dfaaaae784f
Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
463 octave_idx_type n_crit = tcrit.numel (); |
3990 | 464 |
465 if (n_crit > 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
466 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
467 octave_idx_type i_crit = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
468 octave_idx_type i_out = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
469 double next_crit = tcrit(0); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
470 double next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
471 while (i_out < n_out) |
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 bool do_restart = false; |
3990 | 474 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
475 next_out = tout(i_out); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
476 if (i_crit < n_crit) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
477 next_crit = tcrit(i_crit); |
3990 | 478 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
479 octave_idx_type save_output; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
480 double t_out; |
3990 | 481 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
482 if (next_crit == next_out) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
483 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
484 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
485 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
486 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
487 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
488 i_crit++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
489 do_restart = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
490 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
491 else if (next_crit < next_out) |
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 if (i_crit < n_crit) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
494 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
495 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
496 t_out = next_crit; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
497 save_output = 0; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
498 i_crit++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
499 do_restart = true; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
500 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
501 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
502 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
503 clear_stop_time (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
504 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
505 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
506 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
507 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
508 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
509 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
510 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
511 set_stop_time (next_crit); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
512 t_out = next_out; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
513 save_output = 1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
514 i_out++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
515 } |
3990 | 516 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
517 integrate (t_out); |
3990 | 518 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
519 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
520 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
521 retval = DASRT_result (x_out, xdot_out, t_outs); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
522 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
523 } |
3992 | 524 |
3997 | 525 if (istate == 4) |
3990 | 526 t_out = t; |
527 | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
528 if (save_output) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
529 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
530 for (octave_idx_type i = 0; i < n; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
531 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
532 x_out(i_out-1,i) = x(i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
533 xdot_out(i_out-1,i) = xdot(i); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
534 } |
3994 | 535 |
3990 | 536 t_outs(i_out-1) = t_out; |
3994 | 537 |
3997 | 538 if (istate == 4) |
3990 | 539 { |
540 x_out.resize (i_out, n); | |
541 xdot_out.resize (i_out, n); | |
542 t_outs.resize (i_out); | |
543 i_out = n_out; | |
544 } | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
545 } |
3990 | 546 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
547 if (do_restart) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
548 force_restart (); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
549 } |
3990 | 550 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
551 retval = DASRT_result (x_out, xdot_out, t_outs); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
552 } |
3990 | 553 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
554 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
555 retval = integrate (tout); |
3990 | 556 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
557 if (integration_error) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
558 return retval; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
559 } |
3990 | 560 } |
561 | |
562 return retval; | |
563 } | |
564 | |
3995 | 565 std::string |
566 DASRT::error_message (void) const | |
567 { | |
568 std::string retval; | |
569 | |
5765 | 570 std::ostringstream buf; |
571 buf << t; | |
572 std::string t_curr = buf.str (); | |
4043 | 573 |
3997 | 574 switch (istate) |
3995 | 575 { |
3996 | 576 case 1: |
577 retval = "a step was successfully taken in intermediate-output mode."; | |
578 break; | |
579 | |
580 case 2: | |
581 retval = "integration completed by stepping exactly to TOUT"; | |
582 break; | |
583 | |
584 case 3: | |
585 retval = "integration to tout completed by stepping past TOUT"; | |
586 break; | |
587 | |
588 case 4: | |
589 retval = "integration completed by finding one or more roots of G at T"; | |
590 break; | |
591 | |
592 case -1: | |
4043 | 593 retval = std::string ("a large amount of work has been expended (t =") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
594 + t_curr + ")"; |
3996 | 595 break; |
596 | |
597 case -2: | |
598 retval = "the error tolerances are too stringent"; | |
599 break; | |
600 | |
601 case -3: | |
4043 | 602 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
|
603 + t_curr |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
604 + "; solution component i vanished, and atol or atol(i) == 0)"; |
3996 | 605 break; |
606 | |
607 case -6: | |
4043 | 608 retval = std::string ("repeated error test failures on the last attempted step (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
609 + t_curr + ")"; |
3996 | 610 break; |
611 | |
612 case -7: | |
4043 | 613 retval = std::string ("the corrector could not converge (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
614 + t_curr + ")"; |
3996 | 615 break; |
616 | |
617 case -8: | |
4043 | 618 retval = std::string ("the matrix of partial derivatives is singular (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
619 + t_curr + ")"; |
3996 | 620 break; |
621 | |
622 case -9: | |
4043 | 623 retval = std::string ("the corrector could not converge (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
624 + t_curr + "; repeated test failures)"; |
3996 | 625 break; |
626 | |
627 case -10: | |
4043 | 628 retval = std::string ("corrector could not converge because IRES was -1 (t = ") |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
629 + t_curr + ")"; |
3996 | 630 break; |
631 | |
632 case -11: | |
4043 | 633 retval = std::string ("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
|
634 + t_curr + ")"; |
3996 | 635 break; |
636 | |
637 case -12: | |
638 retval = "failed to compute consistent initial conditions"; | |
639 break; | |
640 | |
641 case -33: | |
642 retval = "unrecoverable error (see printed message)"; | |
643 break; | |
644 | |
3995 | 645 default: |
646 retval = "unknown error state"; | |
647 break; | |
648 } | |
649 | |
650 return retval; | |
651 } |