Mercurial > hg > octave-nkf

comparison src/DLD-FUNCTIONS/daspk.cc @ 4140:303b28a7a7e4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-11-01 02:53:13 by jwe]
author jwe
date Fri, 01 Nov 2002 02:53:14 +0000
parents 19a1626b8d57
children b02ada83de67
comparison
equal deleted inserted replaced
4139:02ca908056e9 4140:303b28a7a7e4
44 #include "DASPK-opts.cc" 44 #include "DASPK-opts.cc"
45 45
46 // Global pointer for user defined function required by daspk. 46 // Global pointer for user defined function required by daspk.
47 static octave_function *daspk_fcn; 47 static octave_function *daspk_fcn;
48 48
49 // Global pointer for optional user defined jacobian function.
50 static octave_function *daspk_jac;
51
52 // Have we warned about imaginary values returned from user function?
53 static bool warned_fcn_imaginary = false;
54 static bool warned_jac_imaginary = false;
55
49 // Is this a recursive call? 56 // Is this a recursive call?
50 static int call_depth = 0; 57 static int call_depth = 0;
51 58
52 ColumnVector 59 ColumnVector
53 daspk_user_function (const ColumnVector& x, const ColumnVector& xdot, 60 daspk_user_function (const ColumnVector& x, const ColumnVector& xdot,
97 } 104 }
98 105
99 int tlen = tmp.length (); 106 int tlen = tmp.length ();
100 if (tlen > 0 && tmp(0).is_defined ()) 107 if (tlen > 0 && tmp(0).is_defined ())
101 { 108 {
109 if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
110 {
111 warning ("daspk: ignoring imaginary part returned from user-supplied function");
112 warned_fcn_imaginary = true;
113 }
114
102 retval = ColumnVector (tmp(0).vector_value ()); 115 retval = ColumnVector (tmp(0).vector_value ());
103 116
104 if (tlen > 1) 117 if (tlen > 1)
105 ires = tmp(1).int_value (); 118 ires = tmp(1).int_value ();
119
120 if (error_state || retval.length () == 0)
121 gripe_user_supplied_eval ("daspk");
122 }
123 else
124 gripe_user_supplied_eval ("daspk");
125 }
126
127 return retval;
128 }
129
130 Matrix
131 daspk_user_jacobian (const ColumnVector& x, const ColumnVector& xdot,
132 double t, double cj)
133 {
134 Matrix retval;
135
136 int nstates = x.capacity ();
137
138 assert (nstates == xdot.capacity ());
139
140 octave_value_list args;
141
142 args(3) = cj;
143 args(2) = t;
144
145 if (nstates > 1)
146 {
147 Matrix m1 (nstates, 1);
148 Matrix m2 (nstates, 1);
149 for (int i = 0; i < nstates; i++)
150 {
151 m1 (i, 0) = x (i);
152 m2 (i, 0) = xdot (i);
153 }
154 octave_value state (m1);
155 octave_value deriv (m2);
156 args(1) = deriv;
157 args(0) = state;
158 }
159 else
160 {
161 double d1 = x (0);
162 double d2 = xdot (0);
163 octave_value state (d1);
164 octave_value deriv (d2);
165 args(1) = deriv;
166 args(0) = state;
167 }
168
169 if (daspk_jac)
170 {
171 octave_value_list tmp = daspk_jac->do_multi_index_op (1, args);
172
173 if (error_state)
174 {
175 gripe_user_supplied_eval ("daspk");
176 return retval;
177 }
178
179 int tlen = tmp.length ();
180 if (tlen > 0 && tmp(0).is_defined ())
181 {
182 if (! warned_jac_imaginary && tmp(0).is_complex_type ())
183 {
184 warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function");
185 warned_jac_imaginary = true;
186 }
187
188 retval = tmp(0).matrix_value ();
106 189
107 if (error_state || retval.length () == 0) 190 if (error_state || retval.length () == 0)
108 gripe_user_supplied_eval ("daspk"); 191 gripe_user_supplied_eval ("daspk");
109 } 192 }
110 else 193 else
233 @end deftypefn\n\ 316 @end deftypefn\n\
234 @seealso{dassl}") 317 @seealso{dassl}")
235 { 318 {
236 octave_value_list retval; 319 octave_value_list retval;
237 320
321 warned_fcn_imaginary = false;
322 warned_jac_imaginary = false;
323
238 unwind_protect::begin_frame ("Fdaspk"); 324 unwind_protect::begin_frame ("Fdaspk");
239 325
240 unwind_protect_int (call_depth); 326 unwind_protect_int (call_depth);
241 call_depth++; 327 call_depth++;
242 328
245 331
246 int nargin = args.length (); 332 int nargin = args.length ();
247 333
248 if (nargin > 3 && nargin < 6) 334 if (nargin > 3 && nargin < 6)
249 { 335 {
250 daspk_fcn = extract_function 336 daspk_fcn = 0;
251 (args(0), "daspk", "__daspk_fcn__", 337 daspk_jac = 0;
252 "function res = __daspk_fcn__ (x, xdot, t) res = ", 338
253 "; endfunction"); 339 octave_value f_arg = args(0);
254 340
255 if (! daspk_fcn) 341 switch (f_arg.rows ())
342 {
343 case 1:
344 daspk_fcn = extract_function
345 (args(0), "daspk", "__daspk_fcn__",
346 "function res = __daspk_fcn__ (x, xdot, t) res = ",
347 "; endfunction");
348 break;
349
350 case 2:
351 {
352 string_vector tmp = f_arg.all_strings ();
353
354 if (! error_state)
355 {
356 daspk_fcn = extract_function
357 (tmp(0), "daspk", "__daspk_fcn__",
358 "function res = __daspk_fcn__ (x, xdot, t) res = ",
359 "; endfunction");
360
361 if (daspk_fcn)
362 {
363 daspk_jac = extract_function
364 (tmp(1), "daspk", "__daspk_jac__",
365 "function jac = __daspk_jac__ (x, xdot, t, cj) jac = ",
366 "; endfunction");
367
368 if (! daspk_jac)
369 daspk_fcn = 0;
370 }
371 }
372 }
373 }
374
375 if (error_state || ! daspk_fcn)
256 DASPK_ABORT (); 376 DASPK_ABORT ();
257 377
258 ColumnVector state = ColumnVector (args(1).vector_value ()); 378 ColumnVector state = ColumnVector (args(1).vector_value ());
259 379
260 if (error_state) 380 if (error_state)
286 DASPK_ABORT1 ("x and xdot must have the same size"); 406 DASPK_ABORT1 ("x and xdot must have the same size");
287 407
288 double tzero = out_times (0); 408 double tzero = out_times (0);
289 409
290 DAEFunc func (daspk_user_function); 410 DAEFunc func (daspk_user_function);
411 if (daspk_jac)
412 func.set_jacobian_function (daspk_user_jacobian);
413
291 DASPK dae (state, deriv, tzero, func); 414 DASPK dae (state, deriv, tzero, func);
292 dae.set_options (daspk_opts); 415 dae.set_options (daspk_opts);
293 416
294 Matrix output; 417 Matrix output;
295 Matrix deriv_output; 418 Matrix deriv_output;