Mercurial > octave-nkf

comparison src/corefcn/daspk.cc @ 15039:e753177cde93

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory * __contourc__.cc, __dispatch__.cc, __lin_interpn__.cc, __pchip_deriv__.cc, __qp__.cc, balance.cc, besselj.cc, betainc.cc, bsxfun.cc, cellfun.cc, colloc.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, det.cc, dlmread.cc, dot.cc, eig.cc, fft.cc, fft2.cc, fftn.cc, filter.cc, find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc, hess.cc, hex2num.cc, inv.cc, kron.cc, lookup.cc, lsode.cc, lu.cc, luinc.cc, matrix_type.cc, max.cc, md5sum.cc, mgorth.cc, nproc.cc, pinv.cc, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, spparms.cc, sqrtm.cc, str2double.cc, strfind.cc, sub2ind.cc, svd.cc, syl.cc, time.cc, tril.cc, typecast.cc: Move functions from DLD-FUNCTIONS/ to corefcn/ directory. Include "defun.h", not "defun-dld.h". Change docstring to refer to these as "Built-in Functions". * build-aux/mk-opts.pl: Generate options code with '#include "defun.h"'. Change option docstrings to refer to these as "Built-in Functions". * corefcn/module.mk: List of functions to build in corefcn/ dir. * DLD-FUNCTIONS/config-module.awk: Update to new build system. * DLD-FUNCTIONS/module-files: Remove functions which are now in corefcn/ directory. * src/Makefile.am: Update to build "convenience library" in corefcn/. Octave program now links against all other libraries + corefcn libary. * src/find-defun-files.sh: Strip $srcdir from filename. * src/link-deps.mk: Add REGEX and FFTW link dependencies for liboctinterp. * type.m, which.m: Change failing tests to use 'amd', still a dynamic function, rather than 'dot', which isn't.
author Rik <rik@octave.org>
date Fri, 27 Jul 2012 15:35:00 -0700
parents src/DLD-FUNCTIONS/daspk.cc@5ae9f0f77635
children
comparison
equal deleted inserted replaced
15038:ab18578c2ade 15039:e753177cde93
1 /*
2
3 Copyright (C) 1996-2012 John W. Eaton
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
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
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
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
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
34 #include "defun.h"
35 #include "error.h"
36 #include "gripes.h"
37 #include "oct-obj.h"
38 #include "ov-fcn.h"
39 #include "ov-cell.h"
40 #include "pager.h"
41 #include "unwind-prot.h"
42 #include "utils.h"
43 #include "variables.h"
44
45 #include "DASPK-opts.cc"
46
47 // Global pointer for user defined function required by daspk.
48 static octave_function *daspk_fcn;
49
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
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,
62 double t, octave_idx_type& ires)
63 {
64 ColumnVector retval;
65
66 assert (x.capacity () == xdot.capacity ());
67
68 octave_value_list args;
69
70 args(2) = t;
71 args(1) = xdot;
72 args(0) = x;
73
74 if (daspk_fcn)
75 {
76 octave_value_list tmp = daspk_fcn->do_multi_index_op (1, args);
77
78 if (error_state)
79 {
80 gripe_user_supplied_eval ("daspk");
81 return retval;
82 }
83
84 int tlen = tmp.length ();
85 if (tlen > 0 && tmp(0).is_defined ())
86 {
87 if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
88 {
89 warning ("daspk: ignoring imaginary part returned from user-supplied function");
90 warned_fcn_imaginary = true;
91 }
92
93 retval = ColumnVector (tmp(0).vector_value ());
94
95 if (tlen > 1)
96 ires = tmp(1).int_value ();
97
98 if (error_state || retval.length () == 0)
99 gripe_user_supplied_eval ("daspk");
100 }
101 else
102 gripe_user_supplied_eval ("daspk");
103 }
104
105 return retval;
106 }
107
108 Matrix
109 daspk_user_jacobian (const ColumnVector& x, const ColumnVector& xdot,
110 double t, double cj)
111 {
112 Matrix retval;
113
114 assert (x.capacity () == xdot.capacity ());
115
116 octave_value_list args;
117
118 args(3) = cj;
119 args(2) = t;
120 args(1) = xdot;
121 args(0) = x;
122
123 if (daspk_jac)
124 {
125 octave_value_list tmp = daspk_jac->do_multi_index_op (1, args);
126
127 if (error_state)
128 {
129 gripe_user_supplied_eval ("daspk");
130 return retval;
131 }
132
133 int tlen = tmp.length ();
134 if (tlen > 0 && tmp(0).is_defined ())
135 {
136 if (! warned_jac_imaginary && tmp(0).is_complex_type ())
137 {
138 warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function");
139 warned_jac_imaginary = true;
140 }
141
142 retval = tmp(0).matrix_value ();
143
144 if (error_state || retval.length () == 0)
145 gripe_user_supplied_eval ("daspk");
146 }
147 else
148 gripe_user_supplied_eval ("daspk");
149 }
150
151 return retval;
152 }
153
154 #define DASPK_ABORT() \
155 return retval
156
157 #define DASPK_ABORT1(msg) \
158 do \
159 { \
160 ::error ("daspk: " msg); \
161 DASPK_ABORT (); \
162 } \
163 while (0)
164
165 #define DASPK_ABORT2(fmt, arg) \
166 do \
167 { \
168 ::error ("daspk: " fmt, arg); \
169 DASPK_ABORT (); \
170 } \
171 while (0)
172
173 DEFUN (daspk, args, nargout,
174 "-*- texinfo -*-\n\
175 @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\
176 Solve the set of differential-algebraic equations\n\
177 @tex\n\
178 $$ 0 = f (x, \\dot{x}, t) $$\n\
179 with\n\
180 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\
181 @end tex\n\
182 @ifnottex\n\
183 \n\
184 @example\n\
185 0 = f (x, xdot, t)\n\
186 @end example\n\
187 \n\
188 @noindent\n\
189 with\n\
190 \n\
191 @example\n\
192 x(t_0) = x_0, xdot(t_0) = xdot_0\n\
193 @end example\n\
194 \n\
195 @end ifnottex\n\
196 The solution is returned in the matrices @var{x} and @var{xdot},\n\
197 with each row in the result matrices corresponding to one of the\n\
198 elements in the vector @var{t}. The first element of @var{t}\n\
199 should be @math{t_0} and correspond to the initial state of the\n\
200 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\
201 row of the output @var{x} is @var{x_0} and the first row\n\
202 of the output @var{xdot} is @var{xdot_0}.\n\
203 \n\
204 The first argument, @var{fcn}, is a string, inline, or function handle\n\
205 that names the function @math{f} to call to compute the vector of\n\
206 residuals for the set of equations. It must have the form\n\
207 \n\
208 @example\n\
209 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
210 @end example\n\
211 \n\
212 @noindent\n\
213 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
214 scalar.\n\
215 \n\
216 If @var{fcn} is a two-element string array or a two-element cell array\n\
217 of strings, inline functions, or function handles, the first element names\n\
218 the function @math{f} described above, and the second element names a\n\
219 function to compute the modified Jacobian\n\
220 @tex\n\
221 $$\n\
222 J = {\\partial f \\over \\partial x}\n\
223 + c {\\partial f \\over \\partial \\dot{x}}\n\
224 $$\n\
225 @end tex\n\
226 @ifnottex\n\
227 \n\
228 @example\n\
229 @group\n\
230 df df\n\
231 jac = -- + c ------\n\
232 dx d xdot\n\
233 @end group\n\
234 @end example\n\
235 \n\
236 @end ifnottex\n\
237 \n\
238 The modified Jacobian function must have the form\n\
239 \n\
240 @example\n\
241 @group\n\
242 \n\
243 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
244 \n\
245 @end group\n\
246 @end example\n\
247 \n\
248 The second and third arguments to @code{daspk} specify the initial\n\
249 condition of the states and their derivatives, and the fourth argument\n\
250 specifies a vector of output times at which the solution is desired,\n\
251 including the time corresponding to the initial condition.\n\
252 \n\
253 The set of initial states and derivatives are not strictly required to\n\
254 be consistent. If they are not consistent, you must use the\n\
255 @code{daspk_options} function to provide additional information so\n\
256 that @code{daspk} can compute a consistent starting point.\n\
257 \n\
258 The fifth argument is optional, and may be used to specify a set of\n\
259 times that the DAE solver should not integrate past. It is useful for\n\
260 avoiding difficulties with singularities and points where there is a\n\
261 discontinuity in the derivative.\n\
262 \n\
263 After a successful computation, the value of @var{istate} will be\n\
264 greater than zero (consistent with the Fortran version of @sc{daspk}).\n\
265 \n\
266 If the computation is not successful, the value of @var{istate} will be\n\
267 less than zero and @var{msg} will contain additional information.\n\
268 \n\
269 You can use the function @code{daspk_options} to set optional\n\
270 parameters for @code{daspk}.\n\
271 @seealso{dassl}\n\
272 @end deftypefn")
273 {
274 octave_value_list retval;
275
276 warned_fcn_imaginary = false;
277 warned_jac_imaginary = false;
278
279 unwind_protect frame;
280
281 frame.protect_var (call_depth);
282 call_depth++;
283
284 if (call_depth > 1)
285 DASPK_ABORT1 ("invalid recursive call");
286
287 int nargin = args.length ();
288
289 if (nargin > 3 && nargin < 6)
290 {
291 std::string fcn_name, fname, jac_name, jname;
292 daspk_fcn = 0;
293 daspk_jac = 0;
294
295 octave_value f_arg = args(0);
296
297 if (f_arg.is_cell ())
298 {
299 Cell c = f_arg.cell_value ();
300 if (c.length () == 1)
301 f_arg = c(0);
302 else if (c.length () == 2)
303 {
304 if (c(0).is_function_handle () || c(0).is_inline_function ())
305 daspk_fcn = c(0).function_value ();
306 else
307 {
308 fcn_name = unique_symbol_name ("__daspk_fcn__");
309 fname = "function y = ";
310 fname.append (fcn_name);
311 fname.append (" (x, xdot, t) y = ");
312 daspk_fcn = extract_function
313 (c(0), "daspk", fcn_name, fname, "; endfunction");
314 }
315
316 if (daspk_fcn)
317 {
318 if (c(1).is_function_handle () || c(1).is_inline_function ())
319 daspk_jac = c(1).function_value ();
320 else
321 {
322 jac_name = unique_symbol_name ("__daspk_jac__");
323 jname = "function jac = ";
324 jname.append (jac_name);
325 jname.append (" (x, xdot, t, cj) jac = ");
326 daspk_jac = extract_function
327 (c(1), "daspk", jac_name, jname, "; endfunction");
328
329 if (!daspk_jac)
330 {
331 if (fcn_name.length ())
332 clear_function (fcn_name);
333 daspk_fcn = 0;
334 }
335 }
336 }
337 }
338 else
339 DASPK_ABORT1 ("incorrect number of elements in cell array");
340 }
341
342 if (!daspk_fcn && ! f_arg.is_cell ())
343 {
344 if (f_arg.is_function_handle () || f_arg.is_inline_function ())
345 daspk_fcn = f_arg.function_value ();
346 else
347 {
348 switch (f_arg.rows ())
349 {
350 case 1:
351 do
352 {
353 fcn_name = unique_symbol_name ("__daspk_fcn__");
354 fname = "function y = ";
355 fname.append (fcn_name);
356 fname.append (" (x, xdot, t) y = ");
357 daspk_fcn = extract_function
358 (f_arg, "daspk", fcn_name, fname, "; endfunction");
359 }
360 while (0);
361 break;
362
363 case 2:
364 {
365 string_vector tmp = f_arg.all_strings ();
366
367 if (! error_state)
368 {
369 fcn_name = unique_symbol_name ("__daspk_fcn__");
370 fname = "function y = ";
371 fname.append (fcn_name);
372 fname.append (" (x, xdot, t) y = ");
373 daspk_fcn = extract_function
374 (tmp(0), "daspk", fcn_name, fname, "; endfunction");
375
376 if (daspk_fcn)
377 {
378 jac_name = unique_symbol_name ("__daspk_jac__");
379 jname = "function jac = ";
380 jname.append (jac_name);
381 jname.append (" (x, xdot, t, cj) jac = ");
382 daspk_jac = extract_function
383 (tmp(1), "daspk", jac_name, jname,
384 "; endfunction");
385
386 if (!daspk_jac)
387 {
388 if (fcn_name.length ())
389 clear_function (fcn_name);
390 daspk_fcn = 0;
391 }
392 }
393 }
394 }
395 }
396 }
397 }
398
399 if (error_state || ! daspk_fcn)
400 DASPK_ABORT ();
401
402 ColumnVector state = ColumnVector (args(1).vector_value ());
403
404 if (error_state)
405 DASPK_ABORT1 ("expecting state vector as second argument");
406
407 ColumnVector deriv (args(2).vector_value ());
408
409 if (error_state)
410 DASPK_ABORT1 ("expecting derivative vector as third argument");
411
412 ColumnVector out_times (args(3).vector_value ());
413
414 if (error_state)
415 DASPK_ABORT1 ("expecting output time vector as fourth argument");
416
417 ColumnVector crit_times;
418 int crit_times_set = 0;
419 if (nargin > 4)
420 {
421 crit_times = ColumnVector (args(4).vector_value ());
422
423 if (error_state)
424 DASPK_ABORT1 ("expecting critical time vector as fifth argument");
425
426 crit_times_set = 1;
427 }
428
429 if (state.capacity () != deriv.capacity ())
430 DASPK_ABORT1 ("x and xdot must have the same size");
431
432 double tzero = out_times (0);
433
434 DAEFunc func (daspk_user_function);
435 if (daspk_jac)
436 func.set_jacobian_function (daspk_user_jacobian);
437
438 DASPK dae (state, deriv, tzero, func);
439 dae.set_options (daspk_opts);
440
441 Matrix output;
442 Matrix deriv_output;
443
444 if (crit_times_set)
445 output = dae.integrate (out_times, deriv_output, crit_times);
446 else
447 output = dae.integrate (out_times, deriv_output);
448
449 if (fcn_name.length ())
450 clear_function (fcn_name);
451 if (jac_name.length ())
452 clear_function (jac_name);
453
454 if (! error_state)
455 {
456 std::string msg = dae.error_message ();
457
458 retval(3) = msg;
459 retval(2) = static_cast<double> (dae.integration_state ());
460
461 if (dae.integration_ok ())
462 {
463 retval(1) = deriv_output;
464 retval(0) = output;
465 }
466 else
467 {
468 retval(1) = Matrix ();
469 retval(0) = Matrix ();
470
471 if (nargout < 3)
472 error ("daspk: %s", msg.c_str ());
473 }
474 }
475 }
476 else
477 print_usage ();
478
479 return retval;
480 }