Mercurial > octave
annotate libinterp/corefcn/daspk.cc @ 20725:7164a69472dd
eliminate more uses of error state
* daspk.cc, dasrt.cc, dassl.cc, lsode.cc, quad.cc: Eliminate uses of
error_state. Handle user-supplied functions consistently.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 19 Nov 2015 15:46:40 -0500 |
parents | 7b608fadc663 |
children | 200ae1d650b7 |
rev | line source |
---|---|
3912 | 1 /* |
2 | |
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17787
diff
changeset
|
3 Copyright (C) 1996-2015 John W. Eaton |
3912 | 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. | |
3912 | 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/>. | |
3912 | 20 |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
24 #include <config.h> | |
25 #endif | |
26 | |
27 #include <string> | |
28 | |
29 #include <iomanip> | |
30 #include <iostream> | |
31 | |
32 #include "DASPK.h" | |
33 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14854
diff
changeset
|
34 #include "defun.h" |
3912 | 35 #include "error.h" |
36 #include "gripes.h" | |
37 #include "oct-obj.h" | |
38 #include "ov-fcn.h" | |
5729 | 39 #include "ov-cell.h" |
3912 | 40 #include "pager.h" |
41 #include "unwind-prot.h" | |
42 #include "utils.h" | |
43 #include "variables.h" | |
44 | |
3998 | 45 #include "DASPK-opts.cc" |
46 | |
3912 | 47 // Global pointer for user defined function required by daspk. |
48 static octave_function *daspk_fcn; | |
49 | |
4140 | 50 // Global pointer for optional user defined jacobian function. |
51 static octave_function *daspk_jac; | |
52 | |
53 // Have we warned about imaginary values returned from user function? | |
54 static bool warned_fcn_imaginary = false; | |
55 static bool warned_jac_imaginary = false; | |
56 | |
3912 | 57 // Is this a recursive call? |
58 static int call_depth = 0; | |
59 | |
60 ColumnVector | |
61 daspk_user_function (const ColumnVector& x, const ColumnVector& xdot, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
62 double t, octave_idx_type& ires) |
3912 | 63 { |
64 ColumnVector retval; | |
65 | |
20229
5dfaaaae784f
Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
66 assert (x.numel () == xdot.numel ()); |
3912 | 67 |
68 octave_value_list args; | |
4628 | 69 |
3912 | 70 args(2) = t; |
4628 | 71 args(1) = xdot; |
72 args(0) = x; | |
3912 | 73 |
74 if (daspk_fcn) | |
75 { | |
20725
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
76 octave_value_list tmp; |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
77 |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
78 try |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
79 { |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
80 tmp = daspk_fcn->do_multi_index_op (1, args); |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
81 } |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
82 catch (const octave_execution_exception&) |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
83 { |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
84 gripe_user_supplied_eval ("daspk"); |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
85 throw; |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
86 } |
3912 | 87 |
88 int tlen = tmp.length (); | |
89 if (tlen > 0 && tmp(0).is_defined ()) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
90 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
91 if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
92 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
93 warning ("daspk: ignoring imaginary part returned from user-supplied function"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
94 warned_fcn_imaginary = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
95 } |
4140 | 96 |
20703
85e5efae848a
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20687
diff
changeset
|
97 retval = tmp(0).vector_value (); |
3912 | 98 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
99 if (tlen > 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
100 ires = tmp(1).int_value (); |
3912 | 101 |
20703
85e5efae848a
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20687
diff
changeset
|
102 if (retval.numel () == 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
103 gripe_user_supplied_eval ("daspk"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
104 } |
3912 | 105 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
106 gripe_user_supplied_eval ("daspk"); |
3912 | 107 } |
108 | |
109 return retval; | |
110 } | |
111 | |
4140 | 112 Matrix |
113 daspk_user_jacobian (const ColumnVector& x, const ColumnVector& xdot, | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
114 double t, double cj) |
4140 | 115 { |
116 Matrix retval; | |
117 | |
20229
5dfaaaae784f
Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
118 assert (x.numel () == xdot.numel ()); |
4140 | 119 |
120 octave_value_list args; | |
121 | |
122 args(3) = cj; | |
123 args(2) = t; | |
4628 | 124 args(1) = xdot; |
125 args(0) = x; | |
4140 | 126 |
127 if (daspk_jac) | |
128 { | |
20725
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
129 octave_value_list tmp; |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
130 |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
131 try |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
132 { |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
133 tmp = daspk_jac->do_multi_index_op (1, args); |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
134 } |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
135 catch (const octave_execution_exception&) |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
136 { |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
137 gripe_user_supplied_eval ("daspk"); |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
138 throw; |
7164a69472dd
eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents:
20711
diff
changeset
|
139 } |
4140 | 140 |
141 int tlen = tmp.length (); | |
142 if (tlen > 0 && tmp(0).is_defined ()) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
143 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
144 if (! warned_jac_imaginary && tmp(0).is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
145 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
146 warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
147 warned_jac_imaginary = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
148 } |
4140 | 149 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
150 retval = tmp(0).matrix_value (); |
4140 | 151 |
20703
85e5efae848a
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20687
diff
changeset
|
152 if (retval.numel () == 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
153 gripe_user_supplied_eval ("daspk"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
154 } |
4140 | 155 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
156 gripe_user_supplied_eval ("daspk"); |
4140 | 157 } |
158 | |
159 return retval; | |
160 } | |
161 | |
3912 | 162 #define DASPK_ABORT() \ |
10066
2cd940306a06
make unwind_protect frames local
Jaroslav Hajek <highegg@gmail.com>
parents:
9377
diff
changeset
|
163 return retval |
3912 | 164 |
165 #define DASPK_ABORT1(msg) \ | |
166 do \ | |
167 { \ | |
20428
7ac907da9fba
Use error() rather than ::error() unless explicitly required.
Rik <rik@octave.org>
parents:
20232
diff
changeset
|
168 error ("daspk: " msg); \ |
3912 | 169 DASPK_ABORT (); \ |
170 } \ | |
171 while (0) | |
172 | |
173 #define DASPK_ABORT2(fmt, arg) \ | |
174 do \ | |
175 { \ | |
20428
7ac907da9fba
Use error() rather than ::error() unless explicitly required.
Rik <rik@octave.org>
parents:
20232
diff
changeset
|
176 error ("daspk: " fmt, arg); \ |
3912 | 177 DASPK_ABORT (); \ |
178 } \ | |
179 while (0) | |
180 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14854
diff
changeset
|
181 DEFUN (daspk, args, nargout, |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
182 "-*- texinfo -*-\n\ |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14854
diff
changeset
|
183 @deftypefn {Built-in Function} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} daspk (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\ |
4115 | 184 Solve the set of differential-algebraic equations\n\ |
185 @tex\n\ | |
4852 | 186 $$ 0 = f (x, \\dot{x}, t) $$\n\ |
4115 | 187 with\n\ |
188 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\ | |
189 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
190 @ifnottex\n\ |
4115 | 191 \n\ |
192 @example\n\ | |
4698 | 193 0 = f (x, xdot, t)\n\ |
4115 | 194 @end example\n\ |
195 \n\ | |
10846
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
10840
diff
changeset
|
196 @noindent\n\ |
4115 | 197 with\n\ |
198 \n\ | |
199 @example\n\ | |
200 x(t_0) = x_0, xdot(t_0) = xdot_0\n\ | |
201 @end example\n\ | |
202 \n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
203 @end ifnottex\n\ |
4115 | 204 The solution is returned in the matrices @var{x} and @var{xdot},\n\ |
205 with each row in the result matrices corresponding to one of the\n\ | |
3912 | 206 elements in the vector @var{t}. The first element of @var{t}\n\ |
4115 | 207 should be @math{t_0} and correspond to the initial state of the\n\ |
208 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\ | |
209 row of the output @var{x} is @var{x_0} and the first row\n\ | |
210 of the output @var{xdot} is @var{xdot_0}.\n\ | |
3912 | 211 \n\ |
8787 | 212 The first argument, @var{fcn}, is a string, inline, or function handle\n\ |
213 that names the function @math{f} to call to compute the vector of\n\ | |
214 residuals for the set of equations. It must have the form\n\ | |
3912 | 215 \n\ |
216 @example\n\ | |
217 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ | |
218 @end example\n\ | |
219 \n\ | |
220 @noindent\n\ | |
4115 | 221 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ |
3912 | 222 scalar.\n\ |
223 \n\ | |
8787 | 224 If @var{fcn} is a two-element string array or a two-element cell array\n\ |
9067
8970b4b10e9f
Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
225 of strings, inline functions, or function handles, the first element names\n\ |
8787 | 226 the function @math{f} described above, and the second element names a\n\ |
227 function to compute the modified Jacobian\n\ | |
4115 | 228 @tex\n\ |
229 $$\n\ | |
230 J = {\\partial f \\over \\partial x}\n\ | |
231 + c {\\partial f \\over \\partial \\dot{x}}\n\ | |
232 $$\n\ | |
233 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
234 @ifnottex\n\ |
4144 | 235 \n\ |
236 @example\n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
237 @group\n\ |
4115 | 238 df df\n\ |
239 jac = -- + c ------\n\ | |
240 dx d xdot\n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
241 @end group\n\ |
4115 | 242 @end example\n\ |
10840 | 243 \n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
244 @end ifnottex\n\ |
4115 | 245 \n\ |
246 The modified Jacobian function must have the form\n\ | |
247 \n\ | |
248 @example\n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
249 @group\n\ |
4115 | 250 \n\ |
251 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\ | |
252 \n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
253 @end group\n\ |
4115 | 254 @end example\n\ |
255 \n\ | |
3912 | 256 The second and third arguments to @code{daspk} specify the initial\n\ |
257 condition of the states and their derivatives, and the fourth argument\n\ | |
4115 | 258 specifies a vector of output times at which the solution is desired,\n\ |
3912 | 259 including the time corresponding to the initial condition.\n\ |
260 \n\ | |
261 The set of initial states and derivatives are not strictly required to\n\ | |
4115 | 262 be consistent. If they are not consistent, you must use the\n\ |
263 @code{daspk_options} function to provide additional information so\n\ | |
264 that @code{daspk} can compute a consistent starting point.\n\ | |
3912 | 265 \n\ |
266 The fifth argument is optional, and may be used to specify a set of\n\ | |
267 times that the DAE solver should not integrate past. It is useful for\n\ | |
268 avoiding difficulties with singularities and points where there is a\n\ | |
269 discontinuity in the derivative.\n\ | |
3964 | 270 \n\ |
4115 | 271 After a successful computation, the value of @var{istate} will be\n\ |
10840 | 272 greater than zero (consistent with the Fortran version of @sc{daspk}).\n\ |
4115 | 273 \n\ |
274 If the computation is not successful, the value of @var{istate} will be\n\ | |
275 less than zero and @var{msg} will contain additional information.\n\ | |
276 \n\ | |
3964 | 277 You can use the function @code{daspk_options} to set optional\n\ |
278 parameters for @code{daspk}.\n\ | |
5642 | 279 @seealso{dassl}\n\ |
280 @end deftypefn") | |
3912 | 281 { |
282 octave_value_list retval; | |
283 | |
4140 | 284 warned_fcn_imaginary = false; |
285 warned_jac_imaginary = false; | |
286 | |
10066
2cd940306a06
make unwind_protect frames local
Jaroslav Hajek <highegg@gmail.com>
parents:
9377
diff
changeset
|
287 unwind_protect frame; |
3912 | 288 |
10066
2cd940306a06
make unwind_protect frames local
Jaroslav Hajek <highegg@gmail.com>
parents:
9377
diff
changeset
|
289 frame.protect_var (call_depth); |
3912 | 290 call_depth++; |
291 | |
292 if (call_depth > 1) | |
293 DASPK_ABORT1 ("invalid recursive call"); | |
294 | |
295 int nargin = args.length (); | |
296 | |
297 if (nargin > 3 && nargin < 6) | |
298 { | |
5729 | 299 std::string fcn_name, fname, jac_name, jname; |
4140 | 300 daspk_fcn = 0; |
301 daspk_jac = 0; | |
302 | |
303 octave_value f_arg = args(0); | |
304 | |
5729 | 305 if (f_arg.is_cell ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
306 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
307 Cell c = f_arg.cell_value (); |
20232
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20229
diff
changeset
|
308 if (c.numel () == 1) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
309 f_arg = c(0); |
20232
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
20229
diff
changeset
|
310 else if (c.numel () == 2) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
311 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
312 if (c(0).is_function_handle () || c(0).is_inline_function ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
313 daspk_fcn = c(0).function_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
314 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
315 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
316 fcn_name = unique_symbol_name ("__daspk_fcn__"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
317 fname = "function y = "; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
318 fname.append (fcn_name); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
319 fname.append (" (x, xdot, t) y = "); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
320 daspk_fcn = extract_function |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
321 (c(0), "daspk", fcn_name, fname, "; endfunction"); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
322 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
323 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
324 if (daspk_fcn) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
325 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
326 if (c(1).is_function_handle () || c(1).is_inline_function ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
327 daspk_jac = c(1).function_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
328 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
329 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
330 jac_name = unique_symbol_name ("__daspk_jac__"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
331 jname = "function jac = "; |
14854
5ae9f0f77635
maint: Use Octave coding conventions for coddling parenthis is DLD-FUNCTIONS directory
Rik <octave@nomad.inbox5.com>
parents:
14846
diff
changeset
|
332 jname.append (jac_name); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
333 jname.append (" (x, xdot, t, cj) jac = "); |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
334 daspk_jac = extract_function (c(1), "daspk", jac_name, |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
335 jname, "; endfunction"); |
4140 | 336 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
337 if (!daspk_jac) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
338 { |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
339 if (fcn_name.length ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
340 clear_function (fcn_name); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
341 daspk_fcn = 0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
342 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
343 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
344 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
345 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
346 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
347 DASPK_ABORT1 ("incorrect number of elements in cell array"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
348 } |
3912 | 349 |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
350 if (!daspk_fcn && ! f_arg.is_cell ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
351 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
352 if (f_arg.is_function_handle () || f_arg.is_inline_function ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
353 daspk_fcn = f_arg.function_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
354 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
355 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
356 switch (f_arg.rows ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
357 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
358 case 1: |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
359 do |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
360 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
361 fcn_name = unique_symbol_name ("__daspk_fcn__"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
362 fname = "function y = "; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
363 fname.append (fcn_name); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
364 fname.append (" (x, xdot, t) y = "); |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
365 daspk_fcn = extract_function (f_arg, "daspk", fcn_name, |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
366 fname, "; endfunction"); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
367 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
368 while (0); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
369 break; |
4140 | 370 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
371 case 2: |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
372 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
373 string_vector tmp = f_arg.all_strings (); |
5729 | 374 |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
375 fcn_name = unique_symbol_name ("__daspk_fcn__"); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
376 fname = "function y = "; |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
377 fname.append (fcn_name); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
378 fname.append (" (x, xdot, t) y = "); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
379 daspk_fcn = extract_function (tmp(0), "daspk", fcn_name, |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
380 fname, "; endfunction"); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
381 |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
382 if (daspk_fcn) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
383 { |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
384 jac_name = unique_symbol_name ("__daspk_jac__"); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
385 jname = "function jac = "; |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
386 jname.append (jac_name); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
387 jname.append (" (x, xdot, t, cj) jac = "); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
388 daspk_jac = extract_function (tmp(1), "daspk", |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
389 jac_name, jname, |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
390 "; endfunction"); |
4140 | 391 |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
392 if (!daspk_jac) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
393 { |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
394 if (fcn_name.length ()) |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
395 clear_function (fcn_name); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
396 daspk_fcn = 0; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
397 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
398 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
399 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
400 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
401 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
402 } |
4140 | 403 |
20703
85e5efae848a
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20687
diff
changeset
|
404 if (! daspk_fcn) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
405 DASPK_ABORT (); |
3912 | 406 |
20711
7b608fadc663
Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents:
20703
diff
changeset
|
407 ColumnVector state = args(1).xvector_value ("daspk: initial state X_0 must be a vector"); |
3912 | 408 |
20711
7b608fadc663
Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents:
20703
diff
changeset
|
409 ColumnVector deriv = args(2).xvector_value ("daspk: initial derivatives XDOT_0 must be a vector"); |
3912 | 410 |
20711
7b608fadc663
Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents:
20703
diff
changeset
|
411 ColumnVector out_times = args(3).xvector_value ("daspk: output time variable T must be a vector"); |
3912 | 412 |
413 ColumnVector crit_times; | |
414 int crit_times_set = 0; | |
415 if (nargin > 4) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
416 { |
20711
7b608fadc663
Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents:
20703
diff
changeset
|
417 crit_times = args(4).xvector_value ("daspk: list of critical times T_CRIT must be a vector"); |
3912 | 418 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
419 crit_times_set = 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
420 } |
3912 | 421 |
20229
5dfaaaae784f
Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
422 if (state.numel () != deriv.numel ()) |
20711
7b608fadc663
Make error messages more specific about the variable and problem encountered.
Rik <rik@octave.org>
parents:
20703
diff
changeset
|
423 DASPK_ABORT1 ("X_0 and XDOT_0 must have the same size"); |
3912 | 424 |
425 double tzero = out_times (0); | |
426 | |
427 DAEFunc func (daspk_user_function); | |
4140 | 428 if (daspk_jac) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
429 func.set_jacobian_function (daspk_user_jacobian); |
4140 | 430 |
3912 | 431 DASPK dae (state, deriv, tzero, func); |
4122 | 432 dae.set_options (daspk_opts); |
3912 | 433 |
434 Matrix output; | |
435 Matrix deriv_output; | |
436 | |
437 if (crit_times_set) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
438 output = dae.integrate (out_times, deriv_output, crit_times); |
3912 | 439 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
440 output = dae.integrate (out_times, deriv_output); |
3912 | 441 |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
442 if (fcn_name.length ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
443 clear_function (fcn_name); |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
444 if (jac_name.length ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
445 clear_function (jac_name); |
5729 | 446 |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
447 std::string msg = dae.error_message (); |
3997 | 448 |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
449 retval(3) = msg; |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
450 retval(2) = static_cast<double> (dae.integration_state ()); |
3912 | 451 |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
452 if (dae.integration_ok ()) |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
453 { |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
454 retval(1) = deriv_output; |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
455 retval(0) = output; |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
456 } |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
457 else |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
458 { |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
459 retval(1) = Matrix (); |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
460 retval(0) = Matrix (); |
3997 | 461 |
20556
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
462 if (nargout < 3) |
4bed806ee3d4
eliminate more simple uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20428
diff
changeset
|
463 error ("daspk: %s", msg.c_str ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
10066
diff
changeset
|
464 } |
3912 | 465 } |
466 else | |
5823 | 467 print_usage (); |
3912 | 468 |
469 return retval; | |
470 } |