annotate libinterp/corefcn/daspk.cc @ 22323:bac0d6f07a3e

maint: Update copyright notices for 2016.
author John W. Eaton <jwe@octave.org>
date Wed, 17 Aug 2016 01:05:19 -0400
parents 112b20240c87
children d0562b3159c7
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
1 /*
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
2
22323
bac0d6f07a3e maint: Update copyright notices for 2016.
John W. Eaton <jwe@octave.org>
parents: 21966
diff changeset
3 Copyright (C) 1996-2016 John W. Eaton
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
4
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
6
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
10 option) any later version.
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
11
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
15 for more details.
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
16
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
19 <http://www.gnu.org/licenses/>.
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
20
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
21 */
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
22
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
23 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21200
diff changeset
24 # include "config.h"
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
25 #endif
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
26
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
27 #include <string>
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
28
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
29 #include <iomanip>
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
30 #include <iostream>
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
31
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
32 #include "DASPK.h"
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
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
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
35 #include "error.h"
21100
e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents: 21041
diff changeset
36 #include "errwarn.h"
20940
48b2ad5ee801 maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents: 20918
diff changeset
37 #include "ovl.h"
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
38 #include "ov-fcn.h"
5729
e065f7c18bdc [project @ 2006-04-03 19:03:30 by jwe]
jwe
parents: 5642
diff changeset
39 #include "ov-cell.h"
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
40 #include "pager.h"
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
41 #include "unwind-prot.h"
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
42 #include "utils.h"
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
43 #include "variables.h"
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
44
3998
f6df65db67f9 [project @ 2002-07-24 18:10:39 by jwe]
jwe
parents: 3997
diff changeset
45 #include "DASPK-opts.cc"
f6df65db67f9 [project @ 2002-07-24 18:10:39 by jwe]
jwe
parents: 3997
diff changeset
46
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
47 // Global pointer for user defined function required by daspk.
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
48 static octave_function *daspk_fcn;
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
49
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
50 // Global pointer for optional user defined jacobian function.
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
51 static octave_function *daspk_jac;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
52
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
53 // Have we warned about imaginary values returned from user function?
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
54 static bool warned_fcn_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
55 static bool warned_jac_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
56
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
57 // Is this a recursive call?
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
58 static int call_depth = 0;
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
59
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
60 ColumnVector
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
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
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
63 {
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
64 ColumnVector retval;
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
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
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
67
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
68 octave_value_list args;
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4144
diff changeset
69
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
70 args(2) = t;
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4144
diff changeset
71 args(1) = xdot;
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4144
diff changeset
72 args(0) = x;
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
73
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
74 if (daspk_fcn)
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
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 }
20785
b6f2909e7f94 always throw exception after debugging with debug_on_error
John W. Eaton <jwe@octave.org>
parents: 20756
diff changeset
82 catch (octave_execution_exception& e)
20725
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
83 {
21100
e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents: 21041
diff changeset
84 err_user_supplied_eval (e, "daspk");
20725
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
85 }
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
86
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
87 int tlen = tmp.length ();
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
88 if (tlen == 0 || ! tmp(0).is_defined ())
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
89 err_user_supplied_eval ("daspk");
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
90
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
91 if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
92 {
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
93 warning ("daspk: ignoring imaginary part returned from user-supplied function");
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
94 warned_fcn_imaginary = true;
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
95 }
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
96
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
97 retval = tmp(0).vector_value ();
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
98
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
99 if (tlen > 1)
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
100 ires = tmp(1).idx_type_value ();
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
101
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
102 if (retval.is_empty ())
21100
e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents: 21041
diff changeset
103 err_user_supplied_eval ("daspk");
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
104 }
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
105
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
106 return retval;
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
107 }
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
108
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
109 Matrix
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
110 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
111 double t, double cj)
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
112 {
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
113 Matrix retval;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
114
20229
5dfaaaae784f Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents: 19697
diff changeset
115 assert (x.numel () == xdot.numel ());
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
116
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
117 octave_value_list args;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
118
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
119 args(3) = cj;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
120 args(2) = t;
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4144
diff changeset
121 args(1) = xdot;
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4144
diff changeset
122 args(0) = x;
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
123
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
124 if (daspk_jac)
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
125 {
20725
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
126 octave_value_list tmp;
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
127
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
128 try
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
129 {
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
130 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
131 }
20785
b6f2909e7f94 always throw exception after debugging with debug_on_error
John W. Eaton <jwe@octave.org>
parents: 20756
diff changeset
132 catch (octave_execution_exception& e)
20725
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
133 {
21100
e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents: 21041
diff changeset
134 err_user_supplied_eval (e, "daspk");
20725
7164a69472dd eliminate more uses of error state
John W. Eaton <jwe@octave.org>
parents: 20711
diff changeset
135 }
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
136
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
137 int tlen = tmp.length ();
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
138 if (tlen == 0 || ! tmp(0).is_defined ())
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
139 err_user_supplied_eval ("daspk");
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
140
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
141 if (! warned_jac_imaginary && tmp(0).is_complex_type ())
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
142 {
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
143 warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function");
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
144 warned_jac_imaginary = true;
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
145 }
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
146
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
147 retval = tmp(0).matrix_value ();
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
148
21118
3ac9f47fb04b Invert pattern if/code/else/err_XXX to if !/err_XXX/code.
Rik <rik@octave.org>
parents: 21100
diff changeset
149 if (retval.is_empty ())
21100
e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents: 21041
diff changeset
150 err_user_supplied_eval ("daspk");
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
151 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
152
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
153 return retval;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
154 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
155
15039
e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents: 14854
diff changeset
156 DEFUN (daspk, args, nargout,
21966
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
157 doc: /* -*- texinfo -*-
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
158 @deftypefn {} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} daspk (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
159 Solve the set of differential-algebraic equations
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
160 @tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
161 $$ 0 = f (x, \dot{x}, t) $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
162 with
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
163 $$ x(t_0) = x_0, \dot{x}(t_0) = \dot{x}_0 $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
164 @end tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
165 @ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
166
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
167 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
168 0 = f (x, xdot, t)
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
169 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
170
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
171 @noindent
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
172 with
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
173
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
174 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
175 x(t_0) = x_0, xdot(t_0) = xdot_0
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
176 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
177
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
178 @end ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
179 The solution is returned in the matrices @var{x} and @var{xdot},
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
180 with each row in the result matrices corresponding to one of the
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
181 elements in the vector @var{t}. The first element of @var{t}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
182 should be @math{t_0} and correspond to the initial state of the
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
183 system @var{x_0} and its derivative @var{xdot_0}, so that the first
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
184 row of the output @var{x} is @var{x_0} and the first row
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
185 of the output @var{xdot} is @var{xdot_0}.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
186
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
187 The first argument, @var{fcn}, is a string, inline, or function handle
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
188 that names the function @math{f} to call to compute the vector of
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
189 residuals for the set of equations. It must have the form
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
190
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
191 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
192 @var{res} = f (@var{x}, @var{xdot}, @var{t})
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
193 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
194
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
195 @noindent
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
196 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
197 scalar.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
198
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
199 If @var{fcn} is a two-element string array or a two-element cell array
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
200 of strings, inline functions, or function handles, the first element names
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
201 the function @math{f} described above, and the second element names a
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
202 function to compute the modified Jacobian
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
203 @tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
204 $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
205 J = {\partial f \over \partial x}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
206 + c {\partial f \over \partial \dot{x}}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
207 $$
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
208 @end tex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
209 @ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
210
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
211 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
212 @group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
213 df df
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
214 jac = -- + c ------
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
215 dx d xdot
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
216 @end group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
217 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
218
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
219 @end ifnottex
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
220
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
221 The modified Jacobian function must have the form
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
222
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
223 @example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
224 @group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
225
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
226 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
227
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
228 @end group
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
229 @end example
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
230
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
231 The second and third arguments to @code{daspk} specify the initial
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
232 condition of the states and their derivatives, and the fourth argument
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
233 specifies a vector of output times at which the solution is desired,
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
234 including the time corresponding to the initial condition.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
235
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
236 The set of initial states and derivatives are not strictly required to
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
237 be consistent. If they are not consistent, you must use the
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
238 @code{daspk_options} function to provide additional information so
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
239 that @code{daspk} can compute a consistent starting point.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
240
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
241 The fifth argument is optional, and may be used to specify a set of
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
242 times that the DAE solver should not integrate past. It is useful for
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
243 avoiding difficulties with singularities and points where there is a
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
244 discontinuity in the derivative.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
245
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
246 After a successful computation, the value of @var{istate} will be
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
247 greater than zero (consistent with the Fortran version of @sc{daspk}).
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
248
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
249 If the computation is not successful, the value of @var{istate} will be
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
250 less than zero and @var{msg} will contain additional information.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
251
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
252 You can use the function @code{daspk_options} to set optional
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
253 parameters for @code{daspk}.
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
254 @seealso{dassl}
112b20240c87 move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents: 21743
diff changeset
255 @end deftypefn */)
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
256 {
20918
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
257 int nargin = args.length ();
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
258
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
259 if (nargin < 4 || nargin > 5)
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
260 print_usage ();
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
261
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
262 warned_fcn_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
263 warned_jac_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
264
20918
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
265 octave_value_list retval (4);
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
266
21743
f4d7d0eb5b0c use namespace for unwind_protect class
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
267 octave::unwind_protect frame;
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
268
10066
2cd940306a06 make unwind_protect frames local
Jaroslav Hajek <highegg@gmail.com>
parents: 9377
diff changeset
269 frame.protect_var (call_depth);
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
270 call_depth++;
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
271
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
272 if (call_depth > 1)
20957
9db35d2042be maint: eliminate special cases of statements after error.
Rik <rik@octave.org>
parents: 20940
diff changeset
273 error ("daspk: invalid recursive call");
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
274
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
275 std::string fcn_name, fname, jac_name, jname;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
276 daspk_fcn = 0;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
277 daspk_jac = 0;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
278
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
279 octave_value f_arg = args(0);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
280
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
281 if (f_arg.is_cell ())
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
282 {
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
283 Cell c = f_arg.cell_value ();
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
284 if (c.numel () == 1)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
285 f_arg = c(0);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
286 else if (c.numel () == 2)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
287 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
288 if (c(0).is_function_handle () || c(0).is_inline_function ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
289 daspk_fcn = c(0).function_value ();
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
290 else
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
291 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
292 fcn_name = unique_symbol_name ("__daspk_fcn__");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
293 fname = "function y = ";
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
294 fname.append (fcn_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
295 fname.append (" (x, xdot, t) y = ");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
296 daspk_fcn = extract_function
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
297 (c(0), "daspk", fcn_name, fname, "; endfunction");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
298 }
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
299
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
300 if (daspk_fcn)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
301 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
302 if (c(1).is_function_handle () || c(1).is_inline_function ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
303 daspk_jac = c(1).function_value ();
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
304 else
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
305 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
306 jac_name = unique_symbol_name ("__daspk_jac__");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
307 jname = "function jac = ";
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
308 jname.append (jac_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
309 jname.append (" (x, xdot, t, cj) jac = ");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
310 daspk_jac = extract_function (c(1), "daspk", jac_name,
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
311 jname, "; endfunction");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
312
20918
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
313 if (! daspk_jac)
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
314 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
315 if (fcn_name.length ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
316 clear_function (fcn_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
317 daspk_fcn = 0;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
318 }
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
319 }
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
320 }
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
321 }
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
322 else
20957
9db35d2042be maint: eliminate special cases of statements after error.
Rik <rik@octave.org>
parents: 20940
diff changeset
323 error ("daspk: incorrect number of elements in cell array");
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
324 }
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
325
20918
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
326 if (! daspk_fcn && ! f_arg.is_cell ())
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
327 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
328 if (f_arg.is_function_handle () || f_arg.is_inline_function ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
329 daspk_fcn = f_arg.function_value ();
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
330 else
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
331 {
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
332 switch (f_arg.rows ())
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
333 {
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
334 case 1:
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
335 do
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
336 {
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
337 fcn_name = unique_symbol_name ("__daspk_fcn__");
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
338 fname = "function y = ";
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
339 fname.append (fcn_name);
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
340 fname.append (" (x, xdot, t) y = ");
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
341 daspk_fcn = extract_function (f_arg, "daspk", fcn_name,
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
342 fname, "; endfunction");
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
343 }
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
344 while (0);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
345 break;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
346
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
347 case 2:
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
348 {
20990
fc9cca99b2de Deprecate all_strings, replace with string_vector_value.
Rik <rik@octave.org>
parents: 20974
diff changeset
349 string_vector tmp = f_arg.string_vector_value ();
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
350
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
351 fcn_name = unique_symbol_name ("__daspk_fcn__");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
352 fname = "function y = ";
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
353 fname.append (fcn_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
354 fname.append (" (x, xdot, t) y = ");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
355 daspk_fcn = extract_function (tmp(0), "daspk", fcn_name,
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
356 fname, "; endfunction");
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
357
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
358 if (daspk_fcn)
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
359 {
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
360 jac_name = unique_symbol_name ("__daspk_jac__");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
361 jname = "function jac = ";
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
362 jname.append (jac_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
363 jname.append (" (x, xdot, t, cj) jac = ");
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
364 daspk_jac = extract_function (tmp(1), "daspk",
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
365 jac_name, jname,
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
366 "; endfunction");
5729
e065f7c18bdc [project @ 2006-04-03 19:03:30 by jwe]
jwe
parents: 5642
diff changeset
367
20918
6f0bd96f93c0 maint: Use new C++ archetype in more files.
Rik <rik@octave.org>
parents: 20892
diff changeset
368 if (! daspk_jac)
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
369 {
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
370 if (fcn_name.length ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
371 clear_function (fcn_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
372 daspk_fcn = 0;
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
373 }
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
374 }
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
375 }
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
376 }
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 10066
diff changeset
377 }
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
378 }
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
379
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
380 if (! daspk_fcn)
20957
9db35d2042be maint: eliminate special cases of statements after error.
Rik <rik@octave.org>
parents: 20940
diff changeset
381 return retval;
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
382
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
383 ColumnVector state = args(1).xvector_value ("daspk: initial state X_0 must be a vector");
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
384
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
385 ColumnVector deriv = args(2).xvector_value ("daspk: initial derivatives XDOT_0 must be a vector");
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
386
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
387 ColumnVector out_times = args(3).xvector_value ("daspk: output time variable T must be a vector");
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
388
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
389 ColumnVector crit_times;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
390 int crit_times_set = 0;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
391 if (nargin > 4)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
392 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
393 crit_times = args(4).xvector_value ("daspk: list of critical times T_CRIT must be a vector");
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
394
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
395 crit_times_set = 1;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
396 }
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
397
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
398 if (state.numel () != deriv.numel ())
20957
9db35d2042be maint: eliminate special cases of statements after error.
Rik <rik@octave.org>
parents: 20940
diff changeset
399 error ("daspk: X_0 and XDOT_0 must have the same size");
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
400
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
401 double tzero = out_times (0);
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
402
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
403 DAEFunc func (daspk_user_function);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
404 if (daspk_jac)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
405 func.set_jacobian_function (daspk_user_jacobian);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
406
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
407 DASPK dae (state, deriv, tzero, func);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
408 dae.set_options (daspk_opts);
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
409
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
410 Matrix output;
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
411 Matrix deriv_output;
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
412
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
413 if (crit_times_set)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
414 output = dae.integrate (out_times, deriv_output, crit_times);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
415 else
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
416 output = dae.integrate (out_times, deriv_output);
3997
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3964
diff changeset
417
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
418 if (fcn_name.length ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
419 clear_function (fcn_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
420 if (jac_name.length ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
421 clear_function (jac_name);
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
422
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
423 std::string msg = dae.error_message ();
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
424
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
425 if (dae.integration_ok ())
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
426 {
20892
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
427 retval(0) = output;
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
428 retval(1) = deriv_output;
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
429 }
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
430 else
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
431 {
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
432 if (nargout < 3)
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
433 error ("daspk: %s", msg.c_str ());
20892
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
434
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
435 retval(0) = Matrix ();
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
436 retval(1) = Matrix ();
20802
8bb38ba1bad6 eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents: 20785
diff changeset
437 }
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
438
20892
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
439 retval(2) = static_cast<double> (dae.integration_state ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
440 retval(3) = msg;
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
441
3912
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
442 return retval;
f56cd411adb4 [project @ 2002-04-28 03:12:27 by jwe]
jwe
parents:
diff changeset
443 }