Mercurial > hg > octave-avbm
changeset 4140:303b28a7a7e4
[project @ 2002-11-01 02:53:13 by jwe]
author | jwe |
---|---|
date | Fri, 01 Nov 2002 02:53:14 +0000 |
parents | 02ca908056e9 |
children | 8c710385c572 |
files | src/ChangeLog src/DLD-FUNCTIONS/daspk.cc src/DLD-FUNCTIONS/dasrt.cc src/DLD-FUNCTIONS/dassl.cc src/DLD-FUNCTIONS/fsolve.cc src/DLD-FUNCTIONS/lsode.cc src/DLD-FUNCTIONS/odessa.cc src/DLD-FUNCTIONS/quad.cc |
diffstat | 8 files changed, 350 insertions(+), 15 deletions(-) [+] |
line wrap: on
line diff
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,6 +1,22 @@ 2002-10-31 John W. Eaton <jwe@bevo.che.wisc.edu> - * DLD-FUNCTIONS/fsolve.cc (fsolve_user_function): Print warning if + * DLD-FUNCTIONS/fsolve.cc (fsolve_user_function, fsolve_user_jacobian): + Print warning if user returns complex value. + * DLD-FUNCTIONS/quad.cc (quad_user_function): Likewise. + * DLD-FUNCTIONS/lsode.cc (lsode_user_function, lsode_user_jacobian): + Likewise. + * DLD-FUNCTIONS/dassl.cc (dassl_user_function, dassl_user_jacobian): + Likewise. + * DLD-FUNCTIONS/dasrt.cc (dasrt_user_f, dasrt_user_cf, dasrt_user_j): + Likewise. + * DLD-FUNCTIONS/daspk.cc (daspk_user_function, daspk_user_jacobian): + Likewise. + + * DLD-FUNCTIONS/daspk.cc (daspk_user_jacobian): New function. + (Fdaspk): Handle extracting Jacobian from function argument. + + * DLD-FUNCTIONS/fsolve.cc (fsolve_user_function): New function. + (Ffsolve): Handle extracting Jacobian from function argument. * Makefile.in (%.oct): Depend on octave$(EXEEXT) so that octave will be built before any .oct files.
--- a/src/DLD-FUNCTIONS/daspk.cc +++ b/src/DLD-FUNCTIONS/daspk.cc @@ -46,6 +46,13 @@ // Global pointer for user defined function required by daspk. static octave_function *daspk_fcn; +// Global pointer for optional user defined jacobian function. +static octave_function *daspk_jac; + +// Have we warned about imaginary values returned from user function? +static bool warned_fcn_imaginary = false; +static bool warned_jac_imaginary = false; + // Is this a recursive call? static int call_depth = 0; @@ -99,6 +106,12 @@ int tlen = tmp.length (); if (tlen > 0 && tmp(0).is_defined ()) { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("daspk: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (tlen > 1) @@ -114,6 +127,76 @@ return retval; } +Matrix +daspk_user_jacobian (const ColumnVector& x, const ColumnVector& xdot, + double t, double cj) +{ + Matrix retval; + + int nstates = x.capacity (); + + assert (nstates == xdot.capacity ()); + + octave_value_list args; + + args(3) = cj; + args(2) = t; + + if (nstates > 1) + { + Matrix m1 (nstates, 1); + Matrix m2 (nstates, 1); + for (int i = 0; i < nstates; i++) + { + m1 (i, 0) = x (i); + m2 (i, 0) = xdot (i); + } + octave_value state (m1); + octave_value deriv (m2); + args(1) = deriv; + args(0) = state; + } + else + { + double d1 = x (0); + double d2 = xdot (0); + octave_value state (d1); + octave_value deriv (d2); + args(1) = deriv; + args(0) = state; + } + + if (daspk_jac) + { + octave_value_list tmp = daspk_jac->do_multi_index_op (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("daspk"); + return retval; + } + + int tlen = tmp.length (); + if (tlen > 0 && tmp(0).is_defined ()) + { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + + retval = tmp(0).matrix_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("daspk"); + } + else + gripe_user_supplied_eval ("daspk"); + } + + return retval; +} + #define DASPK_ABORT() \ do \ { \ @@ -235,6 +318,9 @@ { octave_value_list retval; + warned_fcn_imaginary = false; + warned_jac_imaginary = false; + unwind_protect::begin_frame ("Fdaspk"); unwind_protect_int (call_depth); @@ -247,12 +333,46 @@ if (nargin > 3 && nargin < 6) { - daspk_fcn = extract_function - (args(0), "daspk", "__daspk_fcn__", - "function res = __daspk_fcn__ (x, xdot, t) res = ", - "; endfunction"); + daspk_fcn = 0; + daspk_jac = 0; + + octave_value f_arg = args(0); + + switch (f_arg.rows ()) + { + case 1: + daspk_fcn = extract_function + (args(0), "daspk", "__daspk_fcn__", + "function res = __daspk_fcn__ (x, xdot, t) res = ", + "; endfunction"); + break; + + case 2: + { + string_vector tmp = f_arg.all_strings (); - if (! daspk_fcn) + if (! error_state) + { + daspk_fcn = extract_function + (tmp(0), "daspk", "__daspk_fcn__", + "function res = __daspk_fcn__ (x, xdot, t) res = ", + "; endfunction"); + + if (daspk_fcn) + { + daspk_jac = extract_function + (tmp(1), "daspk", "__daspk_jac__", + "function jac = __daspk_jac__ (x, xdot, t, cj) jac = ", + "; endfunction"); + + if (! daspk_jac) + daspk_fcn = 0; + } + } + } + } + + if (error_state || ! daspk_fcn) DASPK_ABORT (); ColumnVector state = ColumnVector (args(1).vector_value ()); @@ -288,6 +408,9 @@ double tzero = out_times (0); DAEFunc func (daspk_user_function); + if (daspk_jac) + func.set_jacobian_function (daspk_user_jacobian); + DASPK dae (state, deriv, tzero, func); dae.set_options (daspk_opts);
--- a/src/DLD-FUNCTIONS/dasrt.cc +++ b/src/DLD-FUNCTIONS/dasrt.cc @@ -48,6 +48,11 @@ static octave_function *dasrt_j; static octave_function *dasrt_cf; +// Have we warned about imaginary values returned from user function? +static bool warned_fcn_imaginary = false; +static bool warned_jac_imaginary = false; +static bool warned_cf_imaginary = false; + // Is this a recursive call? static int call_depth = 0; @@ -91,6 +96,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (error_state || retval.length () == 0) @@ -133,6 +144,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_cf_imaginary && tmp(0).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied constraint function"); + warned_cf_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (error_state || retval.length () == 0) @@ -197,6 +214,12 @@ int tlen = tmp.length (); if (tlen > 0 && tmp(0).is_defined ()) { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + retval = tmp(0).matrix_value (); if (error_state || retval.length () == 0) @@ -366,6 +389,10 @@ { octave_value_list retval; + warned_fcn_imaginary = false; + warned_jac_imaginary = false; + warned_cf_imaginary = false; + unwind_protect::begin_frame ("Fdasrt"); unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/dassl.cc +++ b/src/DLD-FUNCTIONS/dassl.cc @@ -49,6 +49,10 @@ // Global pointer for optional user defined jacobian function. static octave_function *dassl_jac; +// Have we warned about imaginary values returned from user function? +static bool warned_fcn_imaginary = false; +static bool warned_jac_imaginary = false; + // Is this a recursive call? static int call_depth = 0; @@ -102,6 +106,12 @@ int tlen = tmp.length (); if (tlen > 0 && tmp(0).is_defined ()) { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("dassl: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (tlen > 1) @@ -169,6 +179,12 @@ int tlen = tmp.length (); if (tlen > 0 && tmp(0).is_defined ()) { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("dassl: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + retval = tmp(0).matrix_value (); if (error_state || retval.length () == 0) @@ -302,6 +318,9 @@ { octave_value_list retval; + warned_fcn_imaginary = false; + warned_jac_imaginary = false; + unwind_protect::begin_frame ("Fdassl"); unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/fsolve.cc +++ b/src/DLD-FUNCTIONS/fsolve.cc @@ -46,8 +46,12 @@ // Global pointer for user defined function required by hybrd1. static octave_function *fsolve_fcn; +// Global pointer for optional user defined jacobian function. +static octave_function *fsolve_jac; + // Have we warned about imaginary values returned from user function? -static bool warned_imaginary = false; +static bool warned_fcn_imaginary = false; +static bool warned_jac_imaginary = false; // Is this a recursive call? static int call_depth = 0; @@ -91,7 +95,7 @@ { ColumnVector retval; - int n = x.capacity (); + int n = x.length (); octave_value_list args; args.resize (1); @@ -117,10 +121,10 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { - if (! warned_imaginary && tmp(0).is_complex_type ()) + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) { warning ("fsolve: ignoring imaginary part returned from user-supplied function"); - warned_imaginary = true; + warned_fcn_imaginary = true; } retval = ColumnVector (tmp(0).vector_value ()); @@ -135,6 +139,55 @@ return retval; } +Matrix +fsolve_user_jacobian (const ColumnVector& x) +{ + Matrix retval; + + int n = x.length (); + + octave_value_list args; + args.resize (1); + + if (n > 1) + { + Matrix m (n, 1); + for (int i = 0; i < n; i++) + m (i, 0) = x (i); + octave_value vars (m); + args(0) = vars; + } + else + { + double d = x (0); + octave_value vars (d); + args(0) = vars; + } + + if (fsolve_fcn) + { + octave_value_list tmp = fsolve_jac->do_multi_index_op (1, args); + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("fsolve: ignoring imaginary part returned from user-supplied jacobian function"); + warned_fcn_imaginary = true; + } + + retval = tmp(0).matrix_value (); + + if (error_state || retval.length () <= 0) + gripe_user_supplied_eval ("fsolve"); + } + else + gripe_user_supplied_eval ("fsolve"); + } + + return retval; +} + #define FSOLVE_ABORT() \ do \ { \ @@ -172,7 +225,8 @@ { octave_value_list retval; - warned_imaginary = false; + warned_fcn_imaginary = false; + warned_jac_imaginary = false; unwind_protect::begin_frame ("Ffsolve"); @@ -186,10 +240,46 @@ if (nargin == 2 && nargout < 4) { - fsolve_fcn = extract_function (args(0), "fsolve", "__fsolve_fcn__", - "function y = __fsolve_fcn__ (x) y = ", - "; endfunction"); - if (! fsolve_fcn) + fsolve_fcn = 0; + fsolve_jac = 0; + + octave_value f_arg = args(0); + + switch (f_arg.rows ()) + { + case 1: + fsolve_fcn = extract_function + (f_arg, "fsolve", "__fsolve_fcn__", + "function y = __fsolve_fcn__ (x) y = ", + "; endfunction"); + break; + + case 2: + { + string_vector tmp = f_arg.all_strings (); + + if (! error_state) + { + fsolve_fcn = extract_function + (tmp(0), "fsolve", "__fsolve_fcn__", + "function y = __fsolve_fcn__ (x) y = ", + "; endfunction"); + + if (fsolve_fcn) + { + fsolve_jac = extract_function + (tmp(1), "fsolve", "__fsolve_jac__", + "function jac = __fsolve_jac__ (x) jac = ", + "; endfunction"); + + if (! fsolve_jac) + fsolve_fcn = 0; + } + } + } + } + + if (error_state || ! fsolve_fcn) FSOLVE_ABORT (); ColumnVector x (args(1).vector_value ()); @@ -204,6 +294,9 @@ warning ("fsolve: can't compute path output yet"); NLFunc nleqn_fcn (fsolve_user_function); + if (fsolve_jac) + nleqn_fcn.set_jacobian_function (fsolve_user_jacobian); + NLEqn nleqn (x, nleqn_fcn); nleqn.set_options (fsolve_opts);
--- a/src/DLD-FUNCTIONS/lsode.cc +++ b/src/DLD-FUNCTIONS/lsode.cc @@ -51,6 +51,10 @@ // Global pointer for optional user defined jacobian function used by lsode. static octave_function *lsode_jac; +// Have we warned about imaginary values returned from user function? +static bool warned_fcn_imaginary = false; +static bool warned_jac_imaginary = false; + // Is this a recursive call? static int call_depth = 0; @@ -82,6 +86,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("lsode: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (error_state || retval.length () == 0) @@ -122,6 +132,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + retval = tmp(0).matrix_value (); if (error_state || retval.length () == 0) @@ -268,6 +284,9 @@ { octave_value_list retval; + warned_fcn_imaginary = false; + warned_jac_imaginary = false; + unwind_protect::begin_frame ("Flsode"); unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/odessa.cc +++ b/src/DLD-FUNCTIONS/odessa.cc @@ -51,6 +51,11 @@ static octave_function *odessa_j; static octave_function *odessa_b; +// Have we warned about imaginary values returned from user function? +static bool warned_fcn_imaginary = false; +static bool warned_jac_imaginary = false; +static bool warned_b_imaginary = false; + // Is this a recursive call? static int call_depth = 0; @@ -92,6 +97,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied function"); + warned_fcn_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (error_state || retval.length () == 0) @@ -143,6 +154,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_jac_imaginary && tmp(0).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied jacobian function"); + warned_jac_imaginary = true; + } + retval = tmp(0).matrix_value (); if (error_state || retval.length () == 0) @@ -202,6 +219,12 @@ if (tmp.length () > 0 && tmp(0).is_defined ()) { + if (! warned_b_imaginary && tmp(0).is_complex_type ()) + { + warning ("odessa: ignoring imaginary part returned from user-supplied inhomogeneity function"); + warned_b_imaginary = true; + } + retval = ColumnVector (tmp(0).vector_value ()); if (error_state || retval.length () == 0) @@ -376,6 +399,10 @@ { octave_value_list retval; + warned_fcn_imaginary = false; + warned_jac_imaginary = false; + warned_b_imaginary = false; + unwind_protect::begin_frame ("Fodessa"); unwind_protect_int (call_depth);
--- a/src/DLD-FUNCTIONS/quad.cc +++ b/src/DLD-FUNCTIONS/quad.cc @@ -51,6 +51,9 @@ // Global pointer for user defined function required by quadrature functions. static octave_function *quad_fcn; +// Have we warned about imaginary values returned from user function? +static bool warned_imaginary = false; + // Is this a recursive call? static int call_depth = 0; @@ -75,6 +78,12 @@ if (tmp.length () && tmp(0).is_defined ()) { + if (! warned_imaginary && tmp(0).is_complex_type ()) + { + warning ("quad: ignoring imaginary part returned from user-supplied function"); + warned_imaginary = true; + } + retval = tmp(0).double_value (); if (error_state) @@ -156,6 +165,8 @@ { octave_value_list retval; + warned_imaginary = false; + unwind_protect::begin_frame ("Fquad"); unwind_protect_int (call_depth);