Mercurial > hg > octave-max
changeset 2928:295f037b4b3e
[project @ 1997-05-05 05:32:33 by jwe]
line wrap: on
line diff
deleted file mode 100644 --- a/src/Array-oc.cc +++ /dev/null @@ -1,40 +0,0 @@ -/* - -Copyright (C) 1996, 1997 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Instantiate Arrays of octave_child objects. - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "Array.h" -#include "Array.cc" - -#include "sighandlers.h" - -template class Array<octave_child>; - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/src/Array-os.cc +++ /dev/null @@ -1,44 +0,0 @@ -/* - -Copyright (C) 1996, 1997 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Instantiate Arrays of octave_stream objects. - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "Array.h" -#include "Array.cc" - -#include "oct-stream.h" - -template class Array<scanf_format_elt*>; - -template class Array<printf_format_elt*>; - -template class Array<octave_stream*>; - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/src/Array-tc.cc +++ /dev/null @@ -1,52 +0,0 @@ -/* - -Copyright (C) 1996, 1997 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Instantiate Arrays of octave_values. - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "Array.h" -#include "Array.cc" - -#include "ov.h" - -extern template class Array<int>; -extern template class Array2<int>; -extern template class DiagArray2<int>; - -extern template class Array<double>; -extern template class Array2<double>; -extern template class DiagArray2<double>; - -extern template class Array<Complex>; -extern template class Array2<Complex>; -extern template class DiagArray2<Complex>; - -template class Array<octave_value>; - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/balance.cc @@ -0,0 +1,286 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Written by A. S. Hodel <scotte@eng.auburn.edu> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include "CmplxAEPBAL.h" +#include "CmplxAEPBAL.h" +#include "dbleAEPBAL.h" +#include "dbleAEPBAL.h" +#include "dbleGEPBAL.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (balance, args, nargout, + "AA = balance (A [, OPT]) or [[DD,] AA] = balance (A [, OPT])\n\ +\n\ +generalized eigenvalue problem:\n\ +\n\ + [cc, dd, aa, bb] = balance (a, b [, opt])\n\ +\n\ +where OPT is an optional single character argument as follows: \n\ +\n\ + N: no balancing; arguments copied, transformation(s) set to identity\n\ + P: permute argument(s) to isolate eigenvalues where possible\n\ + S: scale to improve accuracy of computed eigenvalues\n\ + B: (default) permute and scale, in that order. Rows/columns\n\ + of a (and b) that are isolated by permutation are not scaled\n\ +\n\ +[DD, AA] = balance (A, OPT) returns aa = dd*a*dd,\n\ +\n\ +[CC, DD, AA, BB] = balance (A, B, OPT) returns AA (BB) = CC*A*DD (CC*B*DD)") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 3 || nargout < 0 || nargout > 4) + { + print_usage ("balance"); + return retval; + } + + string bal_job; + int my_nargin; // # args w/o optional string arg + + // Determine if balancing option is listed. Set my_nargin to the + // number of matrix inputs. + + if (args(nargin-1).is_string ()) + { + bal_job = args(nargin-1).string_value (); + my_nargin = nargin-1; + } + else + { + bal_job = "B"; + my_nargin = nargin; + } + + octave_value arg_a = args(0); + + int a_nr = arg_a.rows (); + int a_nc = arg_a.columns (); + + // Check argument 1 dimensions. + + int arg_is_empty = empty_arg ("balance", a_nr, a_nc); + + if (arg_is_empty < 0) + return retval; + if (arg_is_empty > 0) + return octave_value_list (2, Matrix ()); + + if (a_nr != a_nc) + { + gripe_square_matrix_required ("balance"); + return retval; + } + + // Extract argument 1 parameter for both AEP and GEP. + + Matrix aa; + ComplexMatrix caa; + if (arg_a.is_complex_type ()) + caa = arg_a.complex_matrix_value (); + else + aa = arg_a.matrix_value (); + + if (error_state) + return retval; + + // Treat AEP/GEP cases. + + switch (my_nargin) + { + case 1: + + // Algebraic eigenvalue problem. + + if (arg_a.is_complex_type ()) + { + ComplexAEPBALANCE result (caa, bal_job); + + if (nargout == 0 || nargout == 1) + retval(0) = result.balanced_matrix (); + else + { + retval(1) = result.balanced_matrix (); + retval(0) = result.balancing_matrix (); + } + } + else + { + AEPBALANCE result (aa, bal_job); + + if (nargout == 0 || nargout == 1) + retval(0) = result.balanced_matrix (); + else + { + retval(1) = result.balanced_matrix (); + retval(0) = result.balancing_matrix (); + } + } + break; + + case 2: + { + // Generalized eigenvalue problem. + + // 1st we have to check argument 2 dimensions and type... + + octave_value arg_b = args(1); + + int b_nr = arg_b.rows (); + int b_nc = arg_b.columns (); + + // Check argument 2 dimensions -- must match arg 1. + + if (b_nr != b_nc || b_nr != a_nr) + { + gripe_nonconformant (); + return retval; + } + + // Now, extract the second matrix... + // Extract argument 1 parameter for both AEP and GEP. + + Matrix bb; + ComplexMatrix cbb; + if (arg_b.is_complex_type ()) + cbb = arg_b.complex_matrix_value (); + else + bb = arg_b.matrix_value (); + + if (error_state) + return retval; + + // Both matrices loaded, now let's check what kind of arithmetic: + + if (arg_a.is_complex_type () || arg_b.is_complex_type ()) + { + if (arg_a.is_real_type ()) + caa = aa; + + if (arg_b.is_real_type ()) + cbb = bb; + + // Compute magnitudes of elements for balancing purposes. + // Surely there's a function I can call someplace! + + for (int i = 0; i < a_nr; i++) + for (int j = 0; j < a_nc; j++) + { + aa (i, j) = abs (caa (i, j)); + bb (i, j) = abs (cbb (i, j)); + } + } + + GEPBALANCE result (aa, bb, bal_job); + + if (arg_a.is_complex_type () || arg_b.is_complex_type ()) + { + caa = result.left_balancing_matrix () * caa + * result.right_balancing_matrix (); + + cbb = result.left_balancing_matrix () * cbb + * result.right_balancing_matrix (); + + switch (nargout) + { + case 0: + case 1: + warning ("balance: should use two output arguments"); + retval(0) = caa; + break; + + case 2: + retval(1) = cbb; + retval(0) = caa; + break; + + case 4: + retval(3) = cbb; + retval(2) = caa; + retval(1) = result.right_balancing_matrix (); + retval(0) = result.left_balancing_matrix (); + break; + + default: + error ("balance: invalid number of output arguments"); + break; + } + } + else + { + switch (nargout) + { + case 0: + case 1: + warning ("balance: should use two output arguments"); + retval(0) = result.balanced_a_matrix (); + break; + + case 2: + retval(1) = result.balanced_b_matrix (); + retval(0) = result.balanced_a_matrix (); + break; + + case 4: + retval(3) = result.balanced_b_matrix (); + retval(2) = result.balanced_a_matrix (); + retval(1) = result.right_balancing_matrix (); + retval(0) = result.left_balancing_matrix (); + break; + + default: + error ("balance: invalid number of output arguments"); + break; + } + } + } + break; + + default: + error ("balance requires one (AEP) or two (GEP) numeric arguments"); + break; + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/chol.cc @@ -0,0 +1,103 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "CmplxCHOL.h" +#include "dbleCHOL.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (chol, args, nargout, + "R = chol (X): cholesky factorization") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1 || nargout > 1) + { + print_usage ("chol"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("chol", nr, nc); + + if (arg_is_empty < 0) + return retval; + if (arg_is_empty > 0) + return Matrix (); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int info; + CHOL fact (m, info); + if (info != 0) + error ("chol: matrix not positive definite"); + else + retval = fact.chol_matrix (); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int info; + ComplexCHOL fact (m, info); + if (info != 0) + error ("chol: matrix not positive definite"); + else + retval = fact.chol_matrix (); + } + } + else + { + gripe_wrong_type_arg ("chol", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/colloc.cc @@ -0,0 +1,141 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include "CollocWt.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (colloc, args, , + "[R, A, B, Q] = colloc (N [, \"left\"] [, \"right\"]): collocation weights") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 3) + { + print_usage ("colloc"); + return retval; + } + + if (! args(0).is_scalar_type ()) + { + error ("colloc: first argument must be a scalar"); + return retval; + } + + double tmp = args(0).double_value (); + + if (error_state) + return retval; + + if (xisnan (tmp)) + { + error ("colloc: NaN is invalid as NCOL"); + return retval; + } + + int ncol = NINT (tmp); + if (ncol < 0) + { + error ("colloc: first argument must be non-negative"); + return retval; + } + + int ntot = ncol; + int left = 0; + int right = 0; + + for (int i = 1; i < nargin; i++) + { + if (args(i).is_defined ()) + { + if (! args(i).is_string ()) + { + error ("colloc: expecting string argument"); + return retval; + } + + string s = args(i).string_value (); + + if ((s.length () == 1 && (s[0] == 'R' || s[0] == 'r')) + || s == "right") + { + right = 1; + } + else if ((s.length () == 1 && (s[0] == 'L' || s[0] == 'l')) + || s == "left") + { + left = 1; + } + else + { + error ("colloc: unrecognized argument"); + return retval; + } + } + else + { + error ("colloc: unexpected empty argument"); + return retval; + } + } + + ntot += left + right; + if (ntot < 1) + { + error ("colloc: the total number of roots must be positive"); + return retval; + } + + CollocWt wts (ncol, left, right); + + ColumnVector r = wts.roots (); + Matrix A = wts.first (); + Matrix B = wts.second (); + ColumnVector q = wts.quad_weights (); + + retval(3) = q; + retval(2) = B; + retval(1) = A; + retval(0) = r; + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/dassl.cc @@ -0,0 +1,377 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iostream.h> + +#include "DASSL.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "oct-sym.h" +#include "pager.h" +#include "utils.h" +#include "variables.h" + +// Global pointer for user defined function required by dassl. +static octave_symbol *dassl_fcn; + +static DASSL_options dassl_opts; + +ColumnVector +dassl_user_function (const ColumnVector& x, const ColumnVector& xdot, double t) +{ + ColumnVector retval; + + int nstates = x.capacity (); + + assert (nstates == xdot.capacity ()); + + octave_value_list args; + 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 (dassl_fcn) + { + octave_value_list tmp = dassl_fcn->eval (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("dassl"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).vector_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("dassl"); + } + else + gripe_user_supplied_eval ("dassl"); + } + + return retval; +} + +DEFUN_DLD (dassl, args, , + "dassl (\"function_name\", x_0, xdot_0, t_out)\n\ +dassl (F, X_0, XDOT_0, T_OUT, T_CRIT)\n\ +\n\ +The first argument is the name of the function to call to\n\ +compute the vector of residuals. It must have the form\n\ +\n\ + res = f (x, xdot, t)\n\ +\n\ +where x, xdot, and res are vectors, and t is a scalar.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 4 || nargin > 5) + { + print_usage ("dassl"); + return retval; + } + + dassl_fcn = extract_function + (args(0), "dassl", "__dassl_fcn__", + "function res = __dassl_fcn__ (x, xdot, t) res = ", + "; endfunction"); + + if (! dassl_fcn) + return retval; + + ColumnVector state = args(1).vector_value (); + + if (error_state) + { + error ("dassl: expecting state vector as second argument"); + return retval; + } + + ColumnVector deriv = args(2).vector_value (); + + if (error_state) + { + error ("dassl: expecting derivative vector as third argument"); + return retval; + } + + ColumnVector out_times = args(3).vector_value (); + + if (error_state) + { + error ("dassl: expecting output time vector as fourth argument"); + return retval; + } + + ColumnVector crit_times; + int crit_times_set = 0; + if (nargin > 4) + { + crit_times = args(4).vector_value (); + + if (error_state) + { + error ("dassl: expecting critical time vector as fifth argument"); + return retval; + } + + crit_times_set = 1; + } + + if (state.capacity () != deriv.capacity ()) + { + error ("dassl: x and xdot must have the same size"); + return retval; + } + + double tzero = out_times (0); + + DAEFunc func (dassl_user_function); + DASSL dae (state, deriv, tzero, func); + dae.copy (dassl_opts); + + Matrix output; + Matrix deriv_output; + + if (crit_times_set) + output = dae.integrate (out_times, deriv_output, crit_times); + else + output = dae.integrate (out_times, deriv_output); + + if (! error_state) + { + retval.resize (2); + + retval(0) = output; + retval(1) = deriv_output; + } + + return retval; +} + +typedef void (DASSL_options::*d_set_opt_mf) (double); +typedef double (DASSL_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 3 + +struct DASSL_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static DASSL_OPTIONS dassl_option_table [] = +{ + { "absolute tolerance", + { "absolute", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + DASSL_options::set_absolute_tolerance, + DASSL_options::absolute_tolerance, }, + + { "initial step size", + { "initial", "step", "size", 0, }, + { 1, 0, 0, 0, }, 1, + DASSL_options::set_initial_step_size, + DASSL_options::initial_step_size, }, + + { "maximum step size", + { "maximum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + DASSL_options::set_maximum_step_size, + DASSL_options::maximum_step_size, }, + + { "relative tolerance", + { "relative", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + DASSL_options::set_relative_tolerance, + DASSL_options::relative_tolerance, }, + + { 0, + { 0, 0, 0, 0, }, + { 0, 0, 0, 0, }, 0, + 0, 0, }, +}; + +static void +print_dassl_option_list (ostream& os) +{ + print_usage ("dassl_options", 1); + + os << "\n" + << "Options for dassl include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + DASSL_OPTIONS *list = dassl_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + + double val = (dassl_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_dassl_option (const string& keyword, double val) +{ + DASSL_OPTIONS *list = dassl_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (dassl_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("dassl_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_dassl_option (const string& keyword) +{ + octave_value retval; + + DASSL_OPTIONS *list = dassl_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + double val = (dassl_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + + return retval; + } + list++; + } + + warning ("dassl_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +DEFUN_DLD (dassl_options, args, , + "dassl_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for dassl. Keywords may be abbreviated\n\ +to the shortest match.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin == 0) + { + print_dassl_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_dassl_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_dassl_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("dassl_options"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/det.cc @@ -0,0 +1,129 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "CmplxDET.h" +#include "dbleDET.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (det, args, , + "det (X): determinant of a square matrix") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("det"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + if (nr == 0 && nc == 0) + { + retval = 1.0; + return retval; + } + + int arg_is_empty = empty_arg ("det", nr, nc); + if (arg_is_empty < 0) + return retval; + if (arg_is_empty > 0) + return Matrix (1, 1, 1.0); + + if (nr != nc) + { + gripe_square_matrix_required ("det"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int info; + double rcond = 0.0; + + DET det = m.determinant (info, rcond); + + double d = 0.0; + + if (info == -1) + warning ("det: matrix singular to machine precision, rcond = %g", + rcond); + else + d = det.value (); + + retval = d; + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int info; + double rcond = 0.0; + + ComplexDET det = m.determinant (info, rcond); + + Complex c = 0.0; + + if (info == -1) + warning ("det: matrix singular to machine precision, rcond = %g", + rcond); + else + c = det.value (); + + retval = c; + } + } + else + { + gripe_wrong_type_arg ("det", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/eig.cc @@ -0,0 +1,115 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "EIG.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (eig, args, nargout, + "eig (X) or [V, D] = eig (X): compute eigenvalues and eigenvectors of X") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1 || nargout > 2) + { + print_usage ("eig"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("eig", nr, nc); + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (2, Matrix ()); + + if (nr != nc) + { + gripe_square_matrix_required ("eig"); + return retval; + } + + Matrix tmp; + ComplexMatrix ctmp; + EIG result; + + if (arg.is_real_type ()) + { + tmp = arg.matrix_value (); + + if (error_state) + return retval; + else + result = EIG (tmp); + } + else if (arg.is_complex_type ()) + { + ctmp = arg.complex_matrix_value (); + + if (error_state) + return retval; + else + result = EIG (ctmp); + } + else + { + gripe_wrong_type_arg ("eig", tmp); + return retval; + } + + if (nargout == 0 || nargout == 1) + { + retval(0) = result.eigenvalues (), 1; + } + else + { + // Blame it on Matlab. + + ComplexDiagMatrix d (result.eigenvalues ()); + + retval(1) = d; + retval(0) = result.eigenvectors (); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/expm.cc @@ -0,0 +1,97 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Written by A. S. Hodel <scotte@eng.auburn.edu> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (expm, args, , + "expm (X): matrix exponential, e^A") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("expm"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("expm", nr, nc); + + if (arg_is_empty < 0) + return retval; + if (arg_is_empty > 0) + return Matrix (); + + if (nr != nc) + { + gripe_square_matrix_required ("expm"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (error_state) + return retval; + else + retval = m.expm (); + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (error_state) + return retval; + else + retval = m.expm (); + } + else + { + gripe_wrong_type_arg ("expm", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/fft.cc @@ -0,0 +1,120 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +// This function should be merged with Fifft. + +DEFUN_DLD (fft, args, , + "fft (X [, N]): fast fourier transform of a vector") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2) + { + print_usage ("fft"); + return retval; + } + + octave_value arg = args(0); + + int n_points = arg.rows (); + if (n_points == 1) + n_points = arg.columns (); + + if (nargin == 2) + { + double dval = args(1).double_value (); + if (xisnan (dval)) + error ("fft: NaN is invalid as the N_POINTS"); + else + n_points = NINT (dval); + } + + if (error_state) + return retval; + + if (n_points < 0) + { + error ("fft: number of points must be greater than zero"); + return retval; + } + + int arg_is_empty = empty_arg ("fft", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty || n_points == 0) + return Matrix (); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + m.resize (1, n_points, 0.0); + else + m.resize (n_points, m.columns (), 0.0); + retval = m.fourier (); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + m.resize (1, n_points, 0.0); + else + m.resize (n_points, m.columns (), 0.0); + retval = m.fourier (); + } + } + else + { + gripe_wrong_type_arg ("fft", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/fft2.cc @@ -0,0 +1,126 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +// This function should be merged with Fifft2. + +DEFUN_DLD (fft2, args, , + "fft2 (X [, N] [, M])\n\ +\n\ +two dimensional fast fourier transform of a vector") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 3) + { + print_usage ("fft2"); + return retval; + } + + octave_value arg = args(0); + + int n_rows = arg.rows (); + if (nargin > 1) + { + double dval = args(1).double_value (); + if (xisnan (dval)) + error ("fft2: NaN is invalid as N_ROWS"); + else + n_rows = NINT (dval); + } + + if (error_state) + return retval; + + int n_cols = arg.columns (); + if (nargin > 2) + { + double dval = args(2).double_value (); + if (xisnan (dval)) + error ("fft2: NaN is invalid as N_COLS"); + else + n_cols = NINT (dval); + } + + if (error_state) + return retval; + + if (n_rows < 0 || n_cols < 0) + { + error ("fft2: number of points must be greater than zero"); + return retval; + } + + int arg_is_empty = empty_arg ("fft2", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty || n_rows == 0 || n_cols == 0) + return Matrix (); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + m.resize (n_rows, n_cols, 0.0); + retval = m.fourier2d (); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + m.resize (n_rows, n_cols, 0.0); + retval = m.fourier2d (); + } + } + else + { + gripe_wrong_type_arg ("fft2", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/filter.cc @@ -0,0 +1,295 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Based on Tony Richardson's filter.m. +// +// Originally translated to C++ by KH (Kurt.Hornik@ci.tuwien.ac.at) +// with help from Fritz Leisch and Andreas Weingessel on Oct 20, 1994. +// +// Rewritten to use templates to handle both real and complex cases by +// jwe, Wed Nov 1 19:15:29 1995. + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" +#include "help.h" + +extern MArray<double> +filter (MArray<double>&, MArray<double>&, MArray<double>&); + +extern MArray<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&); + +template <class T> +MArray<T> +filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si) +{ + MArray<T> y; + + int a_len = a.length (); + int b_len = b.length (); + int x_len = x.length (); + + int si_len = si.length (); + + int ab_len = a_len > b_len ? a_len : b_len; + + b.resize (ab_len, 0.0); + + if (si.length () != ab_len - 1) + { + error ("filter: si must be a vector of length max (length (a), length (b)) - 1"); + return y; + } + + T norm = a (0); + + if (norm == 0.0) + { + error ("filter: the first element of a must be non-zero"); + return y; + } + + y.resize (x_len, 0.0); + + if (norm != 1.0) + b = b / norm; + + if (a_len > 1) + { + a.resize (ab_len, 0.0); + + if (norm != 1.0) + a = a / norm; + + for (int i = 0; i < x_len; i++) + { + y (i) = si (0) + b (0) * x (i); + + if (si_len > 1) + { + for (int j = 0; j < si_len - 1; j++) + si (j) = si (j+1) - a (j+1) * y (i) + + b (j+1) * x (i); + + si (si_len-1) = b (si_len) * x (i) + - a (si_len) * y (i); + } + else + si (0) = b (si_len) * x (i) + - a (si_len) * y (i); + } + } + else if (si_len > 0) + { + for (int i = 0; i < x_len; i++) + { + y (i) = si (0) + b (0) * x (i); + + if (si_len > 1) + { + for (int j = 0; j < si_len - 1; j++) + si (j) = si (j+1) + b (j+1) * x (i); + + si (si_len-1) = b (si_len) * x (i); + } + else + si (0) = b (1) * x (i); + } + } + else + y = b (0) * x; + + return y; +} + +extern MArray<double> +filter (MArray<double>&, MArray<double>&, MArray<double>&, + MArray<double>&); + +extern MArray<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, + MArray<Complex>&); + +template <class T> +MArray<T> +filter (MArray<T>& b, MArray<T>& a, MArray<T>& x) +{ + int a_len = a.length (); + int b_len = b.length (); + + int si_len = (a_len > b_len ? a_len : b_len) - 1; + + MArray<T> si (si_len, T (0.0)); + + return filter (b, a, x, si); +} + +DEFUN_DLD (filter, args, , + "usage: [y [, sf]] = filter (b, a, x [, si])\n\ +\n\ +y = filter (b, a, x) returns the solution to the following linear,\n\ +time-invariant difference equation:\n\ +\n\ + a[1] y[n] + ... + a[la] y[n-la+1] = b[1] x[n] + ... + b[lb] x[n-lb+1],\n\ +where la = length (a) and lb = length (b).\n\ +\n\ +[y, sf] = filter (b, a, x, si) sets the initial state of the system, si,\n\ +and returns the final state, sf. The state vector is a column vector\n\ +whose length is equal to the length of the longest coefficient vector\n\ +minus one. If si is not set, the initial state vector is set to all\n\ +zeros.\n\ +\n\ +The particular algorithm employed is known as a transposed Direct Form II\n\ +implementation.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 3 || nargin > 4) + { + print_usage ("filter"); + return retval; + } + + const char *errmsg = "filter: arguments must be vectors"; + + int x_is_vector = (args(2).rows () == 1 || args(2).columns () == 1); + + int si_is_vector = (nargin == 4 + && (args(3).rows () == 1 || args(3).columns () == 1)); + + if (args(0).is_complex_type () + || args(1).is_complex_type () + || args(2).is_complex_type () + || (nargin == 4 && args(3).is_complex_type ())) + { + ComplexColumnVector b = args(0).complex_vector_value (); + ComplexColumnVector a = args(1).complex_vector_value (); + ComplexColumnVector x = args(2).complex_vector_value (); + + if (! error_state) + { + if (nargin == 3) + { + ComplexColumnVector y (filter (b, a, x)); + + if (x_is_vector) + retval (0) = octave_value (y, (args(2).columns () == 1)); + else + retval (0) = y; + } + else + { + ComplexColumnVector si = args(3).complex_vector_value (); + + if (! error_state) + { + ComplexColumnVector y (filter (b, a, x, si)); + + if (si_is_vector) + retval (1) = octave_value (si, (args(3).columns () == 1)); + else + retval (1) = si; + + if (x_is_vector) + retval (0) = octave_value (y, (args(2).columns () == 1)); + else + retval (0) = y; + } + else + error (errmsg); + } + } + else + error (errmsg); + } + else + { + ColumnVector b = args(0).vector_value (); + ColumnVector a = args(1).vector_value (); + ColumnVector x = args(2).vector_value (); + + if (! error_state) + { + if (nargin == 3) + { + ColumnVector y (filter (b, a, x)); + + if (x_is_vector) + retval (0) = octave_value (y, (args(2).columns () == 1)); + else + retval (0) = y; + } + else + { + ColumnVector si = args(3).vector_value (); + + if (! error_state) + { + ColumnVector y (filter (b, a, x, si)); + + if (si_is_vector) + retval (1) = octave_value (si, (args(3).columns () == 1)); + else + retval (1) = si; + + if (x_is_vector) + retval (0) = octave_value (y, (args(2).columns () == 1)); + else + retval (0) = y; + } + else + error (errmsg); + } + } + else + error (errmsg); + } + + return retval; +} + +template MArray<double> +filter (MArray<double>&, MArray<double>&, MArray<double>&, + MArray<double>&); + +template MArray<double> +filter (MArray<double>&, MArray<double>&, MArray<double>&); + +template MArray<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, + MArray<Complex>&); + +template MArray<Complex> +filter (MArray<Complex>&, MArray<Complex>&, MArray <Complex>&); + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/find.cc @@ -0,0 +1,201 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" + +static octave_value_list +find_to_fortran_idx (const ColumnVector i_idx, const ColumnVector j_idx, + const octave_value& val, int nr, int nargout) +{ + octave_value_list retval; + + switch (nargout) + { + case 0: + case 1: + { + int count = i_idx.length (); + ColumnVector tmp (count); + for (int i = 0; i < count; i++) + tmp (i) = nr * (j_idx (i) - 1.0) + i_idx (i); + + // If the original argument was a row vector, force a row + // vector of indices to be returned. + + retval(0) = octave_value (tmp, (nr != 1)); + } + break; + + case 3: + retval(2) = val; + // Fall through! + + case 2: + retval(1) = octave_value (j_idx, 1); + retval(0) = octave_value (i_idx, 1); + + // If you want this to work more like Matlab, use + // + // retval(0) = octave_value (i_idx, (nr != 1)); + // + // instead of the previous statement. + + break; + + default: + panic_impossible (); + break; + } + + return retval; +} + +static octave_value_list +find_nonzero_elem_idx (const Matrix& m, int nargout) +{ + int count = 0; + int m_nr = m.rows (); + int m_nc = m.columns (); + + int i, j; + for (j = 0; j < m_nc; j++) + for (i = 0; i < m_nr; i++) + if (m (i, j) != 0.0) + count++; + + octave_value_list retval (((nargout == 0) ? 1 : nargout), Matrix ()); + + if (count == 0) + return retval; + + ColumnVector i_idx (count); + ColumnVector j_idx (count); + ColumnVector v (count); + + count = 0; + for (j = 0; j < m_nc; j++) + for (i = 0; i < m_nr; i++) + { + double d = m (i, j); + if (d != 0.0) + { + i_idx (count) = i + 1; + j_idx (count) = j + 1; + v (count) = d; + count++; + } + } + + octave_value tmp (v, 1); + return find_to_fortran_idx (i_idx, j_idx, tmp, m_nr, nargout); +} + +static octave_value_list +find_nonzero_elem_idx (const ComplexMatrix& m, int nargout) +{ + int count = 0; + int m_nr = m.rows (); + int m_nc = m.columns (); + + int i, j; + for (j = 0; j < m_nc; j++) + for (i = 0; i < m_nr; i++) + if (m (i, j) != 0.0) + count++; + + octave_value_list retval (((nargout == 0) ? 1 : nargout), Matrix ()); + + if (count == 0) + return retval; + + ColumnVector i_idx (count); + ColumnVector j_idx (count); + ComplexColumnVector v (count); + + count = 0; + for (j = 0; j < m_nc; j++) + for (i = 0; i < m_nr; i++) + { + Complex c = m (i, j); + if (c != 0.0) + { + i_idx (count) = i + 1; + j_idx (count) = j + 1; + v (count) = c; + count++; + } + } + + octave_value tmp (v, 1); + return find_to_fortran_idx (i_idx, j_idx, tmp, m_nr, nargout); +} + +DEFUN_DLD (find, args, nargout, + "find (X) or [I, J, V] = find (X): Return indices of nonzero elements") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1 || nargout > 3) + { + print_usage ("find"); + return retval; + } + + octave_value arg = args(0); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (m, nargout); + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (m, nargout); + } + else + { + gripe_wrong_type_arg ("find", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/fsolve.cc @@ -0,0 +1,336 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iostream.h> + +#include "NLEqn.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-sym.h" +#include "oct-obj.h" +#include "pager.h" +#include "utils.h" +#include "variables.h" + +// Global pointer for user defined function required by hybrd1. +static octave_symbol *fsolve_fcn; + +static NLEqn_options fsolve_opts; + +int +hybrd_info_to_fsolve_info (int info) +{ + switch (info) + { + case -1: + info = -2; + break; + + case 0: + info = -1; + break; + + case 1: + break; + + case 2: + info = 4; + break; + + case 3: + case 4: + case 5: + info = 3; + break; + + default: + panic_impossible (); + break; + } + + return info; +} + +ColumnVector +fsolve_user_function (const ColumnVector& x) +{ + ColumnVector retval; + + int n = x.capacity (); + + 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_fcn->eval (1, args); + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).vector_value (); + + if (error_state || retval.length () <= 0) + gripe_user_supplied_eval ("fsolve"); + } + else + gripe_user_supplied_eval ("fsolve"); + } + + return retval; +} + +DEFUN_DLD (fsolve, args, nargout, + "Solve nonlinear equations using Minpack. Usage:\n\ +\n\ + [X, INFO] = fsolve (F, X0)\n\ +\n\ +Where the first argument is the name of the function to call to\n\ +compute the vector of function values. It must have the form\n\ +\n\ + y = f (x)\n\ +\n\ +where y and x are vectors.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 2 || nargout > 3) + { + print_usage ("fsolve"); + return retval; + } + + fsolve_fcn = extract_function (args(0), "fsolve", "__fsolve_fcn__", + "function y = __fsolve_fcn__ (x) y = ", + "; endfunction"); + if (! fsolve_fcn) + return retval; + + ColumnVector x = args(1).vector_value (); + + if (error_state) + { + error ("fsolve: expecting vector as second argument"); + return retval; + } + + if (nargin > 2) + warning ("fsolve: ignoring extra arguments"); + + if (nargout > 2) + warning ("fsolve: can't compute path output yet"); + + NLFunc foo_fcn (fsolve_user_function); + NLEqn foo (x, foo_fcn); + foo.set_options (fsolve_opts); + + int info; + ColumnVector soln = foo.solve (info); + + info = hybrd_info_to_fsolve_info (info); + + retval.resize (nargout ? nargout : 1); + retval(0) = soln, 1; + + if (nargout > 1) + retval(1) = static_cast<double> (info); + + return retval; +} + +typedef void (NLEqn_options::*d_set_opt_mf) (double); +typedef double (NLEqn_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 1 + +struct NLEQN_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static NLEQN_OPTIONS fsolve_option_table [] = +{ + { "tolerance", + { "tolerance", 0, }, + { 1, 0, }, 1, + NLEqn_options::set_tolerance, + NLEqn_options::tolerance, }, + + { 0, + { 0, 0, }, + { 0, 0, }, 0, + 0, 0, }, +}; + +static void +print_fsolve_option_list (ostream& os) +{ + print_usage ("fsolve_options", 1); + + os << "\n" + << "Options for fsolve include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + NLEQN_OPTIONS *list = fsolve_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + + double val = (fsolve_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_fsolve_option (const string& keyword, double val) +{ + NLEQN_OPTIONS *list = fsolve_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (fsolve_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("fsolve_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_fsolve_option (const string& keyword) +{ + octave_value retval; + + NLEQN_OPTIONS *list = fsolve_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + double val = (fsolve_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + + return retval; + } + list++; + } + + warning ("fsolve_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +DEFUN_DLD (fsolve_options, args, , + "fsolve_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for fsolve. Keywords may be abbreviated\n\ +to the shortest match.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin == 0) + { + print_fsolve_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_fsolve_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_fsolve_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("fsolve_options"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/fsqp.cc @@ -0,0 +1,114 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "FSQP.h" + +#include "defun-dld.h" +#include "error.h" +#include "help.h" +#include "oct-obj.h" + +#ifndef FSQP_MISSING + +// Global pointers for user defined functions required by fsqp. +// static tree *fsqp_objective; +// static tree *fsqp_constraints; + +double +fsqp_objective_function (const ColumnVector&) +{ + return 0.0; +} + +ColumnVector +fsqp_constraint_function (const ColumnVector&) +{ + ColumnVector retval; + return retval; +} + +#endif + +#if defined (FSQP_MISSING) +DEFUN_DLD (fsqp, , , + "This function requires FSQP, which is not freely\n\ +redistributable. For more information, read the file\n\ +libcruft/fsqp/README.MISSING in the source distribution.") +#else +DEFUN_DLD (fsqp, , , + "[X, PHI] = fsqp (X, PHI [, LB, UB] [, LB, A, UB] [, LB, G, UB])\n\ +\n\ +Groups of arguments surrounded in `[]' are optional, but\n\ +must appear in the same relative order shown above.") +#endif +{ +/* + +Handle all of the following: + + 1. fsqp (x, phi) + 2. fsqp (x, phi, lb, ub) + 3. fsqp (x, phi, lb, ub, llb, c, lub) + 4. fsqp (x, phi, lb, ub, llb, c, lub, nllb, g, nlub) + 5. fsqp (x, phi, lb, ub, nllb, g, nlub) + 6. fsqp (x, phi, llb, c, lub, nllb, g, nlub) + 7. fsqp (x, phi, llb, c, lub) + 8. fsqp (x, phi, nllb, g, nlub) + +*/ + + octave_value_list retval; + + error ("fsqp: not implemented yet"); + + return retval; +} + +#if defined (FSQP_MISSING) +DEFUN_DLD (fsqp_options, , , + "This function requires FSQP, which is not freely\n\ +redistributable. For more information, read the file\n\ +libcruft/fsqp/README.MISSING in the source distribution.") +#else +DEFUN_DLD (fsqp_options, , , + "fsqp_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for fsqp. Keywords may be abbreviated\n\ +to the shortest match.") +#endif +{ + octave_value_list retval; + + error ("fsqp_options: not implemented yet"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/getgrent.cc @@ -0,0 +1,234 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#ifdef HAVE_SYS_TYPES_H +#include <sys/types.h> +#endif + +#ifdef HAVE_GRP_H +#include <grp.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-map.h" +#include "ov.h" +#include "oct-obj.h" +#include "utils.h" + +// Group file functions. (Why not?) + +static octave_value +mk_gr_map (struct group *gr) +{ + octave_value retval; + + if (gr) + { + Octave_map m; + + m ["name"] = gr->gr_name; +#if defined (HAVE_GR_PASSWD) + m ["passwd"] = gr->gr_passwd; +#else + m ["passwd"] = ""; +#endif + m ["gid"] = static_cast<double> (gr->gr_gid); + + if (gr->gr_mem) + { + // XXX FIXME XXX -- maybe there should be a string_vector + // constructor that takes a NULL terminated list of C + // strings. + + char **tmp = gr->gr_mem; + + int k = 0; + while (*tmp++) + k++; + + if (k > 0) + { + tmp = gr->gr_mem; + + string_vector members (k); + + for (int i = 0; i < k; i++) + members[i] = tmp[i]; + + m ["mem"] = members; + } + else + m ["mem"] = ""; + } + + retval = m; + } + else + retval = 0.0; + + return retval; +} + +DEFUN_DLD (getgrent, args, , + "getgrent ()\n\ +\n\ +Read an entry from the group-file stream, opening it if necessary.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 0) + { +#ifdef HAVE_GETGRENT + retval = mk_gr_map (getgrent ()); +#else + gripe_not_supported ("getgrent"); +#endif + } + else + print_usage ("getgrent"); + + return retval; +} + +DEFUN_DLD (getgrgid, args, , + "getgrgid (GID)\n\ +\n\ +Search for a group entry with a matching group ID.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { +#ifdef HAVE_GETGRGID + double dval = args(0).double_value (); + + if (! error_state) + { + if (D_NINT (dval) == dval) + { + gid_t gid = static_cast<gid_t> (dval); + + retval = mk_gr_map (getgrgid (gid)); + } + else + error ("getgrgid: argument must be an integer"); + } +#else + gripe_not_supported ("getgrgid"); +#endif + } + else + print_usage ("getgrgid"); + + return retval; +} + +DEFUN_DLD (getgrnam, args, , + "getgrnam (NAME)\n\ +\n\ +Search for group entry with a matching group name.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { +#ifdef HAVE_GETGRNAM + string s = args(0).string_value (); + + if (! error_state) + retval = mk_gr_map (getgrnam (s.c_str ())); +#else + gripe_not_supported ("getgrnam"); +#endif + } + else + print_usage ("getgrnam"); + + return retval; +} + +DEFUN_DLD (setgrent, args, , + "setgrent ()\n\ +\n\ +Rewind the group-file stream.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 0) + { +#ifdef HAVE_SETGRENT + setgrent (); +#else + gripe_not_supported ("setgrent"); +#endif + } + else + print_usage ("setgrent"); + + return retval; +} + +DEFUN_DLD (endgrent, args, , + "endgrent ()\n\ +\n\ +Close the group-file stream.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 0) + { +#ifdef HAVE_ENDGRENT + endgrent (); +#else + gripe_not_supported ("endgrent"); +#endif + } + else + print_usage ("endgrent"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/getpwent.cc @@ -0,0 +1,207 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#ifdef HAVE_SYS_TYPES_H +#include <sys/types.h> +#endif + +#ifdef HAVE_PWD_H +#include <pwd.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-map.h" +#include "ov.h" +#include "oct-obj.h" +#include "utils.h" + +// Password file functions. (Why not?) + +static octave_value +mk_pw_map (struct passwd *pw) +{ + octave_value retval; + + if (pw) + { + Octave_map m; + + m ["name"] = pw->pw_name; + m ["passwd"] = pw->pw_passwd; + m ["uid"] = static_cast<double> (pw->pw_uid); + m ["gid"] = static_cast<double> (pw->pw_gid); + m ["gecos"] = pw->pw_gecos; + m ["dir"] = pw->pw_dir; + m ["shell"] = pw->pw_shell; + + retval = m; + } + else + retval = 0.0; + + return retval; +} + +DEFUN_DLD (getpwent, args, , + "getpwent ()\n\ +\n\ +Read an entry from the password-file stream, opening it if necessary.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 0) + { +#ifdef HAVE_GETPWENT + retval = mk_pw_map (getpwent ()); +#else + gripe_not_supported ("getpwent"); +#endif + } + else + print_usage ("getpwent"); + + return retval; +} + +DEFUN_DLD (getpwuid, args, , + "getpwuid (UID)\n\ +\n\ +Search for a password entry with a matching user ID.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { +#ifdef HAVE_GETPWUID + double dval = args(0).double_value (); + + if (! error_state) + { + if (D_NINT (dval) == dval) + { + uid_t uid = static_cast<uid_t> (dval); + + retval = mk_pw_map (getpwuid (uid)); + } + else + error ("getpwuid: argument must be an integer"); + } +#else + gripe_not_supported ("getpwuid"); +#endif + } + else + print_usage ("getpwuid"); + + return retval; +} + +DEFUN_DLD (getpwnam, args, , + "getpwnam (NAME)\n\ +\n\ +Search for password entry with a matching username.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 1) + { +#ifdef HAVE_GETPWNAM + string s = args(0).string_value (); + + if (! error_state) + retval = mk_pw_map (getpwnam (s.c_str ())); +#else + gripe_not_supported ("getpwnam"); +#endif + } + else + print_usage ("getpwnam"); + + return retval; +} + +DEFUN_DLD (setpwent, args, , + "setpwent ()\n\ +\n\ +Rewind the password-file stream.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 0) + { +#ifdef HAVE_SETPWENT + setpwent (); +#else + gripe_not_supported ("setpwent"); +#endif + } + else + print_usage ("setpwent"); + + return retval; +} + +DEFUN_DLD (endpwent, args, , + "endpwent ()\n\ +\n\ +Close the password-file stream.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 0) + { +#ifdef HAVE_ENDPWENT + endpwent (); +#else + gripe_not_supported ("endpwent"); +#endif + } + else + print_usage ("endpwent"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/getrusage.cc @@ -0,0 +1,168 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "systime.h" + +#ifdef HAVE_SYS_RESOURCE_H +#include <sys/resource.h> +#endif + +#if defined (HAVE_TIMES) && defined (HAVE_SYS_TIMES_H) + +#if defined (HAVE_SYS_PARAM_H) +#include <sys/param.h> +#endif +#include <sys/times.h> + +#if !defined (HZ) +#if defined (CLK_TCK) +#define HZ CLK_TCK +#elif defined (USG) +#define HZ 100 +#else +#define HZ 60 +#endif +#endif + +#endif + +#include "defun-dld.h" +#include "help.h" +#include "oct-map.h" +#include "sysdep.h" +#include "ov.h" +#include "oct-obj.h" +#include "utils.h" + +#ifndef RUSAGE_SELF +#define RUSAGE_SELF 0 +#endif + +// System resource functions. + +DEFUN_DLD (getrusage, , , + "getrusage ()\n\ +\n\ +Return system resource statistics.") +{ + Octave_map m; + Octave_map tv_tmp; + +#if defined (HAVE_GETRUSAGE) + + struct rusage ru; + + getrusage (RUSAGE_SELF, &ru); + + tv_tmp ["sec"] = static_cast<double> (ru.ru_utime.tv_sec); + tv_tmp ["usec"] = static_cast<double> (ru.ru_utime.tv_usec); + m ["utime"] = octave_value (tv_tmp); + + tv_tmp ["sec"] = static_cast<double> (ru.ru_stime.tv_sec); + tv_tmp ["usec"] = static_cast<double> (ru.ru_stime.tv_usec); + m ["stime"] = octave_value (tv_tmp); + +#if ! defined (RUSAGE_TIMES_ONLY) + m ["maxrss"] = static_cast<double> (ru.ru_maxrss); + m ["ixrss"] = static_cast<double> (ru.ru_ixrss); + m ["idrss"] = static_cast<double> (ru.ru_idrss); + m ["isrss"] = static_cast<double> (ru.ru_isrss); + m ["minflt"] = static_cast<double> (ru.ru_minflt); + m ["majflt"] = static_cast<double> (ru.ru_majflt); + m ["nswap"] = static_cast<double> (ru.ru_nswap); + m ["inblock"] = static_cast<double> (ru.ru_inblock); + m ["oublock"] = static_cast<double> (ru.ru_oublock); + m ["msgsnd"] = static_cast<double> (ru.ru_msgsnd); + m ["msgrcv"] = static_cast<double> (ru.ru_msgrcv); + m ["nsignals"] = static_cast<double> (ru.ru_nsignals); + m ["nvcsw"] = static_cast<double> (ru.ru_nvcsw); + m ["nivcsw"] = static_cast<double> (ru.ru_nivcsw); +#endif + +#else +#if defined (HAVE_TIMES) && defined (HAVE_SYS_TIMES_H) + + struct tms t; + + times (&t); + + unsigned long ticks; + unsigned long seconds; + unsigned long fraction; + + ticks = t.tms_utime + t.tms_cutime; + fraction = ticks % HZ; + seconds = ticks / HZ; + + tv_tmp ["sec"] = static_cast<double> (seconds); + tv_tmp ["usec"] = static_cast<double> (fraction * 1e6 / HZ); + m ["utime"] = octave_value (tv_tmp); + + ticks = t.tms_stime + t.tms_cstime; + fraction = ticks % HZ; + seconds = ticks / HZ; + + tv_tmp ["sec"] = static_cast<double> (seconds); + tv_tmp ["usec"] = static_cast<double> (fraction * 1e6 / HZ); + m ["stime"] = octave_value (tv_tmp); + +#else + + tv_tmp ["sec"] = 0.0; + tv_tmp ["usec"] = 0.0; + m ["utime"] = octave_value (tv_tmp); + + tv_tmp ["sec"] = 0.0; + tv_tmp ["usec"] = 0.0; + m ["stime"] = octave_value (tv_tmp); + +#endif + + m ["maxrss"] = octave_NaN; + m ["ixrss"] = octave_NaN; + m ["idrss"] = octave_NaN; + m ["isrss"] = octave_NaN; + m ["minflt"] = octave_NaN; + m ["majflt"] = octave_NaN; + m ["nswap"] = octave_NaN; + m ["inblock"] = octave_NaN; + m ["oublock"] = octave_NaN; + m ["msgsnd"] = octave_NaN; + m ["msgrcv"] = octave_NaN; + m ["nsignals"] = octave_NaN; + m ["nvcsw"] = octave_NaN; + m ["nivcsw"] = octave_NaN; + +#endif + + return octave_value (m); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/givens.cc @@ -0,0 +1,122 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Originally written by A. S. Hodel <scotte@eng.auburn.edu> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "help.h" +#include "oct-obj.h" + +DEFUN_DLD (givens, args, nargout, + "G = givens (X, Y)\n\ +\n\ +compute orthogonal matrix G = [c s; -conj (s) c]\n\ +such that G [x; y] = [*; 0] (x, y scalars)\n\ +\n\ +[c, s] = givens (x, y) returns the (c, s) values themselves.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 2 || nargout > 2) + { + print_usage ("givens"); + return retval; + } + else + { + if (args(0).is_complex_type () || args(1).is_complex_type ()) + { + Complex cx = args(0).complex_value (); + Complex cy = args(1).complex_value (); + + if (! error_state) + { + ComplexMatrix result = Givens (cx, cy); + + if (! error_state) + { + switch (nargout) + { + case 0: + case 1: + retval(0) = result; + break; + + case 2: + retval(1) = result (0, 1); + retval(0) = result (0, 0); + break; + + default: + error ("givens: invalid number of output arguments"); + break; + } + } + } + } + else + { + double x = args(0).double_value (); + double y = args(1).double_value (); + + if (! error_state) + { + Matrix result = Givens (x, y); + + if (! error_state) + { + switch (nargout) + { + case 0: + case 1: + retval(0) = result; + break; + + case 2: + retval(1) = result (0, 1); + retval(0) = result (0, 0); + break; + + default: + error ("givens: invalid number of output arguments"); + break; + } + } + } + } + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/hess.cc @@ -0,0 +1,122 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "CmplxHESS.h" +#include "dbleHESS.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (hess, args, nargout, + "[P, H] = hess (A) or H = hess (A): Hessenberg decomposition") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1 || nargout > 2) + { + print_usage ("hess"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("hess", nr, nc); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (2, Matrix ()); + + if (nr != nc) + { + gripe_square_matrix_required ("hess"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix tmp = arg.matrix_value (); + + if (! error_state) + { + HESS result (tmp); + + if (nargout == 0 || nargout == 1) + { + retval.resize (1); + retval(0) = result.hess_matrix (); + } + else + { + retval.resize (2); + retval(0) = result.unitary_hess_matrix (); + retval(1) = result.hess_matrix (); + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix ctmp = arg.complex_matrix_value (); + + if (! error_state) + { + ComplexHESS result (ctmp); + + if (nargout == 0 || nargout == 1) + { + retval.resize (1); + retval(0) = result.hess_matrix (); + } + else + { + retval.resize (2); + retval(0) = result.unitary_hess_matrix (); + retval(1) = result.hess_matrix (); + } + } + } + else + { + gripe_wrong_type_arg ("hess", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/ifft.cc @@ -0,0 +1,120 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +// This function should be merged with Ffft. + +DEFUN_DLD (ifft, args, , + "ifft (X [, N]): inverse fast fourier transform of a vector") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2) + { + print_usage ("ifft"); + return retval; + } + + octave_value arg = args(0); + + int n_points = arg.rows (); + if (n_points == 1) + n_points = arg.columns (); + + if (nargin == 2) + { + double dval = args(1).double_value (); + if (xisnan (dval)) + error ("fft: NaN is invalid as the N_POINTS"); + else + n_points = NINT (dval); + } + + if (error_state) + return retval; + + if (n_points < 0) + { + error ("ifft: number of points must be greater than zero"); + return retval; + } + + int arg_is_empty = empty_arg ("ifft", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty || n_points == 0) + return Matrix (); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + m.resize (1, n_points, 0.0); + else + m.resize (n_points, m.columns (), 0.0); + retval = m.ifourier (); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + m.resize (1, n_points, 0.0); + else + m.resize (n_points, m.columns (), 0.0); + retval = m.ifourier (); + } + } + else + { + gripe_wrong_type_arg ("ifft", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/ifft2.cc @@ -0,0 +1,126 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +// This function should be merged with Ffft2. + +DEFUN_DLD (ifft2, args, , + "ifft2 (X [, N] [, M])\n\ +\n\ +two dimensional inverse fast fourier transform of a vector") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 3) + { + print_usage ("ifft2"); + return retval; + } + + octave_value arg = args(0); + + int n_rows = arg.rows (); + if (nargin > 1) + { + double dval = args(1).double_value (); + if (xisnan (dval)) + error ("fft2: NaN is invalid as N_ROWS"); + else + n_rows = NINT (dval); + } + + if (error_state) + return retval; + + int n_cols = arg.columns (); + if (nargin > 2) + { + double dval = args(2).double_value (); + if (xisnan (dval)) + error ("fft2: NaN is invalid as N_COLS"); + else + n_cols = NINT (dval); + } + + if (error_state) + return retval; + + if (n_rows < 0 || n_cols < 0) + { + error ("ifft2: number of points must be greater than zero"); + return retval; + } + + int arg_is_empty = empty_arg ("ifft2", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty || n_rows == 0 || n_cols == 0) + return Matrix (); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + m.resize (n_rows, n_cols, 0.0); + retval = m.ifourier2d (); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + m.resize (n_rows, n_cols, 0.0); + retval = m.ifourier2d (); + } + } + else + { + gripe_wrong_type_arg ("ifft2", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/inv.cc @@ -0,0 +1,119 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (inv, args, , + "inv (X): inverse of a square matrix") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("inv"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("inverse", nr, nc); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return Matrix (); + + if (nr != nc) + { + gripe_square_matrix_required ("inverse"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int info; + double rcond = 0.0; + + retval = m.inverse (info, rcond, 1); + + if (info == -1) + warning ("inverse: matrix singular to machine precision,\ + rcond = %g", rcond); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int info; + double rcond = 0.0; + + retval = m.inverse (info, rcond, 1); + + if (info == -1) + warning ("inverse: matrix singular to machine precision,\ + rcond = %g", rcond); + } + } + else + { + gripe_wrong_type_arg ("inv", arg); + } + + return retval; +} + +// XXX FIXME XXX -- this should really be done with an alias, but +// alias_builtin() won't do the right thing if we are actually using +// dynamic linking. + +DEFUN_DLD (inverse, args, nargout, + "inverse (X): inverse of a square matrix") +{ + return Finv (args, nargout); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/log.cc @@ -0,0 +1,267 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "EIG.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +// XXX FIXME XXX -- the next two functions should really be just +// one... + +DEFUN_DLD (logm, args, , + "logm (X): matrix logarithm") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("logm"); + return retval; + } + + octave_value arg = args(0); + + int arg_is_empty = empty_arg ("logm", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return Matrix (); + + if (arg.is_real_scalar ()) + { + double d = arg.double_value (); + if (d > 0.0) + retval(0) = log (d); + else + { + Complex dtmp (d); + retval(0) = log (dtmp); + } + } + else if (arg.is_complex_scalar ()) + { + Complex c = arg.complex_value (); + retval(0) = log (c); + } + else if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("logm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda (i) = log (real (elt)); + else + lambda (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("logm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda (i) = log (real (elt)); + else + lambda (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + } + else + { + gripe_wrong_type_arg ("logm", arg); + } + + return retval; +} + +DEFUN_DLD (sqrtm, args, , + "sqrtm (X): matrix sqrt") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("sqrtm"); + return retval; + } + + octave_value arg = args(0); + + int arg_is_empty = empty_arg ("sqrtm", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return Matrix (); + + if (arg.is_real_scalar ()) + { + double d = arg.double_value (); + if (d > 0.0) + retval(0) = sqrt (d); + else + { + Complex dtmp (d); + retval(0) = sqrt (dtmp); + } + } + else if (arg.is_complex_scalar ()) + { + Complex c = arg.complex_value (); + retval(0) = sqrt (c); + } + else if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("sqrtm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda (i) = sqrt (real (elt)); + else + lambda (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + int nr = m.rows (); + int nc = m.columns (); + + if (nr == 0 || nc == 0 || nr != nc) + gripe_square_matrix_required ("sqrtm"); + else + { + EIG m_eig (m); + ComplexColumnVector lambda (m_eig.eigenvalues ()); + ComplexMatrix Q (m_eig.eigenvectors ()); + + for (int i = 0; i < nr; i++) + { + Complex elt = lambda (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda (i) = sqrt (real (elt)); + else + lambda (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + } + else + { + gripe_wrong_type_arg ("sqrtm", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/lpsolve.cc @@ -0,0 +1,68 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "LPsolve.h" + +#include "defun-dld.h" +#include "error.h" +#include "help.h" +#include "oct-obj.h" + +DEFUN_DLD (lp_solve, , , + "lp_solve (): solve linear programs using lp_solve.") +{ + octave_value_list retval; + + // Force a bad value of inform, and empty matrices for x and phi. + + Matrix m; + retval(2) = -1.0; + retval(1) = m; + retval(0) = m; + + error ("lp_solve: not implemented yet"); + + return retval; +} + +DEFUN_DLD (lp_solve_options, , , + "lp_solve_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for lp_solve. Keywords may be abbreviated\n\ +to the shortest match.") +{ + octave_value_list retval; + + error ("lp_solve_options: not implemented yet"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/lsode.cc @@ -0,0 +1,472 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iostream.h> + +#include "LSODE.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "oct-sym.h" +#include "pager.h" +#include "utils.h" +#include "variables.h" + +// Global pointer for user defined function required by lsode. +static octave_symbol *lsode_fcn; + +// Global pointer for optional user defined jacobian function used by lsode. +static octave_symbol *lsode_jac; + +static LSODE_options lsode_opts; + +ColumnVector +lsode_user_function (const ColumnVector& x, double t) +{ + ColumnVector retval; + + int nstates = x.capacity (); + + octave_value_list args; + args(1) = t; + + Matrix m (nstates, 1); + for (int i = 0; i < nstates; i++) + m (i, 0) = x (i); + octave_value state (m); + args(0) = state; + + if (lsode_fcn) + { + octave_value_list tmp = lsode_fcn->eval (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("lsode"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).vector_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("lsode"); + } + else + gripe_user_supplied_eval ("lsode"); + } + + return retval; +} + +Matrix +lsode_user_jacobian (const ColumnVector& x, double t) +{ + Matrix retval; + + int nstates = x.capacity (); + + octave_value_list args; + args(1) = t; + + Matrix m (nstates, 1); + for (int i = 0; i < nstates; i++) + m (i, 0) = x (i); + octave_value state (m); + args(0) = state; + + if (lsode_jac) + { + octave_value_list tmp = lsode_jac->eval (1, args); + + if (error_state) + { + gripe_user_supplied_eval ("lsode"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).matrix_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("lsode"); + } + else + gripe_user_supplied_eval ("lsode"); + } + + return retval; +} + +DEFUN_DLD (lsode, args, nargout, + "lsode (F, X0, T_OUT, T_CRIT)\n\ +\n\ +The first argument is the name of the function to call to\n\ +compute the vector of right hand sides. It must have the form\n\ +\n\ + xdot = f (x, t)\n\ +\n\ +where xdot and x are vectors and t is a scalar.\n") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 3 || nargin > 4 || nargout > 1) + { + print_usage ("lsode"); + return retval; + } + + octave_value f_arg = args(0); + + switch (f_arg.rows ()) + { + case 1: + lsode_fcn = extract_function + (args(0), "lsode", "__lsode_fcn__", + "function xdot = __lsode_fcn__ (x, t) xdot = ", + "; endfunction"); + break; + + case 2: + { + string_vector tmp = args(0).all_strings (); + + if (! error_state) + { + lsode_fcn = extract_function + (tmp(0), "lsode", "__lsode_fcn__", + "function xdot = __lsode_fcn__ (x, t) xdot = ", + "; endfunction"); + + if (lsode_fcn) + { + lsode_jac = extract_function + (tmp(1), "lsode", "__lsode_jac__", + "function jac = __lsode_jac__ (x, t) jac = ", + "; endfunction"); + + if (! lsode_jac) + lsode_fcn = 0; + } + } + } + break; + + default: + error ("lsode: second arg should be a string or 2-element string array"); + break; + } + + if (error_state || ! lsode_fcn) + return retval; + + ColumnVector state = args(1).vector_value (); + + if (error_state) + { + error ("lsode: expecting state vector as second argument"); + return retval; + } + + ColumnVector out_times = args(2).vector_value (); + + if (error_state) + { + error ("lsode: expecting output time vector as third argument"); + return retval; + } + + ColumnVector crit_times; + + int crit_times_set = 0; + if (nargin > 3) + { + crit_times = args(3).vector_value (); + + if (error_state) + { + error ("lsode: expecting critical time vector as fourth argument"); + return retval; + } + + crit_times_set = 1; + } + + double tzero = out_times (0); + int nsteps = out_times.capacity (); + + ODEFunc func (lsode_user_function); + if (lsode_jac) + func.set_jacobian_function (lsode_user_jacobian); + + LSODE ode (state, tzero, func); + + ode.copy (lsode_opts); + + int nstates = state.capacity (); + Matrix output (nsteps, nstates + 1); + + if (crit_times_set) + output = ode.integrate (out_times, crit_times); + else + output = ode.integrate (out_times); + + if (! error_state) + { + retval.resize (1); + retval(0) = output; + } + + return retval; +} + +typedef void (LSODE_options::*d_set_opt_mf) (double); +typedef void (LSODE_options::*i_set_opt_mf) (int); +typedef double (LSODE_options::*d_get_opt_mf) (void); +typedef int (LSODE_options::*i_get_opt_mf) (void); + +#define MAX_TOKENS 3 + +struct LSODE_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + i_set_opt_mf i_set_fcn; + d_get_opt_mf d_get_fcn; + i_get_opt_mf i_get_fcn; +}; + +static LSODE_OPTIONS lsode_option_table [] = +{ + { "absolute tolerance", + { "absolute", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + LSODE_options::set_absolute_tolerance, 0, + LSODE_options::absolute_tolerance, 0, }, + + { "initial step size", + { "initial", "step", "size", 0, }, + { 1, 0, 0, 0, }, 1, + LSODE_options::set_initial_step_size, 0, + LSODE_options::initial_step_size, 0, }, + + { "maximum step size", + { "maximum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + LSODE_options::set_maximum_step_size, 0, + LSODE_options::maximum_step_size, 0, }, + + { "minimum step size", + { "minimum", "step", "size", 0, }, + { 2, 0, 0, 0, }, 1, + LSODE_options::set_minimum_step_size, 0, + LSODE_options::minimum_step_size, 0, }, + + { "relative tolerance", + { "relative", "tolerance", 0, 0, }, + { 1, 0, 0, 0, }, 1, + LSODE_options::set_relative_tolerance, 0, + LSODE_options::relative_tolerance, 0, }, + + { "step limit", + { "step", "limit", 0, 0, }, + { 1, 0, 0, 0, }, 1, + 0, LSODE_options::set_step_limit, + 0, LSODE_options::step_limit, }, + + { 0, + { 0, 0, 0, 0, }, + { 0, 0, 0, 0, }, 0, + 0, 0, 0, 0, }, +}; + +static void +print_lsode_option_list (ostream& os) +{ + print_usage ("lsode_options", 1); + + os << "\n" + << "Options for lsode include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + LSODE_OPTIONS *list = lsode_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + if (list->d_get_fcn) + { + double val = (lsode_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + } + else + { + int val = (lsode_opts.*list->i_get_fcn) (); + if (val < 0) + os << "infinite"; + else + os << val; + } + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_lsode_option (const string& keyword, double val) +{ + LSODE_OPTIONS *list = lsode_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_set_fcn) + (lsode_opts.*list->d_set_fcn) (val); + else + { + if (xisnan (val)) + { + error ("lsode_options: %s: expecting integer, found NaN", + keyword.c_str ()); + } + else + (lsode_opts.*list->i_set_fcn) (NINT (val)); + } + return; + } + list++; + } + + warning ("lsode_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_lsode_option (const string& keyword) +{ + octave_value retval; + + LSODE_OPTIONS *list = lsode_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_get_fcn) + { + double val = (lsode_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + } + else + { + int val = (lsode_opts.*list->i_get_fcn) (); + if (val < 0) + retval = "infinite"; + else + retval = static_cast<double> (val); + } + + return retval; + } + list++; + } + + warning ("lsode_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +DEFUN_DLD (lsode_options, args, , + "lsode_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for lsode. Keywords may be abbreviated\n\ +to the shortest match.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin == 0) + { + print_lsode_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_lsode_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_lsode_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("lsode_options"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/lu.cc @@ -0,0 +1,140 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "CmplxLU.h" +#include "dbleLU.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (lu, args, nargout, + "[L, U, P] = lu (A): LU factorization") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1 || nargout > 3) + { + print_usage ("lu"); + return retval; + } + + octave_value arg = args(0); + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("lu", nr, nc); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (3, Matrix ()); + + if (nr != nc) + { + gripe_square_matrix_required ("lu"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + LU fact (m); + + switch (nargout) + { + case 0: + case 1: + case 2: + { + Matrix P = fact.P (); + Matrix L = P.transpose () * fact.L (); + retval(1) = fact.U (); + retval(0) = L; + } + break; + + case 3: + default: + retval(2) = fact.P (); + retval(1) = fact.U (); + retval(0) = fact.L (); + break; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + ComplexLU fact (m); + + switch (nargout) + { + case 0: + case 1: + case 2: + { + ComplexMatrix P = fact.P (); + ComplexMatrix L = P.transpose () * fact.L (); + retval(1) = fact.U (); + retval(0) = L; + } + break; + + case 3: + default: + retval(2) = fact.P (); + retval(1) = fact.U (); + retval(0) = fact.L (); + break; + } + } + } + else + { + gripe_wrong_type_arg ("lu", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/minmax.cc @@ -0,0 +1,797 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "lo-ieee.h" +#include "oct-math.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" + +#ifndef MAX +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#endif + +#ifndef MIN +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#endif + +// XXX FIXME XXX -- it would be nice to share code among the min/max +// functions below. + +static Matrix +min (double d, const Matrix& m) +{ + int nr = m.rows (); + int nc = m.columns (); + + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double m_elem = m (i, j); + result (i, j) = MIN (d, m_elem); + } + + return result; +} + +static Matrix +min (const Matrix& m, double d) +{ + int nr = m.rows (); + int nc = m.columns (); + + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double m_elem = m (i, j); + result (i, j) = MIN (m_elem, d); + } + + return result; +} + +static ComplexMatrix +min (const Complex& c, const ComplexMatrix& m) +{ + int nr = m.rows (); + int nc = m.columns (); + + ComplexMatrix result (nr, nc); + + double abs_c = abs (c); + + for (int j = 0; j < nc; j++) + { + for (int i = 0; i < nr; i++) + { + double abs_m_elem = abs (m (i, j)); + if (abs_c < abs_m_elem) + result (i, j) = c; + else + result (i, j) = m (i, j); + } + } + + return result; +} + +static ComplexMatrix +min (const ComplexMatrix& m, const Complex& c) +{ + int nr = m.rows (); + int nc = m.columns (); + + ComplexMatrix result (nr, nc); + + double abs_c = abs (c); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double abs_m_elem = abs (m (i, j)); + if (abs_m_elem < abs_c) + result (i, j) = m (i, j); + else + result (i, j) = c; + } + + return result; +} + +static Matrix +min (const Matrix& a, const Matrix& b) +{ + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) + { + error ("two-arg min expecting args of same size"); + return Matrix (); + } + + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double a_elem = a (i, j); + double b_elem = b (i, j); + result (i, j) = MIN (a_elem, b_elem); + } + + return result; +} + +static ComplexMatrix +min (const ComplexMatrix& a, const ComplexMatrix& b) +{ + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) + { + error ("two-arg min expecting args of same size"); + return ComplexMatrix (); + } + + ComplexMatrix result (nr, nc); + + for (int j = 0; j < nc; j++) + { + int columns_are_real_only = 1; + for (int i = 0; i < nr; i++) + if (imag (a (i, j)) != 0.0 && imag (b (i, j)) != 0.0) + { + columns_are_real_only = 0; + break; + } + + if (columns_are_real_only) + { + for (int i = 0; i < nr; i++) + { + double a_elem = real (a (i, j)); + double b_elem = real (b (i, j)); + if (a_elem < b_elem) + result (i, j) = a_elem; + else + result (i, j) = b_elem; + } + } + else + { + for (int i = 0; i < nr; i++) + { + double abs_a_elem = abs (a (i, j)); + double abs_b_elem = abs (b (i, j)); + if (abs_a_elem < abs_b_elem) + result (i, j) = a (i, j); + else + result (i, j) = b (i, j); + } + } + } + + return result; +} + +static Matrix +max (double d, const Matrix& m) +{ + int nr = m.rows (); + int nc = m.columns (); + + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double m_elem = m (i, j); + result (i, j) = MAX (d, m_elem); + } + + return result; +} + +static Matrix +max (const Matrix& m, double d) +{ + int nr = m.rows (); + int nc = m.columns (); + + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double m_elem = m (i, j); + result (i, j) = MAX (m_elem, d); + } + + return result; +} + +static ComplexMatrix +max (const Complex& c, const ComplexMatrix& m) +{ + int nr = m.rows (); + int nc = m.columns (); + + ComplexMatrix result (nr, nc); + + double abs_c = abs (c); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double abs_m_elem = abs (m (i, j)); + if (abs_c > abs_m_elem) + result (i, j) = c; + else + result (i, j) = m (i, j); + } + + return result; +} + +static ComplexMatrix +max (const ComplexMatrix& m, const Complex& c) +{ + int nr = m.rows (); + int nc = m.columns (); + + ComplexMatrix result (nr, nc); + + double abs_c = abs (c); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double abs_m_elem = abs (m (i, j)); + if (abs_m_elem > abs_c) + result (i, j) = m (i, j); + else + result (i, j) = c; + } + + return result; +} + +static Matrix +max (const Matrix& a, const Matrix& b) +{ + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) + { + error ("two-arg max expecting args of same size"); + return Matrix (); + } + + Matrix result (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + { + double a_elem = a (i, j); + double b_elem = b (i, j); + result (i, j) = MAX (a_elem, b_elem); + } + + return result; +} + +static ComplexMatrix +max (const ComplexMatrix& a, const ComplexMatrix& b) +{ + int nr = a.rows (); + int nc = a.columns (); + if (nr != b.rows () || nc != b.columns ()) + { + error ("two-arg max expecting args of same size"); + return ComplexMatrix (); + } + + ComplexMatrix result (nr, nc); + + for (int j = 0; j < nc; j++) + { + int columns_are_real_only = 1; + for (int i = 0; i < nr; i++) + if (imag (a (i, j)) != 0.0 && imag (b (i, j)) != 0.0) + { + columns_are_real_only = 0; + break; + } + + if (columns_are_real_only) + { + for (int i = 0; i < nr; i++) + { + double a_elem = real (a (i, j)); + double b_elem = real (b (i, j)); + if (a_elem > b_elem) + result (i, j) = a_elem; + else + result (i, j) = b_elem; + } + } + else + { + for (int i = 0; i < nr; i++) + { + double abs_a_elem = abs (a (i, j)); + double abs_b_elem = abs (b (i, j)); + if (abs_a_elem > abs_b_elem) + result (i, j) = a (i, j); + else + result (i, j) = b (i, j); + } + } + } + + return result; +} + +DEFUN_DLD (min, args, nargout, + "min (X): minimum value(s) of a vector (matrix)") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2 || nargout > 2) + { + print_usage ("min"); + return retval; + } + + octave_value arg1; + octave_value arg2; + + switch (nargin) + { + case 2: + arg2 = args(1); + // Fall through... + + case 1: + arg1 = args(0); + break; + + default: + panic_impossible (); + break; + } + + if (nargin == 1 && (nargout == 1 || nargout == 0)) + { + if (arg1.is_real_type ()) + { + Matrix m = arg1.matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + retval(0) = m.row_min (); + else + retval(0) = octave_value (m.column_min (), 0); + } + } + else if (arg1.is_complex_type ()) + { + ComplexMatrix m = arg1.complex_matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + retval(0) = m.row_min (); + else + retval(0) = octave_value (m.column_min (), 0); + } + } + else + gripe_wrong_type_arg ("min", arg1); + } + else if (nargin == 1 && nargout == 2) + { + Array<int> index; + + if (arg1.is_real_type ()) + { + Matrix m = arg1.matrix_value (); + + if (! error_state) + { + retval.resize (2); + + if (m.rows () == 1) + retval(0) = m.row_min (index); + else + retval(0) = octave_value (m.column_min (index), 0); + } + } + else if (arg1.is_complex_type ()) + { + ComplexMatrix m = arg1.complex_matrix_value (); + + if (! error_state) + { + retval.resize (2); + + if (m.rows () == 1) + retval(0) = m.row_min (index); + else + retval(0) = octave_value (m.column_min (index), 0); + } + } + else + gripe_wrong_type_arg ("min", arg1); + + int len = index.length (); + + if (len > 0) + { + RowVector idx (len); + + for (int i = 0; i < len; i++) + { + int tmp = index.elem (i) + 1; + idx.elem (i) = (tmp <= 0) + ? octave_NaN : static_cast<double> (tmp); + } + + retval(1) = octave_value (idx, 0); + } + } + else if (nargin == 2) + { + int arg1_is_scalar = arg1.is_scalar_type (); + int arg2_is_scalar = arg2.is_scalar_type (); + + int arg1_is_complex = arg1.is_complex_type (); + int arg2_is_complex = arg2.is_complex_type (); + + if (arg1_is_scalar) + { + if (arg1_is_complex || arg2_is_complex) + { + Complex c1 = arg1.complex_value (); + ComplexMatrix m2 = arg2.complex_matrix_value (); + if (! error_state) + { + ComplexMatrix result = min (c1, m2); + if (! error_state) + retval(0) = result; + } + } + else + { + double d1 = arg1.double_value (); + Matrix m2 = arg2.matrix_value (); + + if (! error_state) + { + Matrix result = min (d1, m2); + if (! error_state) + retval(0) = result; + } + } + } + else if (arg2_is_scalar) + { + if (arg1_is_complex || arg2_is_complex) + { + ComplexMatrix m1 = arg1.complex_matrix_value (); + + if (! error_state) + { + Complex c2 = arg2.complex_value (); + ComplexMatrix result = min (m1, c2); + if (! error_state) + retval(0) = result; + } + } + else + { + Matrix m1 = arg1.matrix_value (); + + if (! error_state) + { + double d2 = arg2.double_value (); + Matrix result = min (m1, d2); + if (! error_state) + retval(0) = result; + } + } + } + else + { + if (arg1_is_complex || arg2_is_complex) + { + ComplexMatrix m1 = arg1.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix m2 = arg2.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix result = min (m1, m2); + if (! error_state) + retval(0) = result; + } + } + } + else + { + Matrix m1 = arg1.matrix_value (); + + if (! error_state) + { + Matrix m2 = arg2.matrix_value (); + + if (! error_state) + { + Matrix result = min (m1, m2); + if (! error_state) + retval(0) = result; + } + } + } + } + } + else + panic_impossible (); + + return retval; +} + +DEFUN_DLD (max, args, nargout, + "max (X): maximum value(s) of a vector (matrix)") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2 || nargout > 2) + { + print_usage ("max"); + return retval; + } + + octave_value arg1; + octave_value arg2; + + switch (nargin) + { + case 2: + arg2 = args(1); + // Fall through... + + case 1: + arg1 = args(0); + break; + + default: + panic_impossible (); + break; + } + + if (nargin == 1 && (nargout == 1 || nargout == 0)) + { + if (arg1.is_real_type ()) + { + Matrix m = arg1.matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + retval(0) = m.row_max (); + else + retval(0) = octave_value (m.column_max (), 0); + } + } + else if (arg1.is_complex_type ()) + { + ComplexMatrix m = arg1.complex_matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + retval(0) = m.row_max (); + else + retval(0) = octave_value (m.column_max (), 0); + } + } + else + gripe_wrong_type_arg ("max", arg1); + } + else if (nargin == 1 && nargout == 2) + { + Array<int> index; + + if (arg1.is_real_type ()) + { + Matrix m = arg1.matrix_value (); + + if (! error_state) + { + retval.resize (2); + + if (m.rows () == 1) + retval(0) = m.row_max (index); + else + retval(0) = octave_value (m.column_max (index), 0); + } + } + else if (arg1.is_complex_type ()) + { + ComplexMatrix m = arg1.complex_matrix_value (); + + if (! error_state) + { + retval.resize (2); + + if (m.rows () == 1) + retval(0) = m.row_max (index); + else + retval(0) = octave_value (m.column_max (index), 0); + } + } + else + gripe_wrong_type_arg ("max", arg1); + + int len = index.length (); + + if (len > 0) + { + RowVector idx (len); + + for (int i = 0; i < len; i++) + { + int tmp = index.elem (i) + 1; + idx.elem (i) = (tmp <= 0) + ? octave_NaN : static_cast<double> (tmp); + } + + retval(1) = octave_value (idx, 0); + } + } + else if (nargin == 2) + { + int arg1_is_scalar = arg1.is_scalar_type (); + int arg2_is_scalar = arg2.is_scalar_type (); + + int arg1_is_complex = arg1.is_complex_type (); + int arg2_is_complex = arg2.is_complex_type (); + + if (arg1_is_scalar) + { + if (arg1_is_complex || arg2_is_complex) + { + Complex c1 = arg1.complex_value (); + ComplexMatrix m2 = arg2.complex_matrix_value (); + if (! error_state) + { + ComplexMatrix result = max (c1, m2); + if (! error_state) + retval(0) = result; + } + } + else + { + double d1 = arg1.double_value (); + Matrix m2 = arg2.matrix_value (); + + if (! error_state) + { + Matrix result = max (d1, m2); + if (! error_state) + retval(0) = result; + } + } + } + else if (arg2_is_scalar) + { + if (arg1_is_complex || arg2_is_complex) + { + ComplexMatrix m1 = arg1.complex_matrix_value (); + + if (! error_state) + { + Complex c2 = arg2.complex_value (); + ComplexMatrix result = max (m1, c2); + if (! error_state) + retval(0) = result; + } + } + else + { + Matrix m1 = arg1.matrix_value (); + + if (! error_state) + { + double d2 = arg2.double_value (); + Matrix result = max (m1, d2); + if (! error_state) + retval(0) = result; + } + } + } + else + { + if (arg1_is_complex || arg2_is_complex) + { + ComplexMatrix m1 = arg1.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix m2 = arg2.complex_matrix_value (); + + if (! error_state) + { + ComplexMatrix result = max (m1, m2); + if (! error_state) + retval(0) = result; + } + } + } + else + { + Matrix m1 = arg1.matrix_value (); + + if (! error_state) + { + Matrix m2 = arg2.matrix_value (); + + if (! error_state) + { + Matrix result = max (m1, m2); + if (! error_state) + retval(0) = result; + } + } + } + } + } + else + panic_impossible (); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/npsol.cc @@ -0,0 +1,853 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iostream.h> + +#include "NPSOL.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "oct-sym.h" +#include "pager.h" +#include "utils.h" +#include "variables.h" + +#ifndef NPSOL_MISSING + +// Global pointers for user defined functions required by npsol. +static octave_symbol *npsol_objective; +static octave_symbol *npsol_constraints; + +static NPSOL_options npsol_opts; + +double +npsol_objective_function (const ColumnVector& x) +{ + int n = x.capacity (); + + octave_value decision_vars; + if (n > 1) + { + Matrix m (n, 1); + for (int i = 0; i < n; i++) + m (i, 0) = x (i); + decision_vars = m; + } + else + { + double d = x (0); + decision_vars = d; + } + + octave_value_list args; + args(0) = decision_vars; + + static double retval; + retval = 0.0; + + octave_value objective_value; + if (npsol_objective) + { + octave_value_list tmp = npsol_objective->eval (1, args); + + if (error_state) + { + error ("npsol: error evaluating objective function"); + npsol_objective_error = 1; // XXX FIXME XXX + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + objective_value = tmp(0); + else + { + error ("npsol: error evaluating objective function"); + npsol_objective_error = 1; // XXX FIXME XXX + return retval; + } + } + + if (objective_value.is_real_matrix ()) + { + Matrix m = objective_value.matrix_value (); + if (m.rows () == 1 && m.columns () == 1) + retval = m (0, 0); + else + { + gripe_user_returned_invalid ("npsol_objective"); + npsol_objective_error = 1; // XXX FIXME XXX + } + } + else if (objective_value.is_real_scalar ()) + { + retval = objective_value.double_value (); + } + else + { + gripe_user_returned_invalid ("npsol_objective"); + npsol_objective_error = 1; // XXX FIXME XXX + } + + return retval; +} + +ColumnVector +npsol_constraint_function (const ColumnVector& x) +{ + ColumnVector retval; + + int n = x.capacity (); + + octave_value decision_vars; + if (n > 1) + { + Matrix m (n, 1); + for (int i = 0; i < n; i++) + m (i, 0) = x (i); + decision_vars = m; + } + else + { + double d = x (0); + decision_vars = d; + } + + octave_value_list args; + args(0) = decision_vars; + + if (npsol_constraints) + { + octave_value_list tmp = npsol_constraints->eval (1, args); + + if (error_state) + { + error ("npsol: error evaluating constraints"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).vector_value (); + + if (error_state || retval.length () <= 0) + error ("npsol: error evaluating constraints"); + } + else + error ("npsol: error evaluating constraints"); + } + + return retval; +} + +static int +linear_constraints_ok (const ColumnVector& x, const ColumnVector& llb, + const Matrix& c, const ColumnVector& lub, + char *warn_for, int warn) +{ + int x_len = x.capacity (); + int llb_len = llb.capacity (); + int lub_len = lub.capacity (); + int c_rows = c.rows (); + int c_cols = c.columns (); + + int ok = 1; + if (warn) + { + if (c_rows == 0 || c_cols == 0 || llb_len == 0 || lub_len == 0) + { + ok = 0; + error ("%s: linear constraints must have nonzero dimensions", + warn_for); + } + else if (x_len != c_cols || llb_len != lub_len || llb_len != c_rows) + { + ok = 0; + error ("%s: linear constraints have inconsistent dimensions", + warn_for); + } + } + + return ok; +} + +static int +nonlinear_constraints_ok (const ColumnVector& x, const ColumnVector& nllb, + NLFunc::nonlinear_fcn g, const ColumnVector& nlub, + char *warn_for, int warn) +{ + int nllb_len = nllb.capacity (); + int nlub_len = nlub.capacity (); + ColumnVector c = (*g) (x); + int c_len = c.capacity (); + + int ok = 1; + if (warn) + { + if (nllb_len == 0 || nlub_len == 0 || c_len == 0) + { + ok = 0; + error ("%s: nonlinear constraints have nonzero dimensions", + warn_for); + } + else if (nllb_len != nlub_len || nllb_len != c_len) + { + ok = 0; + error ("%s: nonlinear constraints have inconsistent dimensions", + warn_for); + } + } + return ok; +} + +#endif + +#if defined (NPSOL_MISSING) +DEFUN_DLD (npsol, , , + "This function requires NPSOL, which is not freely\n\ +redistributable. For more information, read the file\n\ +libcruft/npsol/README.MISSING in the source distribution.") +#else +DEFUN_DLD (npsol, args, nargout, + "[X, OBJ, INFO, LAMBDA] = npsol (X, PHI [, LB, UB] [, A_LB, A, A_UB]\n\ + [, G_LB, G, G_UB])\n\ +\n\ +Groups of arguments surrounded in `[]' are optional, but\n\ +must appear in the same relative order shown above.\n\ +\n\ +The second argument is a string containing the name of the objective\n\ +function to call. The objective function must be of the form\n\ +\n\ + y = phi (x)\n\ +\n\ +where x is a vector and y is a scalar.\n\ +\n\ +The argument G is a string containing the name of the function that\n\ +defines the nonlinear constraints. It must be of the form\n\ +\n\ + y = g (x)\n\ +\n\ +where x is a vector and y is a vector.") +#endif +{ +/* + +Handle all of the following: + + 1. npsol (x, phi) + 2. npsol (x, phi, lb, ub) + 3. npsol (x, phi, lb, ub, llb, c, lub) + 4. npsol (x, phi, lb, ub, llb, c, lub, nllb, g, nlub) + 5. npsol (x, phi, lb, ub, nllb, g, nlub) + 6. npsol (x, phi, llb, c, lub, nllb, g, nlub) + 7. npsol (x, phi, llb, c, lub) + 8. npsol (x, phi, nllb, g, nlub) + +*/ + + octave_value_list retval; + +#if defined (NPSOL_MISSING) + + // Force a bad value of inform, and empty matrices for x, phi, and + // lambda. + + retval.resize (4, Matrix ()); + + retval(2) = -1.0; + + print_usage ("npsol"); + +#else + + int nargin = args.length (); + + if (nargin < 2 || nargin == 3 || nargin == 6 || nargin == 9 + || nargin > 10 || nargout > 4) + { + print_usage ("npsol"); + return retval; + } + + ColumnVector x = args(0).vector_value (); + + if (error_state || x.capacity () == 0) + { + error ("npsol: expecting vector as first argument"); + return retval; + } + + npsol_objective = extract_function + (args(1), "npsol", "__npsol_obj__", + "function phi = __npsol_obj__ (x) phi = ", + "; endfunction"); + + if (! npsol_objective) + return retval; + + Objective func (npsol_objective_function); + + ColumnVector soln; + + Bounds bounds; + if (nargin == 4 || nargin == 7 || nargin == 10) + { + ColumnVector lb = args(2).vector_value (); + ColumnVector ub = args(3).vector_value (); + + int lb_len = lb.capacity (); + int ub_len = ub.capacity (); + + if (error_state || lb_len != ub_len || lb_len != x.capacity ()) + { + error ("npsol: lower and upper bounds and decision variable vector"); + error ("must all have the same number of elements"); + return retval; + } + + bounds.resize (lb_len); + bounds.set_lower_bounds (lb); + bounds.set_upper_bounds (ub); + } + + double objf; + ColumnVector lambda; + int inform; + + if (nargin == 2) + { + // 1. npsol (x, phi) + + NPSOL nlp (x, func); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + + goto solved; + } + + if (nargin == 4) + { + // 2. npsol (x, phi, lb, ub) + + NPSOL nlp (x, func, bounds); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + + goto solved; + } + + npsol_constraints = 0; + if (nargin == 5 || nargin == 7 || nargin == 8 || nargin == 10) + npsol_constraints = extract_function + (args(nargin-2), "npsol", "__npsol_constr__", + "function y = __npsol_constr__ (x) y = ", + "; endfunction"); + + if (nargin == 7 || nargin == 5) + { + if (! npsol_constraints) + { + ColumnVector lub = args(nargin-1).vector_value (); + ColumnVector llb = args(nargin-3).vector_value (); + + if (error_state || llb.capacity () == 0 || lub.capacity () == 0) + { + error ("npsol: bounds for linear constraints must be vectors"); + return retval; + } + + Matrix c = args(nargin-2).matrix_value (); + + if (error_state) + { + error ("npsol: invalid linear constraint matrix"); + return retval; + } + + if (! linear_constraints_ok (x, llb, c, lub, "npsol", 1)) + return retval; + + LinConst linear_constraints (llb, c, lub); + + if (nargin == 5) + { + // 7. npsol (x, phi, llb, c, lub) + + NPSOL nlp (x, func, linear_constraints); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + } + else + { + // 3. npsol (x, phi, lb, ub, llb, c, lub) + + NPSOL nlp (x, func, bounds, linear_constraints); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + } + goto solved; + } + else + { + ColumnVector nlub = args(nargin-1).vector_value (); + ColumnVector nllb = args(nargin-3).vector_value (); + + if (error_state + || (! nonlinear_constraints_ok + (x, nllb, npsol_constraint_function, nlub, "npsol", 1))) + return retval; + + NLFunc const_func (npsol_constraint_function); + NLConst nonlinear_constraints (nllb, const_func, nlub); + + if (nargin == 5) + { + // 8. npsol (x, phi, nllb, g, nlub) + + NPSOL nlp (x, func, nonlinear_constraints); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + } + else + { + // 5. npsol (x, phi, lb, ub, nllb, g, nlub) + + NPSOL nlp (x, func, bounds, nonlinear_constraints); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + } + goto solved; + } + } + + if (nargin == 8 || nargin == 10) + { + if (! npsol_constraints) + { + // Produce error message. + + is_valid_function (args(nargin-2), "npsol", 1); + } + else + { + ColumnVector nlub = args(nargin-1).vector_value (); + ColumnVector nllb = args(nargin-3).vector_value (); + + if (error_state + || (! nonlinear_constraints_ok + (x, nllb, npsol_constraint_function, nlub, "npsol", 1))) + return retval; + + NLFunc const_func (npsol_constraint_function); + NLConst nonlinear_constraints (nllb, const_func, nlub); + + ColumnVector lub = args(nargin-4).vector_value (); + ColumnVector llb = args(nargin-6).vector_value (); + + if (error_state || llb.capacity () == 0 || lub.capacity () == 0) + { + error ("npsol: bounds for linear constraints must be vectors"); + return retval; + } + + Matrix c = args(nargin-5).matrix_value (); + + if (error_state) + { + error ("npsol: invalid linear constraint matrix"); + return retval; + } + + if (! linear_constraints_ok (x, llb, c, lub, "npsol", 1)) + return retval; + + LinConst linear_constraints (llb, c, lub); + + if (nargin == 8) + { + // 6. npsol (x, phi, llb, c, lub, nllb, g, nlub) + + NPSOL nlp (x, func, linear_constraints, + nonlinear_constraints); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + } + else + { + // 4. npsol (x, phi, lb, ub, llb, c, lub, nllb, g, nlub) + + NPSOL nlp (x, func, bounds, linear_constraints, + nonlinear_constraints); + nlp.set_options (npsol_opts); + soln = nlp.minimize (objf, inform, lambda); + } + goto solved; + } + } + + return retval; + + solved: + + retval.resize (nargout ? nargout : 1); + retval(0) = soln, 1; + if (nargout > 1) + retval(1) = objf; + if (nargout > 2) + retval(2) = static_cast<double> (inform); + if (nargout > 3) + retval(3) = lambda; + +#endif + + return retval; +} + +#ifndef NPSOL_MISSING + +typedef void (NPSOL_options::*d_set_opt_mf) (double); +typedef void (NPSOL_options::*i_set_opt_mf) (int); +typedef double (NPSOL_options::*d_get_opt_mf) (void); +typedef int (NPSOL_options::*i_get_opt_mf) (void); + +#define MAX_TOKENS 5 + +struct NPSOL_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + i_set_opt_mf i_set_fcn; + d_get_opt_mf d_get_fcn; + i_get_opt_mf i_get_fcn; +}; + +static NPSOL_OPTIONS npsol_option_table [] = +{ + { "central difference interval", + { "central", "difference", "interval", 0, 0, 0, }, + { 2, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_central_difference_interval, 0, + NPSOL_options::central_difference_interval, 0, }, + + { "crash tolerance", + { "crash", "tolerance", 0, 0, 0, 0, }, + { 2, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_crash_tolerance, 0, + NPSOL_options::crash_tolerance, 0, }, + + { "derivative level", + { "derivative", "level", 0, 0, 0, 0, }, + { 1, 0, 0, 0, 0, 0, }, 1, + 0, NPSOL_options::set_derivative_level, + 0, NPSOL_options::derivative_level, }, + + { "difference interval", + { "difference", "interval", 0, 0, 0, 0, }, + { 3, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_difference_interval, 0, + NPSOL_options::difference_interval, 0, }, + + { "function precision", + { "function", "precision", 0, 0, 0, 0, }, + { 2, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_function_precision, 0, + NPSOL_options::function_precision, 0, }, + + { "infinite bound size", + { "infinite", "bound", "size", 0, 0, 0, }, + { 1, 1, 0, 0, 0, 0, }, 2, + NPSOL_options::set_infinite_bound, 0, + NPSOL_options::infinite_bound, 0, }, + + { "infinite step size", + { "infinite", "step", "size", 0, 0, 0, }, + { 1, 1, 0, 0, 0, 0, }, 2, + NPSOL_options::set_infinite_step, 0, + NPSOL_options::infinite_step, 0, }, + + { "linear feasibility tolerance", + { "linear", "feasibility", "tolerance", 0, 0, 0, }, + { 5, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_linear_feasibility_tolerance, 0, + NPSOL_options::linear_feasibility_tolerance, 0, }, + + { "linesearch tolerance", + { "linesearch", "tolerance", 0, 0, 0, 0, }, + { 5, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_linesearch_tolerance, 0, + NPSOL_options::linesearch_tolerance, 0, }, + + { "major iteration limit", + { "major", "iteration", "limit", 0, 0, 0, }, + { 2, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_major_iteration_limit, + 0, NPSOL_options::major_iteration_limit, }, + + { "minor iteration limit", + { "minor", "iteration", "limit", 0, 0, 0, }, + { 2, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_minor_iteration_limit, + 0, NPSOL_options::minor_iteration_limit, }, + + { "major print level", + { "major", "print", "level", 0, 0, 0, }, + { 2, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_major_print_level, + 0, NPSOL_options::major_print_level, }, + + { "minor print level", + { "minor", "print", "level", 0, 0, 0, }, + { 2, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_minor_print_level, + 0, NPSOL_options::minor_print_level, }, + + { "nonlinear feasibility tolerance", + { "nonlinear", "feasibility", "tolerance", 0, 0, }, + { 1, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_nonlinear_feasibility_tolerance, 0, + NPSOL_options::nonlinear_feasibility_tolerance, 0, }, + + { "optimality tolerance", + { "optimality", "tolerance", 0, 0, 0, 0, }, + { 1, 0, 0, 0, 0, 0, }, 1, + NPSOL_options::set_optimality_tolerance, 0, + NPSOL_options::optimality_tolerance, 0, }, + + { "start objective check at variable", + { "start", "objective", "check", "at", "variable", 0, }, + { 3, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_start_objective_check, + 0, NPSOL_options::start_objective_check, }, + + { "start constraint check at variable", + { "start", "constraint", "check", "at", "variable", 0, }, + { 3, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_start_constraint_check, + 0, NPSOL_options::start_constraint_check, }, + + { "stop objective check at variable", + { "stop", "objective", "check", "at", "variable", 0, }, + { 3, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_stop_objective_check, + 0, NPSOL_options::stop_objective_check, }, + + { "stop constraint check at variable", + { "stop", "constraint", "check", "at", "variable", 0, }, + { 3, 1, 0, 0, 0, 0, }, 2, + 0, NPSOL_options::set_stop_constraint_check, + 0, NPSOL_options::stop_constraint_check, }, + + { "verify level", + { "verify", "level", 0, 0, 0, 0, }, + { 1, 0, 0, 0, 0, 0, }, 1, + 0, NPSOL_options::set_verify_level, + 0, NPSOL_options::verify_level, }, + + { 0, + { 0, 0, 0, 0, 0, 0, }, + { 0, 0, 0, 0, 0, 0, }, 0, + 0, 0, 0, 0, }, +}; + +static void +print_npsol_option_list (ostream& os) +{ + print_usage ("npsol_options", 1); + + os << "\n" + << "Options for npsol include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + NPSOL_OPTIONS *list = npsol_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + if (list->d_get_fcn) + { + double val = (npsol_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + } + else + { + int val = (npsol_opts.*list->i_get_fcn) (); + if (val < 0) + os << "depends on problem size"; + else + os << val; + } + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_npsol_option (const string& keyword, double val) +{ + NPSOL_OPTIONS *list = npsol_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_set_fcn) + (npsol_opts.*list->d_set_fcn) (val); + else + { + if (xisnan (val)) + { + error ("npsol_options: %s: expecting integer, found NaN", + keyword.c_str ()); + } + else + (npsol_opts.*list->i_set_fcn) (NINT (val)); + } + return; + } + list++; + } + + warning ("npsol_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_npsol_option (const string& keyword) +{ + octave_value retval; + + NPSOL_OPTIONS *list = npsol_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_get_fcn) + { + double val = (npsol_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + } + else + { + int val = (npsol_opts.*list->i_get_fcn) (); + if (val < 0) + retval = "depends on problem size"; + else + retval = static_cast<double> (val); + } + + return retval; + } + list++; + } + + warning ("npsol_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +#endif + +#if defined (NPSOL_MISSING) +DEFUN_DLD (npsol_options, , , + "This function requires NPSOL, which is not freely\n\ +redistributable. For more information, read the file\n\ +libcruft/npsol/README.MISSING in the source distribution.") +#else +DEFUN_DLD (npsol_options, args, , + "npsol_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for npsol. Keywords may be abbreviated\n\ +to the shortest match.") +#endif +{ + octave_value_list retval; + +#if defined (NPSOL_MISSING) + + print_usage ("npsol_options"); + +#else + + int nargin = args.length (); + + if (nargin == 0) + { + print_npsol_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_npsol_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_npsol_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("npsol_options"); + +#endif + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/pinv.cc @@ -0,0 +1,96 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (pinv, args, , + "pinv ( [, tol])\n\ +Returns the pseudoinverse of X; singular values less than tol are ignored.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2) + { + print_usage ("pinv"); + return retval; + } + + octave_value arg = args(0); + + double tol = 0.0; + if (nargin == 2) + tol = args(1).double_value (); + + if (error_state) + return retval; + + if (tol < 0.0) + { + error ("pinv: tol must be greater than zero"); + return retval; + } + + int arg_is_empty = empty_arg ("pinv", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return Matrix (); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + retval = m.pseudo_inverse (tol); + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + retval = m.pseudo_inverse (tol); + } + else + { + gripe_wrong_type_arg ("pinv", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/qpsol.cc @@ -0,0 +1,487 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iostream.h> + +#include "QPSOL.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "oct-sym.h" +#include "pager.h" +#include "utils.h" +#include "variables.h" + +#ifndef QPSOL_MISSING + +// XXX FIXME XXX -- this is duplicated in npsol.cc + +static int +linear_constraints_ok (const ColumnVector& x, const ColumnVector& llb, + const Matrix& c, const ColumnVector& lub, + char *warn_for, int warn) +{ + int x_len = x.capacity (); + int llb_len = llb.capacity (); + int lub_len = lub.capacity (); + int c_rows = c.rows (); + int c_cols = c.columns (); + + int ok = 1; + if (warn) + { + if (c_rows == 0 || c_cols == 0 || llb_len == 0 || lub_len == 0) + { + ok = 0; + error ("%s: linear constraints must have nonzero dimensions", + warn_for); + } + else if (x_len != c_cols || llb_len != lub_len || llb_len != c_rows) + { + ok = 0; + error ("%s: linear constraints have inconsistent dimensions", + warn_for); + } + } + + return ok; +} + +static QPSOL_options qpsol_opts; + +#endif + +#if defined (QPSOL_MISSING) +DEFUN_DLD (qpsol, , , + "This function requires QPSOL, which is not freely\n\ +redistributable. For more information, read the file\n\ +libcruft/qpsol/README.MISSING in the source distribution.") +#else +DEFUN_DLD (qpsol, args, nargout, + "[X, OBJ, INFO, LAMBDA] = qpsol (X, H, C [, LB, UB] [, A_LB, A, A_UB])\n\ +\n\ +Groups of arguments surrounded in `[]' are optional, but\n\ +must appear in the same relative order shown above.") +#endif +{ +/* + +Handle all of the following: + + 1. qpsol (x, H, c) + 2. qpsol (x, H, c, lb, ub) + 3. qpsol (x, H, c, lb, ub, llb, A, lub) + 4. qpsol (x, H, c, llb, A, lub) + +*/ + + octave_value_list retval; + +#if defined (QPSOL_MISSING) + + // Force a bad value of inform, and empty matrices for x, phi, and + // lambda. + + retval.resize (4, Matrix ()); + + retval(2) = -1.0; + + print_usage ("qpsol"); + +#else + + int nargin = args.length (); + + if (nargin < 3 || nargin == 4 || nargin == 7 || nargin > 8 + || nargout > 4) + { + print_usage ("qpsol"); + return retval; + } + + ColumnVector x = args(0).vector_value (); + + if (error_state || x.capacity () == 0) + { + error ("qpsol: expecting vector as first argument"); + return retval; + } + + Matrix H = args(1).matrix_value (); + + if (error_state || H.rows () != H.columns () || H.rows () != x.capacity ()) + { + error ("qpsol: H must be a square matrix consistent with the size of x"); + return retval; + } + + ColumnVector c = args(2).vector_value (); + + if (error_state || c.capacity () != x.capacity ()) + { + error ("qpsol: c must be a vector the same size as x"); + return retval; + } + + Bounds bounds; + if (nargin == 5 || nargin == 8) + { + ColumnVector lb = args(3).vector_value (); + ColumnVector ub = args(4).vector_value (); + + int lb_len = lb.capacity (); + int ub_len = ub.capacity (); + + if (error_state || lb_len != ub_len || lb_len != x.capacity ()) + { + error ("qpsol: lower and upper bounds and decision variable vector"); + error ("must all have the same number of elements"); + return retval; + } + + bounds.resize (lb_len); + bounds.set_lower_bounds (lb); + bounds.set_upper_bounds (ub); + } + + ColumnVector soln; + double objf; + ColumnVector lambda; + int inform; + + if (nargin == 3) + { + // 1. qpsol (x, H, c) + + QPSOL qp (x, H, c); + qp.set_options (qpsol_opts); + soln = qp.minimize (objf, inform, lambda); + + goto solved; + } + + if (nargin == 5) + { + // 2. qpsol (x, H, c, lb, ub) + + QPSOL qp (x, H, c, bounds); + qp.set_options (qpsol_opts); + soln = qp.minimize (objf, inform, lambda); + + goto solved; + } + + if (nargin == 6 || nargin == 8) + { + ColumnVector lub = args(nargin-1).vector_value (); + ColumnVector llb = args(nargin-3).vector_value (); + + if (error_state || llb.capacity () == 0 || lub.capacity () == 0) + { + error ("qpsol: bounds for linear constraints must be vectors"); + return retval; + } + + Matrix A = args(nargin-2).matrix_value (); + + if (error_state) + { + error ("qpsol: invalid linear constraint matrix"); + return retval; + } + + if (! linear_constraints_ok (x, llb, A, lub, "qpsol", 1)) + return retval; + + LinConst linear_constraints (llb, A, lub); + + if (nargin == 8) + { + // 3. qpsol (x, H, c, lb, ub, llb, A, lub) + + QPSOL qp (x, H, c, bounds, linear_constraints); + qp.set_options (qpsol_opts); + soln = qp.minimize (objf, inform, lambda); + } + else + { + // 4. qpsol (x, H, c, llb, A, lub) + + QPSOL qp (x, H, c, linear_constraints); + qp.set_options (qpsol_opts); + soln = qp.minimize (objf, inform, lambda); + } + goto solved; + } + + return retval; + + solved: + + retval.resize (nargout ? nargout : 1); + retval(0) = soln, 1; + if (nargout > 1) + retval(1) = objf; + if (nargout > 2) + retval(2) = static_cast<double> (inform); + if (nargout > 3) + retval(3) = lambda; + +#endif + + return retval; +} + +#ifndef QPSOL_MISSING + +typedef void (QPSOL_options::*d_set_opt_mf) (double); +typedef void (QPSOL_options::*i_set_opt_mf) (int); +typedef double (QPSOL_options::*d_get_opt_mf) (void); +typedef int (QPSOL_options::*i_get_opt_mf) (void); + +#define MAX_TOKENS 2 + +struct QPSOL_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + i_set_opt_mf i_set_fcn; + d_get_opt_mf d_get_fcn; + i_get_opt_mf i_get_fcn; +}; + +static QPSOL_OPTIONS qpsol_option_table [] = +{ + { "feasibility tolerance", + { "feasibility", "tolerance", 0, }, + { 1, 0, 0, }, 1, + QPSOL_options::set_feasibility_tolerance, 0, + QPSOL_options::feasibility_tolerance, 0, }, + + { "infinite bound", + { "infinite", "bound", 0, }, + { 2, 0, 0, }, 1, + QPSOL_options::set_infinite_bound, 0, + QPSOL_options::infinite_bound, 0, }, + + { "iteration limit", + { "iteration", "limit", 0, }, + { 2, 0, 0, }, 1, + 0, QPSOL_options::set_iteration_limit, + 0, QPSOL_options::iteration_limit, }, + + { "print level", + { "print", "level", 0, }, + { 1, 0, 0, }, 1, + 0, QPSOL_options::set_print_level, + 0, QPSOL_options::print_level, }, + + { 0, + { 0, 0, 0, }, + { 0, 0, 0, }, 0, + 0, 0, 0, 0, }, +}; + +static void +print_qpsol_option_list (ostream& os) +{ + print_usage ("qpsol_options", 1); + + os << "\n" + << "Options for qpsol include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + QPSOL_OPTIONS *list = qpsol_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + if (list->d_get_fcn) + { + double val = (qpsol_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + } + else + { + int val = (qpsol_opts.*list->i_get_fcn) (); + if (val < 0) + os << "depends on problem size"; + else + os << val; + } + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_qpsol_option (const string& keyword, double val) +{ + QPSOL_OPTIONS *list = qpsol_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_set_fcn) + (qpsol_opts.*list->d_set_fcn) (val); + else + { + if (xisnan (val)) + { + error ("qpsol_options: %s: expecting integer, found NaN", + keyword.c_str ()); + } + else + (qpsol_opts.*list->i_set_fcn) (NINT (val)); + } + return; + } + list++; + } + + warning ("qpsol_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_qpsol_option (const string& keyword) +{ + octave_value retval; + + QPSOL_OPTIONS *list = qpsol_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + if (list->d_get_fcn) + { + double val = (qpsol_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + } + else + { + int val = (qpsol_opts.*list->i_get_fcn) (); + if (val < 0) + retval = "depends on problem size"; + else + retval = static_cast<double> (val); + } + + return retval; + } + list++; + } + + warning ("qpsol_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +#endif + +#if defined (QPSOL_MISSING) +DEFUN_DLD (qpsol_options, , , + "This function requires QPSOL, which is not freely\n\ +redistributable. For more information, read the file\n\ +libcruft/qpsol/README.MISSING in the source distribution.") +#else +DEFUN_DLD (qpsol_options, args, , + "qpsol_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for qpsol. Keywords may be abbreviated\n\ +to the shortest match.") +#endif +{ + octave_value_list retval; + +#if defined (QPSOL_MISSING) + + print_usage ("qpsol"); + +#else + + int nargin = args.length (); + + if (nargin == 0) + { + print_qpsol_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_qpsol_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_qpsol_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("qpsol_options"); + +#endif + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/qr.cc @@ -0,0 +1,133 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "CmplxQR.h" +#include "CmplxQRP.h" +#include "dbleQR.h" +#include "dbleQRP.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (qr, args, nargout, + "[Q, R] = qr (X): form Q unitary and R upper triangular such\n\ + that Q * R = X\n\ +\n\ +[Q, R] = qr (X, 0): form the economy decomposition such that if X is\n\ + m by n then only the first n columns of Q are\n\ + computed.\n\ +\n\ +[Q, R, P] = qr (X): form QRP factorization of X where\n\ + P is a permutation matrix such that\n\ + A * P = Q * R\n\ +\n\ +[Q, R, P] = qr (X, 0): form the economy decomposition with \n\ + permutation vector P such that Q * R = X (:, P)\n\ +\n\ +qr (X) alone returns the output of the LAPACK routine dgeqrf, such\n\ +that R = triu (qr (X))") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1 && nargin != 2 || nargout > 3) + { + print_usage ("qr"); + return retval; + } + + octave_value arg = args(0); + + int arg_is_empty = empty_arg ("qr", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (3, Matrix ()); + + QR::type type = nargout == 1 ? QR::raw + : (nargin == 2 ? QR::economy : QR::std); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + if (nargout < 3) + { + QR fact (m, type); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + else + { + QRP fact (m, type); + retval(2) = fact.P (); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + if (nargout < 3) + { + ComplexQR fact (m, type); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + else + { + ComplexQRP fact (m, type); + retval(2) = fact.P (); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + } + } + else + { + gripe_wrong_type_arg ("qr", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/quad.cc @@ -0,0 +1,406 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include <iostream.h> + +#include "Quad.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-sym.h" +#include "pager.h" +#include "oct-obj.h" +#include "utils.h" +#include "variables.h" + +#if defined (quad) +#undef quad +#endif + +// Global pointer for user defined function required by quadrature functions. +static octave_symbol *quad_fcn; + +static Quad_options quad_opts; + +double +quad_user_function (double x) +{ + double retval = 0.0; + + octave_value_list args; + args(0) = x; + + if (quad_fcn) + { + octave_value_list tmp = quad_fcn->eval (1, args); + + if (error_state) + { + quad_integration_error = 1; // XXX FIXME XXX + gripe_user_supplied_eval ("quad"); + return retval; + } + + if (tmp.length () && tmp(0).is_defined ()) + { + retval = tmp(0).double_value (); + + if (error_state) + { + quad_integration_error = 1; // XXX FIXME XXX + gripe_user_supplied_eval ("quad"); + } + } + else + { + quad_integration_error = 1; // XXX FIXME XXX + gripe_user_supplied_eval ("quad"); + } + } + + return retval; +} + +DEFUN_DLD (quad, args, nargout, + "[V, IER, NFUN] = quad (F, A, B [, TOL] [, SING])\n\ +\n\ +Where the first argument is the name of the function to call to\n\ +compute the value of the integrand. It must have the form\n\ +\n\ + y = f (x)\n\ +\n\ +where y and x are scalars.\n\ +\n\ +The second and third arguments are limits of integration. Either or\n\ +both may be infinite.\n\ +\n\ +The optional argument tol is a vector that specifies the desired\n\ +accuracy of the result. The first element of the vector is the desired\n\ +absolute tolerance, and the second element is the desired relative\n\ +tolerance.\n\ +\n\ +The optional argument @var{sing} is a vector of values at which the\n\ +integrand is singular.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 3 || nargin > 5 || nargout > 4) + { + print_usage ("quad"); + return retval; + } + + quad_fcn = extract_function (args(0), "quad", "__quad_fcn__", + "function y = __quad_fcn__ (x) y = ", + "; endfunction"); + if (! quad_fcn) + return retval; + + double a = args(1).double_value (); + + if (error_state) + { + error ("quad: expecting second argument to be a scalar"); + return retval; + } + + double b = args(2).double_value (); + + if (error_state) + { + error ("quad: expecting third argument to be a scalar"); + return retval; + } + + int indefinite = 0; + IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; + double bound = 0.0; + if (xisinf (a) && xisinf (b)) + { + indefinite = 1; + indef_type = IndefQuad::doubly_infinite; + } + else if (xisinf (a)) + { + indefinite = 1; + bound = b; + indef_type = IndefQuad::neg_inf_to_bound; + } + else if (xisinf (b)) + { + indefinite = 1; + bound = a; + indef_type = IndefQuad::bound_to_inf; + } + + int ier = 0; + int nfun = 0; + double abserr = 0.0; + double val = 0.0; + double abstol = 1e-6; + double reltol = 1e-6; + ColumnVector tol (2); + ColumnVector sing; + int have_sing = 0; + switch (nargin) + { + case 5: + if (indefinite) + { + error("quad: singularities not allowed on infinite intervals"); + return retval; + } + + have_sing = 1; + + sing = args(4).vector_value (); + + if (error_state) + { + error ("quad: expecting vector of singularities as fourth argument"); + return retval; + } + + case 4: + tol = args(3).vector_value (); + + if (error_state) + { + error ("quad: expecting vector of tolerances as fifth argument"); + return retval; + } + + switch (tol.capacity ()) + { + case 2: + reltol = tol (1); + + case 1: + abstol = tol (0); + break; + + default: + error ("quad: expecting tol to contain no more than two values"); + return retval; + } + + case 3: + if (indefinite) + { + IndefQuad iq (quad_user_function, bound, indef_type, abstol, reltol); + iq.set_options (quad_opts); + val = iq.integrate (ier, nfun, abserr); + } + else + { + if (have_sing) + { + DefQuad dq (quad_user_function, a, b, sing, abstol, reltol); + dq.set_options (quad_opts); + val = dq.integrate (ier, nfun, abserr); + } + else + { + DefQuad dq (quad_user_function, a, b, abstol, reltol); + dq.set_options (quad_opts); + val = dq.integrate (ier, nfun, abserr); + } + } + break; + + default: + panic_impossible (); + break; + } + + retval(3) = abserr; + retval(2) = static_cast<double> (nfun); + retval(1) = static_cast<double> (ier); + retval(0) = val; + + return retval; +} + +typedef void (Quad_options::*d_set_opt_mf) (double); +typedef double (Quad_options::*d_get_opt_mf) (void); + +#define MAX_TOKENS 2 + +struct QUAD_OPTIONS +{ + const char *keyword; + const char *kw_tok[MAX_TOKENS + 1]; + int min_len[MAX_TOKENS + 1]; + int min_toks_to_match; + d_set_opt_mf d_set_fcn; + d_get_opt_mf d_get_fcn; +}; + +static QUAD_OPTIONS quad_option_table [] = +{ + { "absolute tolerance", + { "absolute", "tolerance", 0, }, + { 1, 0, 0, }, 1, + Quad_options::set_absolute_tolerance, + Quad_options::absolute_tolerance, }, + + { "relative tolerance", + { "relative", "tolerance", 0, }, + { 1, 0, 0, }, 1, + Quad_options::set_relative_tolerance, + Quad_options::relative_tolerance, }, + + { 0, + { 0, 0, 0, }, + { 0, 0, 0, }, 0, + 0, 0, }, +}; + +static void +print_quad_option_list (ostream& os) +{ + print_usage ("quad_options", 1); + + os << "\n" + << "Options for quad include:\n\n" + << " keyword value\n" + << " ------- -----\n\n"; + + QUAD_OPTIONS *list = quad_option_table; + + const char *keyword; + while ((keyword = list->keyword) != 0) + { + os.form (" %-40s ", keyword); + + double val = (quad_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + + os << "\n"; + list++; + } + + os << "\n"; +} + +static void +set_quad_option (const string& keyword, double val) +{ + QUAD_OPTIONS *list = quad_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + (quad_opts.*list->d_set_fcn) (val); + + return; + } + list++; + } + + warning ("quad_options: no match for `%s'", keyword.c_str ()); +} + +static octave_value_list +show_quad_option (const string& keyword) +{ + octave_value retval; + + QUAD_OPTIONS *list = quad_option_table; + + while (list->keyword != 0) + { + if (keyword_almost_match (list->kw_tok, list->min_len, keyword, + list->min_toks_to_match, MAX_TOKENS)) + { + return (quad_opts.*list->d_get_fcn) (); + } + list++; + } + + warning ("quad_options: no match for `%s'", keyword.c_str ()); + + return retval; +} + +DEFUN_DLD (quad_options, args, , + "quad_options (KEYWORD, VALUE)\n\ +\n\ +Set or show options for quad. Keywords may be abbreviated\n\ +to the shortest match.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin == 0) + { + print_quad_option_list (octave_stdout); + return retval; + } + else if (nargin == 1 || nargin == 2) + { + string keyword = args(0).string_value (); + + if (! error_state) + { + if (nargin == 1) + return show_quad_option (keyword); + else + { + double val = args(1).double_value (); + + if (! error_state) + { + set_quad_option (keyword, val); + return retval; + } + } + } + } + + print_usage ("quad_options"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/qzval.cc @@ -0,0 +1,73 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Written by A. S. Hodel <scotte@eng.auburn.edu> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <cfloat> + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" + +DEFUN_DLD (qzval, args, , + "X = qzval (A, B)\n\ +\n\ +compute generalized eigenvalues of the matrix pencil (A - lambda B).\n\ +A and B must be real matrices.") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 2) + { + octave_value arg_a = args(0); + octave_value arg_b = args(1); + + Matrix a = arg_a.matrix_value (); + Matrix b = arg_b.matrix_value (); + + if (! error_state) + { + ComplexColumnVector tmp = Qzval (a, b); + + if (! error_state) + retval = tmp; + } + } + else + print_usage ("qzval"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/rand.cc @@ -0,0 +1,410 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <ctime> + +#include <string> + +#include "f77-fcn.h" +#include "lo-mappers.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "unwind-prot.h" +#include "utils.h" + +// Possible distributions of random numbers. This was handled with an +// enum, but unwind_protecting that doesn't work so well. +#define uniform_dist 1 +#define normal_dist 2 + +// Current distribution of random numbers. +static int current_distribution = uniform_dist; + +// Has the seed been set yet? +static int initialized = 0; + +extern "C" +{ + int *F77_FCN (dgennor, DGENNOR) (const double&, const double&, + double&); + + int *F77_FCN (dgenunf, DGENUNF) (const double&, const double&, + double&); + + int *F77_FCN (setall, SETALL) (const int&, const int&); + + int *F77_FCN (getsd, GETSD) (int&, int&); + + int *F77_FCN (setsd, SETSD) (const int&, const int&); + + int *F77_FCN (setcgn, SETCGN) (const int&); +} + +static double +curr_rand_seed (void) +{ + union d2i { double d; int i[2]; }; + union d2i u; + F77_FCN (getsd, GETSD) (u.i[0], u.i[1]); + return u.d; +} + +static int +force_to_fit_range (int i, int lo, int hi) +{ + assert (hi > lo && lo >= 0 && hi > lo); + + i = i > 0 ? i : -i; + + if (i < lo) + i = lo; + else if (i > hi) + i = i % hi; + + return i; +} + +static void +set_rand_seed (double val) +{ + union d2i { double d; int i[2]; }; + union d2i u; + u.d = val; + int i0 = force_to_fit_range (u.i[0], 1, 2147483563); + int i1 = force_to_fit_range (u.i[1], 1, 2147483399); + F77_FCN (setsd, SETSD) (i0, i1); +} + +static char * +curr_rand_dist (void) +{ + if (current_distribution == uniform_dist) + return "uniform"; + else if (current_distribution == normal_dist) + return "normal"; + else + { + panic_impossible (); + return 0; + } +} + +// Make the random number generator give us a different sequence every +// time we start octave unless we specifically set the seed. The +// technique used below will cycle monthly, but it it does seem to +// work ok to give fairly different seeds each time Octave starts. + +static void +do_initialization (void) +{ + time_t now; + struct tm *tm; + + time (&now); + tm = localtime (&now); + + int hour = tm->tm_hour + 1; + int minute = tm->tm_min + 1; + int second = tm->tm_sec + 1; + + int s0 = tm->tm_mday * hour * minute * second; + int s1 = hour * minute * second; + + s0 = force_to_fit_range (s0, 1, 2147483563); + s1 = force_to_fit_range (s1, 1, 2147483399); + + F77_FCN (setall, SETALL) (s0, s1); + + initialized = 1; +} + +static octave_value_list +do_rand (const octave_value_list& args, int nargin) +{ + octave_value_list retval; + + int n = 0; + int m = 0; + + if (nargin == 0) + { + n = 1; + m = 1; + + goto gen_matrix; + } + else if (nargin == 1) + { + octave_value tmp = args(0); + + if (tmp.is_string ()) + { + string s_arg = tmp.string_value (); + + if (s_arg == "dist") + { + retval(0) = curr_rand_dist (); + } + else if (s_arg == "seed") + { + retval(0) = curr_rand_seed (); + } + else if (s_arg == "uniform") + { + current_distribution = uniform_dist; + + F77_FCN (setcgn, SETCGN) (uniform_dist); + } + else if (s_arg == "normal") + { + current_distribution = normal_dist; + + F77_FCN (setcgn, SETCGN) (normal_dist); + } + else + error ("rand: unrecognized string argument"); + } + else if (tmp.is_scalar_type ()) + { + double dval = tmp.double_value (); + + if (xisnan (dval)) + { + error ("rand: NaN is invalid a matrix dimension"); + } + else + { + m = n = NINT (tmp.double_value ()); + + if (! error_state) + goto gen_matrix; + } + } + else if (tmp.is_range ()) + { + Range r = tmp.range_value (); + n = 1; + m = r.nelem (); + goto gen_matrix; + } + else if (tmp.is_matrix_type ()) + { + // XXX FIXME XXX -- this should probably use the function + // from data.cc. + + Matrix a = args(0).matrix_value (); + + if (error_state) + return retval; + + n = a.rows (); + m = a.columns (); + + if (n == 1 && m == 2) + { + n = NINT (a (0, 0)); + m = NINT (a (0, 1)); + } + else if (n == 2 && m == 1) + { + n = NINT (a (0, 0)); + m = NINT (a (1, 0)); + } + else + warning ("rand (A): use rand (size (A)) instead"); + + goto gen_matrix; + } + else + { + gripe_wrong_type_arg ("rand", tmp); + return retval; + } + } + else if (nargin == 2) + { + if (args(0).is_string ()) + { + if (args(0).string_value () == "seed") + { + double d = args(1).double_value (); + + if (! error_state) + set_rand_seed (d); + } + } + else + { + double dval = args(0).double_value (); + + if (xisnan (dval)) + { + error ("rand: NaN is invalid as a matrix dimension"); + } + else + { + n = NINT (dval); + + if (! error_state) + { + m = NINT (args(1).double_value ()); + + if (! error_state) + goto gen_matrix; + } + } + } + } + + return retval; + + gen_matrix: + + if (n == 0 || m == 0) + { + Matrix m; + retval.resize (1, m); + } + else if (n > 0 && m > 0) + { + Matrix rand_mat (n, m); + for (int j = 0; j < m; j++) + for (int i = 0; i < n; i++) + { + double val; + switch (current_distribution) + { + case uniform_dist: + F77_FCN (dgenunf, DGENUNF) (0.0, 1.0, val); + rand_mat (i, j) = val; + break; + + case normal_dist: + F77_FCN (dgennor, DGENNOR) (0.0, 1.0, val); + rand_mat (i, j) = val; + break; + + default: + panic_impossible (); + break; + } + } + + retval(0) = rand_mat; + } + else + error ("rand: invalid negative argument"); + + return retval; +} + +DEFUN_DLD (rand, args, nargout, + "rand -- generate a random value from a uniform distribution\n\ +\n\ +rand (N) -- generate N x N matrix\n\ +rand (size (A)) -- generate matrix the size of A\n\ +rand (N, M) -- generate N x M matrix\n\ +rand (SEED) -- get current seed\n\ +rand (SEED, N) -- set seed\n\ +\n\ +See also: randn") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin > 2 || nargout > 1) + print_usage ("rand"); + else + { + if (! initialized) + do_initialization (); + + retval = do_rand (args, nargin); + } + + return retval; +} + +static void +reset_rand_generator (void *) +{ + F77_FCN (setcgn, SETCGN) (current_distribution); +} + +DEFUN_DLD (randn, args, nargout, + "randn -- generate a random value from a normal distribution\n\ +\n\ +randn (N) -- generate N x N matrix\n\ +randn (size (A)) -- generate matrix the size of A\n\ +randn (N, M) -- generate N x M matrix\n\ +randn (SEED) -- get current seed\n\ +randn (SEED, N) -- set seed\n\ +\n\ +See also: rand") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin > 2 || nargout > 1) + print_usage ("randn"); + else + { + if (! initialized) + do_initialization (); + + begin_unwind_frame ("randn"); + + // This relies on the fact that elements are popped from the + // unwind stack in the reverse of the order they are pushed + // (i.e. current_distribution will be reset before calling + // reset_rand_generator()). + + add_unwind_protect (reset_rand_generator, 0); + unwind_protect_int (current_distribution); + + current_distribution = normal_dist; + + F77_FCN (setcgn, SETCGN) (normal_dist); + + retval = do_rand (args, nargin); + + run_unwind_frame ("randn"); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/schur.cc @@ -0,0 +1,153 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include "CmplxSCHUR.h" +#include "dbleSCHUR.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (schur, args, nargout, + "[U, S] = schur (A) or S = schur (A)\n\ +\n\ +or, for ordered Schur:\n\ +\n\ + [U, S] = schur (A, TYPE) or S = schur (A, TYPE)\n\ +where TYPE is a string that begins with one of the following\n\ +characters:\n\ +\n\ + A = continuous time poles\n\ + D = discrete time poles\n\ + U = unordered schur (default)") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2 || nargout > 2) + { + print_usage ("schur"); + return retval; + } + + octave_value arg = args(0); + + string ord; + + if (nargin == 2) + { + ord = args(1).string_value (); + + if (error_state) + { + error ("schur: expecting string as second argument"); + return retval; + } + } + + char ord_char = ord.empty () ? 'U' : ord[0]; + + if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D' + && ord_char != 'u' && ord_char != 'a' && ord_char != 'd') + { + warning ("schur: incorrect ordered schur argument `%c'", + ord.c_str ()); + return retval; + } + + int nr = arg.rows (); + int nc = arg.columns (); + + int arg_is_empty = empty_arg ("schur", nr, nc); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (2, Matrix ()); + + if (nr != nc) + { + gripe_square_matrix_required ("schur"); + return retval; + } + + if (arg.is_real_type ()) + { + Matrix tmp = arg.matrix_value (); + + if (! error_state) + { + SCHUR result (tmp, ord); + + if (nargout == 0 || nargout == 1) + { + retval(0) = result.schur_matrix (); + } + else + { + retval(1) = result.schur_matrix (); + retval(0) = result.unitary_matrix (); + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix ctmp = arg.complex_matrix_value (); + + if (! error_state) + { + ComplexSCHUR result (ctmp, ord); + + if (nargout == 0 || nargout == 1) + { + retval(0) = result.schur_matrix (); + } + else + { + retval(1) = result.schur_matrix (); + retval(0) = result.unitary_matrix (); + } + } + } + else + { + gripe_wrong_type_arg ("schur", arg); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/sort.cc @@ -0,0 +1,382 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" + +// This is algorithm 5.2.4L from Knuth, Volume 3. + +// XXX FIXME XXX -- there is way too much duplicated code here given +// that the sort algorithms are all the same, and only the type of the +// data and the comparison changes... +// +// Maybe some cpp abuse will make it better. + +static Array<int> +create_index_array (int n) +{ + Array<int> l (n+2); + + l (0) = 1; + + for (int i = 1; i < n - 1; i++) + l (i) = -(i+2); + + l (n-1) = 0; + l (n) = 0; + l (n+1) = 2; + + return l; +} + +#define SORT_INIT_PHASE(n) \ + int s = 0; \ + int t = n + 1; \ + int p = l (s); \ + int q = l (t); \ + if (q == 0) \ + break + +#define SORT_COMMON_CODE \ + p = -p; \ + q = -q; \ + if (q == 0) \ + { \ + l (s) = (l (s) < 0) \ + ? ((p < 0) ? p : -p) \ + : ((p >= 0) ? p : -p); \ + l (t) = 0; \ + break; \ + } \ + +#define SORT_REORDER_PHASE_ONE \ + l (s) = (l (s) < 0) \ + ? ((q < 0) ? q : -q) \ + : ((q >= 0) ? q : -q); \ + s = q; \ + q = l (q); \ + if (q <= 0) \ + { \ + l (s) = p; \ + s = t; \ + do \ + { \ + t = p; \ + p = l (p); \ + } \ + while (p > 0); \ + SORT_COMMON_CODE; \ + } \ + +#define SORT_REORDER_PHASE_TWO \ + l (s) = (l (s) < 0) \ + ? ((p < 0) ? p : -p) \ + : ((p >= 0) ? p : -p); \ + s = p; \ + p = l (p); \ + if (p <= 0) \ + { \ + l (s) = q; \ + s = t; \ + do \ + { \ + t = q; \ + q = l (q); \ + } \ + while (q > 0); \ + SORT_COMMON_CODE; \ + } + +#define DO_SORT(n, condition) \ + while (1) \ + { \ + SORT_INIT_PHASE(n); \ + while (1) \ + { \ + if (condition) \ + { \ + SORT_REORDER_PHASE_ONE; \ + } \ + else \ + { \ + SORT_REORDER_PHASE_TWO; \ + } \ + } \ + } + +#define VECTOR_CREATE_RETURN_VALUES(vs, v) \ + int k = l (0); \ + idx (0) = k; \ + vs (0) = v (k-1); \ + for (int i = 1; i < n; i++) \ + { \ + k = l (static_cast<int> (idx (i-1))); \ + idx (i) = k; \ + vs (i) = v (k-1); \ + } + +#define MATRIX_CREATE_RETURN_VALUES(ms, m) \ + int k = l (0); \ + idx (0, j) = k; \ + ms (0, j) = m (k-1, j); \ + for (int i = 1; i < nr; i++) \ + { \ + k = l (static_cast<int> (idx (i-1, j))); \ + idx (i, j) = k; \ + ms (i, j) = m (k-1, j); \ + } + +static octave_value_list +mx_sort (const Matrix& m) +{ + octave_value_list retval; + + int nr = m.rows (); + int nc = m.columns (); + + Matrix ms (nr, nc); + Matrix idx (nr, nc); + + if (nr == 1 && nc > 0) + { + retval (1) = Matrix (nr, nc, 1.0); + retval (0) = m; + + return retval; + } + else if (nr > 1 && nc > 0) + { + for (int j = 0; j < nc; j++) + { + Array<int> l = create_index_array (nr); + + DO_SORT (nr, (m (p-1, j) > m (q-1, j))); + + MATRIX_CREATE_RETURN_VALUES (ms, m); + } + } + + retval (1) = idx; + retval (0) = ms; + + return retval; +} + +static octave_value_list +mx_sort (const RowVector& v) +{ + octave_value_list retval; + + int n = v.capacity (); + + RowVector vs (n); + RowVector idx (n); + + if (n == 1) + { + retval (1) = RowVector (n, 1.0); + retval (0) = v; + + return retval; + } + else if (n > 1) + { + Array<int> l = create_index_array (n); + + DO_SORT (n, (v (p-1) > v (q-1))); + + VECTOR_CREATE_RETURN_VALUES (vs, v); + } + + retval (1) = octave_value (idx, 0); + retval (0) = octave_value (vs, 0); + + return retval; +} + +static octave_value_list +mx_sort (const ComplexMatrix& cm) +{ + octave_value_list retval; + + int nr = cm.rows (); + int nc = cm.columns (); + + ComplexMatrix cms (nr, nc); + Matrix idx (nr, nc); + + if (nr == 1 && nc > 0) + { + retval (1) = Matrix (nr, nc, 1.0); + retval (0) = cm; + + return retval; + } + else if (nr > 1 && nc > 0) + { + for (int j = 0; j < nc; j++) + { + Array<int> l = create_index_array (nr); + + int all_elts_real = 1; + for (int i = 0; i < nr; i++) + if (imag (cm (i, j)) != 0.0) + { + all_elts_real = 0; + break; + } + + DO_SORT (nr, ((all_elts_real + && real (cm (p-1, j)) > real (cm (q-1, j))) + || abs (cm (p-1, j)) > abs (cm (q-1, j)))); + + MATRIX_CREATE_RETURN_VALUES (cms, cm); + } + } + + retval (1) = idx; + retval (0) = cms; + + return retval; +} + +static octave_value_list +mx_sort (ComplexRowVector& cv) +{ + octave_value_list retval; + + int n = cv.capacity (); + + ComplexRowVector cvs (n); + RowVector idx (n); + + if (n == 1) + { + retval (1) = RowVector (n, 1.0); + retval (0) = cv; + + return retval; + } + else if (n > 1) + { + Array<int> l = create_index_array (n); + + int all_elts_real = 1; + for (int i = 0; i < n; i++) + if (imag (cv (i)) != 0.0) + { + all_elts_real = 0; + break; + } + + DO_SORT (n, ((all_elts_real + && real (cv (p-1)) > real (cv (q-1))) + || abs (cv (p-1)) > abs (cv (q-1)))); + + VECTOR_CREATE_RETURN_VALUES (cvs, cv); + } + + retval (1) = octave_value (idx, 0); + retval (0) = octave_value (cvs, 0); + + return retval; +} + +DEFUN_DLD (sort, args, nargout, + "[S, I] = sort (X)\n\ +\n\ +sort the columns of X, optionally return sort index") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 1) + { + print_usage ("sort"); + return retval; + } + + int return_idx = nargout > 1; + if (return_idx) + retval.resize (2); + else + retval.resize (1); + + octave_value arg = args(0); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + if (m.rows () == 1) + { + int nc = m.columns (); + RowVector v (nc); + for (int i = 0; i < nc; i++) + v (i) = m (0, i); + + retval = mx_sort (v); + } + else + retval = mx_sort (m); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix cm = arg.complex_matrix_value (); + + if (! error_state) + { + if (cm.rows () == 1) + { + int nc = cm.columns (); + ComplexRowVector cv (nc); + for (int i = 0; i < nc; i++) + cv (i) = cm (0, i); + + retval = mx_sort (cv); + } + else + retval = mx_sort (cm); + } + } + else + gripe_wrong_type_arg ("sort", arg); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/svd.cc @@ -0,0 +1,141 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "CmplxSVD.h" +#include "dbleSVD.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "pr-output.h" +#include "utils.h" + +DEFUN_DLD (svd, args, nargout, + "S = svd (X) or [U, S, V] = svd (X [, 0])\n\ +\n\ +Compute the singular value decomposition of X. Given a second input\n\ +argument, an `economy' sized factorization is computed that omits\n\ +unnecessary rows and columns of U and V.\n\ +\n\ +X may not contain any Inf or NaN values.") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin < 1 || nargin > 2 || nargout == 2 || nargout > 3) + { + print_usage ("svd"); + return retval; + } + + octave_value arg = args(0); + + int arg_is_empty = empty_arg ("svd", arg.rows (), arg.columns ()); + + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (3, Matrix ()); + + SVD::type type = ((nargout == 0 || nargout == 1) + ? SVD::sigma_only + : (nargin == 2) ? SVD::economy : SVD::std); + + if (arg.is_real_type ()) + { + Matrix tmp = arg.matrix_value (); + + if (! error_state) + { + if (tmp.any_element_is_inf_or_nan ()) + { + error ("svd: cannot take SVD of matrix containing Inf or\ + NaN values"); + return retval; + } + + SVD result (tmp, type); + + DiagMatrix sigma = result.singular_values (); + + if (nargout == 0 || nargout == 1) + { + retval(0) = octave_value (sigma.diag (), 1); + } + else + { + retval(2) = result.right_singular_matrix (); + retval(1) = sigma; + retval(0) = result.left_singular_matrix (); + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix ctmp = arg.complex_matrix_value (); + + if (! error_state) + { + if (ctmp.any_element_is_inf_or_nan ()) + { + error ("svd: cannot take SVD of matrix containing Inf or\ + NaN values"); + return retval; + } + + ComplexSVD result (ctmp, type); + + DiagMatrix sigma = result.singular_values (); + + if (nargout == 0 || nargout == 1) + { + retval(0) = octave_value (sigma.diag (), 1); + } + else + { + retval(2) = result.right_singular_matrix (); + retval(1) = sigma; + retval(0) = result.left_singular_matrix (); + } + } + } + else + { + gripe_wrong_type_arg ("svd", arg); + return retval; + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/syl.cc @@ -0,0 +1,138 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +// Written by A. S. Hodel <scotte@eng.auburn.edu> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (syl, args, nargout, + "X = syl (A, B, C): solve the Sylvester equation A X + X B + C = 0") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 3 || nargout > 1) + { + print_usage ("syl"); + return retval; + } + + octave_value arg_a = args(0); + octave_value arg_b = args(1); + octave_value arg_c = args(2); + + int a_nr = arg_a.rows (); + int a_nc = arg_a.columns (); + + int b_nr = arg_b.rows (); + int b_nc = arg_b.columns (); + + int c_nr = arg_c.rows (); + int c_nc = arg_c.columns (); + + int arg_a_is_empty = empty_arg ("syl", a_nr, a_nc); + int arg_b_is_empty = empty_arg ("syl", b_nr, b_nc); + int arg_c_is_empty = empty_arg ("syl", c_nr, c_nc); + + if (arg_a_is_empty > 0 && arg_b_is_empty > 0 && arg_c_is_empty > 0) + return Matrix (); + else if (arg_a_is_empty || arg_b_is_empty || arg_c_is_empty) + return retval; + + // Arguments are not empty, so check for correct dimensions. + + if (a_nr != a_nc || b_nr != b_nc) + { + gripe_square_matrix_required ("syl: first two parameters:"); + return retval; + } + else if (a_nr != c_nr || b_nr != c_nc) + { + gripe_nonconformant (); + return retval; + } + + // Dimensions look o.k., let's solve the problem. + + if (arg_a.is_complex_type () + || arg_b.is_complex_type () + || arg_c.is_complex_type ()) + { + // Do everything in complex arithmetic; + + ComplexMatrix ca = arg_a.complex_matrix_value (); + + if (error_state) + return retval; + + ComplexMatrix cb = arg_b.complex_matrix_value (); + + if (error_state) + return retval; + + ComplexMatrix cc = arg_c.complex_matrix_value (); + + if (error_state) + return retval; + + retval = Sylvester (ca, cb, cc); + } + else + { + // Do everything in real arithmetic. + + Matrix ca = arg_a.matrix_value (); + + if (error_state) + return retval; + + Matrix cb = arg_b.matrix_value (); + + if (error_state) + return retval; + + Matrix cc = arg_c.matrix_value (); + + if (error_state) + return retval; + + retval = Sylvester (ca, cb, cc); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/time.cc @@ -0,0 +1,314 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <string> + +#include<iostream.h> + +#include "defun-dld.h" +#include "error.h" +#include "help.h" +#include "oct-map.h" +#include "systime.h" +#include "ov.h" +#include "oct-obj.h" +#include "utils.h" + +// Date and time functions. + +static Octave_map +mk_tm_map (struct tm *tm, double fraction) +{ + Octave_map m; + + m ["usec"] = fraction * 1e6; + m ["sec"] = static_cast<double> (tm->tm_sec); + m ["min"] = static_cast<double> (tm->tm_min); + m ["hour"] = static_cast<double> (tm->tm_hour); + m ["mday"] = static_cast<double> (tm->tm_mday); + m ["mon"] = static_cast<double> (tm->tm_mon); + m ["year"] = static_cast<double> (tm->tm_year); + m ["wday"] = static_cast<double> (tm->tm_wday); + m ["yday"] = static_cast<double> (tm->tm_yday); + m ["isdst"] = static_cast<double> (tm->tm_isdst); + +#if defined (HAVE_TM_ZONE) + m ["zone"] = tm->tm_zone; +#elif defined (HAVE_TZNAME) + if (tm->tm_isdst == 0 || tm->tm_isdst == 1) + m ["zone"] = tzname[tm->tm_isdst]; +#endif + + return m; +} + +static struct tm* +extract_tm (Octave_map &m, double& fraction) +{ + static struct tm tm; + + fraction = (m ["usec"] . double_value ()) / 1e6; + tm.tm_sec = static_cast<int> (m ["sec"] . double_value ()); + tm.tm_min = static_cast<int> (m ["min"] . double_value ()); + tm.tm_hour = static_cast<int> (m ["hour"] . double_value ()); + tm.tm_mday = static_cast<int> (m ["mday"] . double_value ()); + tm.tm_mon = static_cast<int> (m ["mon"] . double_value ()); + tm.tm_year = static_cast<int> (m ["year"] . double_value ()); + tm.tm_wday = static_cast<int> (m ["wday"] . double_value ()); + tm.tm_yday = static_cast<int> (m ["yday"] . double_value ()); + tm.tm_isdst = static_cast<int> (m ["isdst"] . double_value ()); + +#if defined (HAVE_TM_ZONE) + string tstr = m ["zone"] . string_value (); + tm.tm_zone = tstr.c_str (); +#endif + + return &tm; +} + +DEFUN_DLD (time, , , + "time ()\n\ +\n\ +Return current time. On Unix systems, this is the number of\n\ +seconds since the epoch.") +{ + time_t now; + double fraction = 0.0; + +#if defined (HAVE_GETTIMEOFDAY) + + struct timeval tp; + +#if defined (GETTIMEOFDAY_NO_TZ) + gettimeofday (&tp); +#else + gettimeofday (&tp, 0); +#endif + + now = tp.tv_sec; + + fraction = tp.tv_usec / 1e6; + +#else + + now = time (0); + +#endif + + return static_cast<double> (now) + fraction; +} + +DEFUN_DLD (gmtime, args, , + "gmtime (TIME)\n\ +\n\ +Given a value returned from time(), return a structure like that\n\ +returned from localtime() but with values corresponding to\n\ +Coordinated Universal Time (UTC).") +{ + octave_value_list retval; + + if (args.length () == 1) + { + double tmp = args(0).double_value (); + + if (! error_state) + { + time_t timeval = static_cast<int> (tmp); + double ip; + double fraction = modf (tmp, &ip); + + retval = octave_value (mk_tm_map (gmtime (&timeval), fraction)); + } + } + else + print_usage ("gmtime"); + + return retval; +} + +DEFUN_DLD (localtime, args, , + "localtime (TIME)\n\ +\n\ +Given a value returned from time(), return a structure with\n\ +the following elements:\n\ +\n\ + usec : microseconds after the second (0, 999999)\n\ + sec : seconds after the minute (0, 61)\n\ + min : minutes after the hour (0, 59)\n\ + hour : hours since midnight (0, 23)\n\ + mday : day of the month (1, 31)\n\ + mon : months since January (0, 11)\n\ + year : years since 1900\n\ + wday : days since Sunday (0, 6)\n\ + yday : days since January 1 (0, 365)\n\ + isdst : Daylight Savings Time flag\n\ + zone : Time zone") +{ + octave_value_list retval; + + if (args.length () == 1) + { + double tmp = args(0).double_value (); + + if (! error_state) + { + time_t timeval = static_cast<int> (tmp); + double ip; + double fraction = modf (tmp, &ip); + + retval = octave_value (mk_tm_map (localtime (&timeval), fraction)); + } + } + else + print_usage ("localtime"); + + return retval; +} + +DEFUN_DLD (mktime, args, , + "mktime (TMSTRUCT)") +{ + octave_value_list retval; + + if (args.length () == 1 && args(0).is_map ()) + { + Octave_map map = args(0).map_value (); + + double fraction; + + struct tm *tm = extract_tm (map, fraction); + + if (! error_state) + retval = static_cast<double> (mktime (tm)) + fraction; + } + else + print_usage ("mktime"); + + return retval; +} + +DEFUN_DLD (strftime, args, , + "strftime (FMT, TMSTRUCT)\n\ +\n\ +Performs `%' substitutions similar to those in printf. Except where\n\ +noted, substituted fields have a fixed size; numeric fields are\n\ +padded if necessary. Padding is with zeros by default; for fields\n\ +that display a single number, padding can be changed or inhibited by\n\ +following the `%' with one of the modifiers described below.\n\ +Unknown field specifiers are copied as normal characters. All other\n\ +characters are copied to the output without change.\n\ +\n\ +Supports a superset of the ANSI C field specifiers.\n\ +\n\ +Literal character fields:\n\ +\n\ + % %\n\ + n newline\n\ + t tab\n\ +\n\ +Numeric modifiers (a nonstandard extension):\n\ +\n\ + - do not pad the field\n\ + _ pad the field with spaces\n\ +\n\ +Time fields:\n\ +\n\ + %H hour (00..23)\n\ + %I hour (01..12)\n\ + %k hour ( 0..23)\n\ + %l hour ( 1..12)\n\ + %M minute (00..59)\n\ + %p locale's AM or PM\n\ + %r time, 12-hour (hh:mm:ss [AP]M)\n\ + %R time, 24-hour (hh:mm)\n\ + %s time in seconds since 00:00:00, Jan 1, 1970 (a nonstandard extension)\n\ + %S second (00..61)\n\ + %T time, 24-hour (hh:mm:ss)\n\ + %X locale's time representation (%H:%M:%S)\n\ + %Z time zone (EDT), or nothing if no time zone is determinable\n\ + %z offset from GMT\n\ +\n\ +Date fields:\n\ +\n\ + %a locale's abbreviated weekday name (Sun..Sat)\n\ + %A locale's full weekday name, variable length (Sunday..Saturday)\n\ + %b locale's abbreviated month name (Jan..Dec)\n\ + %B locale's full month name, variable length (January..December)\n\ + %c locale's date and time (Sat Nov 04 12:02:33 EST 1989)\n\ + %C century (00..99)\n\ + %d day of month (01..31)\n\ + %e day of month ( 1..31)\n\ + %D date (mm/dd/yy)\n\ + %h same as %b\n\ + %j day of year (001..366)\n\ + %m month (01..12)\n\ + %U week number of year with Sunday as first day of week (00..53)\n\ + %w day of week (0..6)\n\ + %W week number of year with Monday as first day of week (00..53)\n\ + %x locale's date representation (mm/dd/yy)\n\ + %y last two digits of year (00..99)\n\ + %Y year (1970...)") +{ + octave_value_list retval; + + if (args.length () == 2 && args(0).is_string () && args(1).is_map ()) + { + string fmt = args(0).string_value (); + + Octave_map map = args(1).map_value (); + + double fraction; + + struct tm *tm = extract_tm (map, fraction); + + if (! error_state) + { + const char *fmt_str = fmt.c_str (); + + size_t bufsize = strftime (0, (size_t) UINT_MAX, fmt_str, tm); + + char *buf = new char [++bufsize]; + + buf[0] = '\0'; + + strftime (buf, bufsize, fmt_str, tm); + + retval = buf; + + delete [] buf; + } + } + else + print_usage ("strftime"); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
deleted file mode 100644 --- a/src/Map-fnc.cc +++ /dev/null @@ -1,49 +0,0 @@ -/* - -Copyright (C) 1996, 1997 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Instantiate Maps of file_name_cache_elts. - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <string> - -#include "Map.h" -#include "Map.cc" - -#include "str-vec.h" - -#include "fn-cache.h" - -template class Map<file_name_cache_elt>; -template class CHNode<file_name_cache_elt>; -template class CHMap<file_name_cache_elt>; - -template static int goodCHptr (CHNode<file_name_cache_elt> *t); -template static int CHptr_to_index (CHNode<file_name_cache_elt> *t); - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
deleted file mode 100644 --- a/src/Map-tc.cc +++ /dev/null @@ -1,45 +0,0 @@ -/* - -Copyright (C) 1996, 1997 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Instantiate Maps of octave_values. - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "Map.h" -#include "Map.cc" - -#include "ov.h" - -template class Map<octave_value>; -template class CHNode<octave_value>; -template class CHMap<octave_value>; - -template static int goodCHptr (CHNode<octave_value> *t); -template static int CHptr_to_index (CHNode<octave_value> *t); - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-b-b.cc @@ -0,0 +1,55 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-bool.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// bool by bool ops. + +DEFBINOP_OP (eq, bool, bool, ==) +DEFBINOP_OP (ne, bool, bool, !=) + +void +install_b_b_ops (void) +{ + INSTALL_BINOP (eq, octave_bool, octave_bool, eq); + INSTALL_BINOP (ne, octave_bool, octave_bool, ne); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-bm-bm.cc @@ -0,0 +1,55 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-bool-mat.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// bool matrix by bool matrix ops. + +DEFBINOP_OP (eq, bool_matrix, bool_matrix, ==) +DEFBINOP_OP (ne, bool_matrix, bool_matrix, !=) + +void +install_bm_bm_ops (void) +{ + INSTALL_BINOP (eq, octave_bool_matrix, octave_bool_matrix, eq); + INSTALL_BINOP (ne, octave_bool_matrix, octave_bool_matrix, ne); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cm-cm.cc @@ -0,0 +1,107 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-cx-mat.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex matrix by complex matrix ops. + +DEFBINOP_OP (add, complex_matrix, complex_matrix, +) +DEFBINOP_OP (sub, complex_matrix, complex_matrix, -) +DEFBINOP_OP (mul, complex_matrix, complex_matrix, *) +DEFBINOP_FN (div, complex_matrix, complex_matrix, xdiv) + +DEFBINOPX (pow, complex_matrix, complex_matrix) +{ + error ("can't do A ^ B for A and B both matrices"); + return octave_value (); +} + +DEFBINOP_FN (ldiv, complex_matrix, complex_matrix, xleftdiv) + +DEFBINOP_FN (lt, complex_matrix, complex_matrix, mx_el_lt) +DEFBINOP_FN (le, complex_matrix, complex_matrix, mx_el_le) +DEFBINOP_FN (eq, complex_matrix, complex_matrix, mx_el_eq) +DEFBINOP_FN (ge, complex_matrix, complex_matrix, mx_el_ge) +DEFBINOP_FN (gt, complex_matrix, complex_matrix, mx_el_gt) +DEFBINOP_FN (ne, complex_matrix, complex_matrix, mx_el_ne) + +DEFBINOP_FN (el_mul, complex_matrix, complex_matrix, product) +DEFBINOP_FN (el_div, complex_matrix, complex_matrix, quotient) +DEFBINOP_FN (el_pow, complex_matrix, complex_matrix, elem_xpow) + +DEFBINOP (el_ldiv, complex_matrix, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_complex_matrix&); + + return octave_value (quotient (v2.complex_matrix_value (), + v1.complex_matrix_value ())); +} + +DEFBINOP_FN (el_and, complex_matrix, complex_matrix, mx_el_and) +DEFBINOP_FN (el_or, complex_matrix, complex_matrix, mx_el_or) + +DEFASSIGNOP_FN (assign, complex_matrix, complex_matrix, assign) + +void +install_cm_cm_ops (void) +{ + INSTALL_BINOP (add, octave_complex_matrix, octave_complex_matrix, add); + INSTALL_BINOP (sub, octave_complex_matrix, octave_complex_matrix, sub); + INSTALL_BINOP (mul, octave_complex_matrix, octave_complex_matrix, mul); + INSTALL_BINOP (div, octave_complex_matrix, octave_complex_matrix, div); + INSTALL_BINOP (pow, octave_complex_matrix, octave_complex_matrix, pow); + INSTALL_BINOP (ldiv, octave_complex_matrix, octave_complex_matrix, ldiv); + INSTALL_BINOP (lt, octave_complex_matrix, octave_complex_matrix, lt); + INSTALL_BINOP (le, octave_complex_matrix, octave_complex_matrix, le); + INSTALL_BINOP (eq, octave_complex_matrix, octave_complex_matrix, eq); + INSTALL_BINOP (ge, octave_complex_matrix, octave_complex_matrix, ge); + INSTALL_BINOP (gt, octave_complex_matrix, octave_complex_matrix, gt); + INSTALL_BINOP (ne, octave_complex_matrix, octave_complex_matrix, ne); + INSTALL_BINOP (el_mul, octave_complex_matrix, octave_complex_matrix, el_mul); + INSTALL_BINOP (el_div, octave_complex_matrix, octave_complex_matrix, el_div); + INSTALL_BINOP (el_pow, octave_complex_matrix, octave_complex_matrix, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex_matrix, octave_complex_matrix, el_ldiv); + INSTALL_BINOP (el_and, octave_complex_matrix, octave_complex_matrix, el_and); + INSTALL_BINOP (el_or, octave_complex_matrix, octave_complex_matrix, el_or); + + INSTALL_ASSIGNOP (asn_eq, octave_complex_matrix, octave_complex_matrix, assign); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cm-cs.cc @@ -0,0 +1,129 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-cx-mat.h" +#include "ov-complex.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex matrix by complex scalar ops. + +DEFBINOP_OP (add, complex_matrix, complex, +) +DEFBINOP_OP (sub, complex_matrix, complex, -) +DEFBINOP_OP (mul, complex_matrix, complex, *) + +DEFBINOP (div, complex_matrix, complex) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_complex&); + + Complex d = v2.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_matrix_value () / d); +} + +DEFBINOP_FN (pow, complex_matrix, complex, xpow) + +DEFBINOP (ldiv, complex_matrix, complex) +{ + BINOP_NONCONFORMANT ("operator \\"); +} + +DEFBINOP_FN (lt, complex_matrix, complex, mx_el_lt) +DEFBINOP_FN (le, complex_matrix, complex, mx_el_le) +DEFBINOP_FN (eq, complex_matrix, complex, mx_el_eq) +DEFBINOP_FN (ge, complex_matrix, complex, mx_el_ge) +DEFBINOP_FN (gt, complex_matrix, complex, mx_el_gt) +DEFBINOP_FN (ne, complex_matrix, complex, mx_el_ne) + +DEFBINOP_OP (el_mul, complex_matrix, complex, *) + +DEFBINOP (el_div, complex_matrix, complex) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_complex&); + + Complex d = v2.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_matrix_value () / d); +} + +DEFBINOP_FN (el_pow, complex_matrix, complex, elem_xpow) + +DEFBINOP (el_ldiv, complex_matrix, complex) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_complex&); + + return x_el_div (v2.complex_value (), v1.complex_matrix_value ()); +} + +DEFBINOP_FN (el_and, complex_matrix, complex, mx_el_and) +DEFBINOP_FN (el_or, complex_matrix, complex, mx_el_or) + +DEFASSIGNOP_FN (assign, complex_matrix, complex, assign) + +void +install_cm_cs_ops (void) +{ + INSTALL_BINOP (add, octave_complex_matrix, octave_complex, add); + INSTALL_BINOP (sub, octave_complex_matrix, octave_complex, sub); + INSTALL_BINOP (mul, octave_complex_matrix, octave_complex, mul); + INSTALL_BINOP (div, octave_complex_matrix, octave_complex, div); + INSTALL_BINOP (pow, octave_complex_matrix, octave_complex, pow); + INSTALL_BINOP (ldiv, octave_complex_matrix, octave_complex, ldiv); + INSTALL_BINOP (lt, octave_complex_matrix, octave_complex, lt); + INSTALL_BINOP (le, octave_complex_matrix, octave_complex, le); + INSTALL_BINOP (eq, octave_complex_matrix, octave_complex, eq); + INSTALL_BINOP (ge, octave_complex_matrix, octave_complex, ge); + INSTALL_BINOP (gt, octave_complex_matrix, octave_complex, gt); + INSTALL_BINOP (ne, octave_complex_matrix, octave_complex, ne); + INSTALL_BINOP (el_mul, octave_complex_matrix, octave_complex, el_mul); + INSTALL_BINOP (el_div, octave_complex_matrix, octave_complex, el_div); + INSTALL_BINOP (el_pow, octave_complex_matrix, octave_complex, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex_matrix, octave_complex, el_ldiv); + INSTALL_BINOP (el_and, octave_complex_matrix, octave_complex, el_and); + INSTALL_BINOP (el_or, octave_complex_matrix, octave_complex, el_or); + + INSTALL_ASSIGNOP (asn_eq, octave_complex_matrix, octave_complex, assign); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cm-m.cc @@ -0,0 +1,116 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "mx-cm-m.h" +#include "mx-m-cm.h" + +#include "gripes.h" +#include "ov.h" +#include "ov-cx-mat.h" +#include "ov-re-mat.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex matrix by matrix ops. + +DEFBINOP_OP (add, complex_matrix, matrix, +) +DEFBINOP_OP (sub, complex_matrix, matrix, -) +DEFBINOP_OP (mul, complex_matrix, matrix, *) + +DEFBINOP (div, complex_matrix, matrix) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_matrix&); + + return xdiv (v1.complex_matrix_value (), v2.matrix_value ()); +} + +DEFBINOPX (pow, complex_matrix, matrix) +{ + error ("can't do A ^ B for A and B both matrices"); + return octave_value (); +} + +DEFBINOP_FN (ldiv, complex_matrix, matrix, xleftdiv) + +DEFBINOP_FN (lt, complex_matrix, matrix, mx_el_lt) +DEFBINOP_FN (le, complex_matrix, matrix, mx_el_le) +DEFBINOP_FN (eq, complex_matrix, matrix, mx_el_eq) +DEFBINOP_FN (ge, complex_matrix, matrix, mx_el_ge) +DEFBINOP_FN (gt, complex_matrix, matrix, mx_el_gt) +DEFBINOP_FN (ne, complex_matrix, matrix, mx_el_ne) + +DEFBINOP_FN (el_mul, complex_matrix, matrix, product) +DEFBINOP_FN (el_div, complex_matrix, matrix, quotient) +DEFBINOP_FN (el_pow, complex_matrix, matrix, elem_xpow) + +DEFBINOP (el_ldiv, complex_matrix, matrix) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_matrix&); + + return quotient (v2.matrix_value (), v1.complex_matrix_value ()); +} + +DEFBINOP_FN (el_and, complex_matrix, matrix, mx_el_and) +DEFBINOP_FN (el_or, complex_matrix, matrix, mx_el_or) + +DEFASSIGNOP_FN (assign, complex_matrix, matrix, assign) + +void +install_cm_m_ops (void) +{ + INSTALL_BINOP (add, octave_complex_matrix, octave_matrix, add); + INSTALL_BINOP (sub, octave_complex_matrix, octave_matrix, sub); + INSTALL_BINOP (mul, octave_complex_matrix, octave_matrix, mul); + INSTALL_BINOP (div, octave_complex_matrix, octave_matrix, div); + INSTALL_BINOP (pow, octave_complex_matrix, octave_matrix, pow); + INSTALL_BINOP (ldiv, octave_complex_matrix, octave_matrix, ldiv); + INSTALL_BINOP (lt, octave_complex_matrix, octave_matrix, lt); + INSTALL_BINOP (le, octave_complex_matrix, octave_matrix, le); + INSTALL_BINOP (eq, octave_complex_matrix, octave_matrix, eq); + INSTALL_BINOP (ge, octave_complex_matrix, octave_matrix, ge); + INSTALL_BINOP (gt, octave_complex_matrix, octave_matrix, gt); + INSTALL_BINOP (ne, octave_complex_matrix, octave_matrix, ne); + INSTALL_BINOP (el_mul, octave_complex_matrix, octave_matrix, el_mul); + INSTALL_BINOP (el_div, octave_complex_matrix, octave_matrix, el_div); + INSTALL_BINOP (el_pow, octave_complex_matrix, octave_matrix, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex_matrix, octave_matrix, el_ldiv); + INSTALL_BINOP (el_and, octave_complex_matrix, octave_matrix, el_and); + INSTALL_BINOP (el_or, octave_complex_matrix, octave_matrix, el_or); + + INSTALL_ASSIGNOP (asn_eq, octave_complex_matrix, octave_matrix, assign); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cm-s.cc @@ -0,0 +1,131 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "mx-cm-s.h" + +#include "gripes.h" +#include "ov.h" +#include "ov-cx-mat.h" +#include "ov-scalar.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex matrix by scalar ops. + +DEFBINOP_OP (add, complex_matrix, scalar, +) +DEFBINOP_OP (sub, complex_matrix, scalar, -) +DEFBINOP_OP (mul, complex_matrix, scalar, *) + +DEFBINOP (div, complex_matrix, scalar) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_scalar&); + + double d = v2.double_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_matrix_value () / d); +} + +DEFBINOP_FN (pow, complex_matrix, scalar, xpow) + +DEFBINOP (ldiv, complex_matrix, scalar) +{ + BINOP_NONCONFORMANT ("operator \\"); +} + +DEFBINOP_FN (lt, complex_matrix, scalar, mx_el_lt) +DEFBINOP_FN (le, complex_matrix, scalar, mx_el_le) +DEFBINOP_FN (eq, complex_matrix, scalar, mx_el_eq) +DEFBINOP_FN (ge, complex_matrix, scalar, mx_el_ge) +DEFBINOP_FN (gt, complex_matrix, scalar, mx_el_gt) +DEFBINOP_FN (ne, complex_matrix, scalar, mx_el_ne) + +DEFBINOP_OP (el_mul, complex_matrix, scalar, *) + +DEFBINOP (el_div, complex_matrix, scalar) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_scalar&); + + double d = v2.double_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_matrix_value () / d); +} + +DEFBINOP_FN (el_pow, complex_matrix, scalar, elem_xpow) + +DEFBINOP (el_ldiv, complex_matrix, scalar) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_scalar&); + + return x_el_div (v2.double_value (), v1.complex_matrix_value ()); +} + +DEFBINOP_FN (el_and, complex_matrix, scalar, mx_el_and) +DEFBINOP_FN (el_or, complex_matrix, scalar, mx_el_or) + +DEFASSIGNOP_FN (assign, complex_matrix, scalar, assign) + +void +install_cm_s_ops (void) +{ + INSTALL_BINOP (add, octave_complex_matrix, octave_scalar, add); + INSTALL_BINOP (sub, octave_complex_matrix, octave_scalar, sub); + INSTALL_BINOP (mul, octave_complex_matrix, octave_scalar, mul); + INSTALL_BINOP (div, octave_complex_matrix, octave_scalar, div); + INSTALL_BINOP (pow, octave_complex_matrix, octave_scalar, pow); + INSTALL_BINOP (ldiv, octave_complex_matrix, octave_scalar, ldiv); + INSTALL_BINOP (lt, octave_complex_matrix, octave_scalar, lt); + INSTALL_BINOP (le, octave_complex_matrix, octave_scalar, le); + INSTALL_BINOP (eq, octave_complex_matrix, octave_scalar, eq); + INSTALL_BINOP (ge, octave_complex_matrix, octave_scalar, ge); + INSTALL_BINOP (gt, octave_complex_matrix, octave_scalar, gt); + INSTALL_BINOP (ne, octave_complex_matrix, octave_scalar, ne); + INSTALL_BINOP (el_mul, octave_complex_matrix, octave_scalar, el_mul); + INSTALL_BINOP (el_div, octave_complex_matrix, octave_scalar, el_div); + INSTALL_BINOP (el_pow, octave_complex_matrix, octave_scalar, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex_matrix, octave_scalar, el_ldiv); + INSTALL_BINOP (el_and, octave_complex_matrix, octave_scalar, el_and); + INSTALL_BINOP (el_or, octave_complex_matrix, octave_scalar, el_or); + + INSTALL_ASSIGNOP (asn_eq, octave_complex_matrix, octave_scalar, assign); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cs-cm.cc @@ -0,0 +1,129 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-complex.h" +#include "ov-cx-mat.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex scalar by complex matrix ops. + +DEFBINOP_OP (add, complex, complex_matrix, +) +DEFBINOP_OP (sub, complex, complex_matrix, -) +DEFBINOP_OP (mul, complex, complex_matrix, *) + +DEFBINOP (div, complex, complex_matrix) +{ + BINOP_NONCONFORMANT ("operator /"); +} + +DEFBINOP_FN (pow, complex, complex_matrix, xpow) + +DEFBINOP (ldiv, complex, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex_matrix&); + + Complex d = v1.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.complex_matrix_value () / d); +} + +DEFBINOP_FN (lt, complex, complex_matrix, mx_el_lt) +DEFBINOP_FN (le, complex, complex_matrix, mx_el_le) +DEFBINOP_FN (eq, complex, complex_matrix, mx_el_eq) +DEFBINOP_FN (ge, complex, complex_matrix, mx_el_ge) +DEFBINOP_FN (gt, complex, complex_matrix, mx_el_gt) +DEFBINOP_FN (ne, complex, complex_matrix, mx_el_ne) + +DEFBINOP_OP (el_mul, complex, complex_matrix, *) +DEFBINOP_FN (el_div, complex, complex_matrix, x_el_div) +DEFBINOP_FN (el_pow, complex, complex_matrix, elem_xpow) + +DEFBINOP (el_ldiv, complex, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex_matrix&); + + Complex d = v1.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.complex_matrix_value () / d); +} + +DEFBINOP_FN (el_and, complex, complex_matrix, mx_el_and) +DEFBINOP_FN (el_or, complex, complex_matrix, mx_el_or) + +DEFCONV (complex_matrix_conv, complex, complex_matrix) +{ + CAST_CONV_ARG (const octave_complex&); + + return new octave_complex_matrix (v.complex_matrix_value ()); +} + +void +install_cs_cm_ops (void) +{ + INSTALL_BINOP (add, octave_complex, octave_complex_matrix, add); + INSTALL_BINOP (sub, octave_complex, octave_complex_matrix, sub); + INSTALL_BINOP (mul, octave_complex, octave_complex_matrix, mul); + INSTALL_BINOP (div, octave_complex, octave_complex_matrix, div); + INSTALL_BINOP (pow, octave_complex, octave_complex_matrix, pow); + INSTALL_BINOP (ldiv, octave_complex, octave_complex_matrix, ldiv); + INSTALL_BINOP (lt, octave_complex, octave_complex_matrix, lt); + INSTALL_BINOP (le, octave_complex, octave_complex_matrix, le); + INSTALL_BINOP (eq, octave_complex, octave_complex_matrix, eq); + INSTALL_BINOP (ge, octave_complex, octave_complex_matrix, ge); + INSTALL_BINOP (gt, octave_complex, octave_complex_matrix, gt); + INSTALL_BINOP (ne, octave_complex, octave_complex_matrix, ne); + INSTALL_BINOP (el_mul, octave_complex, octave_complex_matrix, el_mul); + INSTALL_BINOP (el_div, octave_complex, octave_complex_matrix, el_div); + INSTALL_BINOP (el_pow, octave_complex, octave_complex_matrix, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex, octave_complex_matrix, el_ldiv); + INSTALL_BINOP (el_and, octave_complex, octave_complex_matrix, el_and); + INSTALL_BINOP (el_or, octave_complex, octave_complex_matrix, el_or); + + INSTALL_ASSIGNCONV (octave_complex, octave_complex_matrix, octave_complex_matrix); + + INSTALL_WIDENOP (octave_complex, octave_complex_matrix, complex_matrix_conv); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cs-cs.cc @@ -0,0 +1,194 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-complex.h" +#include "ov-cx-mat.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex scalar by complex scalar ops. + +DEFBINOP_OP (add, complex, complex, +) +DEFBINOP_OP (sub, complex, complex, -) +DEFBINOP_OP (mul, complex, complex, *) + +DEFBINOP (div, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + Complex d = v2.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_value () / d); +} + +DEFBINOP_FN (pow, complex, complex, xpow) + +DEFBINOP (ldiv, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + Complex d = v1.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.complex_value () / d); +} + +DEFBINOP (lt, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return real (v1.complex_value ()) < real (v2.complex_value ()); +} + +DEFBINOP (le, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return real (v1.complex_value ()) <= real (v2.complex_value ()); +} + +DEFBINOP (eq, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return v1.complex_value () == v2.complex_value (); +} + +DEFBINOP (ge, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return real (v1.complex_value ()) >= real (v2.complex_value ()); +} + +DEFBINOP (gt, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return real (v1.complex_value ()) > real (v2.complex_value ()); +} + +DEFBINOP (ne, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return v1.complex_value () != v2.complex_value (); +} + +DEFBINOP_OP (el_mul, complex, complex, *) + +DEFBINOP (el_div, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + Complex d = v2.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_value () / d); +} + +DEFBINOP_FN (el_pow, complex, complex, xpow) + +DEFBINOP (el_ldiv, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + Complex d = v1.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.complex_value () / d); +} + +DEFBINOP (el_and, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return v1.complex_value () != 0.0 && v2.complex_value () != 0.0; +} + +DEFBINOP (el_or, complex, complex) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_complex&); + + return v1.complex_value () != 0.0 || v2.complex_value () != 0.0; +} + +DEFCONV (complex_matrix_conv, complex, complex_matrix) +{ + CAST_CONV_ARG (const octave_complex&); + + return new octave_complex_matrix (v.complex_matrix_value ()); +} + +void +install_cs_cs_ops (void) +{ + INSTALL_BINOP (add, octave_complex, octave_complex, add); + INSTALL_BINOP (sub, octave_complex, octave_complex, sub); + INSTALL_BINOP (mul, octave_complex, octave_complex, mul); + INSTALL_BINOP (div, octave_complex, octave_complex, div); + INSTALL_BINOP (pow, octave_complex, octave_complex, pow); + INSTALL_BINOP (ldiv, octave_complex, octave_complex, ldiv); + INSTALL_BINOP (lt, octave_complex, octave_complex, lt); + INSTALL_BINOP (le, octave_complex, octave_complex, le); + INSTALL_BINOP (eq, octave_complex, octave_complex, eq); + INSTALL_BINOP (ge, octave_complex, octave_complex, ge); + INSTALL_BINOP (gt, octave_complex, octave_complex, gt); + INSTALL_BINOP (ne, octave_complex, octave_complex, ne); + INSTALL_BINOP (el_mul, octave_complex, octave_complex, el_mul); + INSTALL_BINOP (el_div, octave_complex, octave_complex, el_div); + INSTALL_BINOP (el_pow, octave_complex, octave_complex, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex, octave_complex, el_ldiv); + INSTALL_BINOP (el_and, octave_complex, octave_complex, el_and); + INSTALL_BINOP (el_or, octave_complex, octave_complex, el_or); + + INSTALL_ASSIGNCONV (octave_complex, octave_complex, octave_complex_matrix); + + INSTALL_WIDENOP (octave_complex, octave_complex_matrix, complex_matrix_conv); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cs-m.cc @@ -0,0 +1,133 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "mx-cs-m.h" +#include "mx-m-cs.h" + +#include "gripes.h" +#include "ov.h" +#include "ov-complex.h" +#include "ov-cx-mat.h" +#include "ov-re-mat.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex scalar by matrix ops. + +DEFBINOP_OP (add, complex, matrix, +) +DEFBINOP_OP (sub, complex, matrix, -) +DEFBINOP_OP (mul, complex, matrix, *) + +DEFBINOP (div, complex, matrix) +{ + BINOP_NONCONFORMANT ("operator /"); +} + +DEFBINOP_FN (pow, complex, matrix, xpow) + +DEFBINOP (ldiv, complex, matrix) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_matrix&); + + Complex d = v1.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.matrix_value () / d); +} + +DEFBINOP_FN (lt, complex, matrix, mx_el_lt) +DEFBINOP_FN (le, complex, matrix, mx_el_le) +DEFBINOP_FN (eq, complex, matrix, mx_el_eq) +DEFBINOP_FN (ge, complex, matrix, mx_el_ge) +DEFBINOP_FN (gt, complex, matrix, mx_el_gt) +DEFBINOP_FN (ne, complex, matrix, mx_el_ne) + +DEFBINOP_OP (el_mul, complex, matrix, *) +DEFBINOP_FN (el_div, complex, matrix, x_el_div) +DEFBINOP_FN (el_pow, complex, matrix, elem_xpow) + +DEFBINOP (el_ldiv, complex, matrix) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_matrix&); + + Complex d = v1.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.matrix_value () / d); +} + +DEFBINOP_FN (el_and, complex, matrix, mx_el_and) +DEFBINOP_FN (el_or, complex, matrix, mx_el_or) + +DEFCONV (complex_matrix_conv, complex, complex_matrix) +{ + CAST_CONV_ARG (const octave_complex&); + + return new octave_complex_matrix (v.complex_matrix_value ()); +} + +void +install_cs_m_ops (void) +{ + INSTALL_BINOP (add, octave_complex, octave_matrix, add); + INSTALL_BINOP (sub, octave_complex, octave_matrix, sub); + INSTALL_BINOP (mul, octave_complex, octave_matrix, mul); + INSTALL_BINOP (div, octave_complex, octave_matrix, div); + INSTALL_BINOP (pow, octave_complex, octave_matrix, pow); + INSTALL_BINOP (ldiv, octave_complex, octave_matrix, ldiv); + INSTALL_BINOP (lt, octave_complex, octave_matrix, lt); + INSTALL_BINOP (le, octave_complex, octave_matrix, le); + INSTALL_BINOP (eq, octave_complex, octave_matrix, eq); + INSTALL_BINOP (ge, octave_complex, octave_matrix, ge); + INSTALL_BINOP (gt, octave_complex, octave_matrix, gt); + INSTALL_BINOP (ne, octave_complex, octave_matrix, ne); + INSTALL_BINOP (el_mul, octave_complex, octave_matrix, el_mul); + INSTALL_BINOP (el_div, octave_complex, octave_matrix, el_div); + INSTALL_BINOP (el_pow, octave_complex, octave_matrix, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex, octave_matrix, el_ldiv); + INSTALL_BINOP (el_and, octave_complex, octave_matrix, el_and); + INSTALL_BINOP (el_or, octave_complex, octave_matrix, el_or); + + INSTALL_ASSIGNCONV (octave_complex, octave_matrix, octave_complex_matrix); + + INSTALL_WIDENOP (octave_complex, octave_complex_matrix, complex_matrix_conv); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-cs-s.cc @@ -0,0 +1,195 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "ov.h" +#include "ov-complex.h" +#include "ov-cx-mat.h" +#include "ov-scalar.h" +#include "ov-typeinfo.h" +#include "ops.h" +#include "xdiv.h" +#include "xpow.h" + +// complex scalar by scalar ops. + +DEFBINOP_OP (add, complex, scalar, +) +DEFBINOP_OP (sub, complex, scalar, -) +DEFBINOP_OP (mul, complex, scalar, *) + +DEFBINOP (div, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + double d = v2.double_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_value () / d); +} + +DEFBINOP_FN (pow, complex, scalar, xpow) + +DEFBINOP (ldiv, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + double d = v1.double_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.complex_value () / d); +} + +DEFBINOP (lt, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return real (v1.complex_value ()) < v2.double_value (); +} + +DEFBINOP (le, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return real (v1.complex_value ()) <= v2.double_value (); +} + +DEFBINOP (eq, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return v1.complex_value () == v2.double_value (); +} + +DEFBINOP (ge, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return real (v1.complex_value ()) >= v2.double_value (); +} + +DEFBINOP (gt, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return real (v1.complex_value ()) > v2.double_value (); +} + +DEFBINOP (ne, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return v1.complex_value () != v2.double_value (); +} + +DEFBINOP_OP (el_mul, complex, scalar, *) + +DEFBINOP (el_div, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + double d = v2.double_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.complex_value () / d); +} + +DEFBINOP_FN (el_pow, complex, scalar, xpow) + +DEFBINOP (el_ldiv, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + double d = v1.double_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v2.complex_value () / d); +} + +DEFBINOP (el_and, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return v1.complex_value () != 0.0 && v2.double_value (); +} + +DEFBINOP (el_or, complex, scalar) +{ + CAST_BINOP_ARGS (const octave_complex&, const octave_scalar&); + + return v1.complex_value () != 0.0 || v2.double_value (); +} + +DEFCONV (complex_matrix_conv, complex, complex_matrix) +{ + CAST_CONV_ARG (const octave_complex&); + + return new octave_complex_matrix (v.complex_matrix_value ()); +} + +void +install_cs_s_ops (void) +{ + INSTALL_BINOP (add, octave_complex, octave_scalar, add); + INSTALL_BINOP (sub, octave_complex, octave_scalar, sub); + INSTALL_BINOP (mul, octave_complex, octave_scalar, mul); + INSTALL_BINOP (div, octave_complex, octave_scalar, div); + INSTALL_BINOP (pow, octave_complex, octave_scalar, pow); + INSTALL_BINOP (ldiv, octave_complex, octave_scalar, ldiv); + INSTALL_BINOP (lt, octave_complex, octave_scalar, lt); + INSTALL_BINOP (le, octave_complex, octave_scalar, le); + INSTALL_BINOP (eq, octave_complex, octave_scalar, eq); + INSTALL_BINOP (ge, octave_complex, octave_scalar, ge); + INSTALL_BINOP (gt, octave_complex, octave_scalar, gt); + INSTALL_BINOP (ne, octave_complex, octave_scalar, ne); + INSTALL_BINOP (el_mul, octave_complex, octave_scalar, el_mul); + INSTALL_BINOP (el_div, octave_complex, octave_scalar, el_div); + INSTALL_BINOP (el_pow, octave_complex, octave_scalar, el_pow); + INSTALL_BINOP (el_ldiv, octave_complex, octave_scalar, el_ldiv); + INSTALL_BINOP (el_and, octave_complex, octave_scalar, el_and); + INSTALL_BINOP (el_or, octave_complex, octave_scalar, el_or); + + INSTALL_ASSIGNCONV (octave_complex, octave_scalar, octave_complex_matrix); + + INSTALL_WIDENOP (octave_complex, octave_complex_matrix, complex_matrix_conv); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-fil-b.cc @@ -0,0 +1,78 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <iostream.h> + +#include "mach-info.h" + +#include "error.h" +#include "oct-stream.h" +#include "ops.h" +#include "ov.h" +#include "ov-file.h" +#include "ov-bool.h" +#include "ov-typeinfo.h" + +// file by bool ops. + +DEFBINOP (lshift, file, bool) +{ + CAST_BINOP_ARGS (const octave_file&, const octave_bool&); + + octave_stream *oct_stream = v1.stream_value (); + + if (oct_stream) + { + ostream *osp = oct_stream->output_stream (); + + if (osp) + { + ostream& os = *osp; + + v2.print_raw (os); + } + else + error ("invalid file specified for binary operator `<<'"); + } + + return octave_value (oct_stream, v1.stream_number ()); +} + +void +install_fil_b_ops (void) +{ + INSTALL_BINOP (lshift, octave_file, octave_bool, lshift); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-fil-bm.cc @@ -0,0 +1,78 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <iostream.h> + +#include "mach-info.h" + +#include "error.h" +#include "oct-stream.h" +#include "ops.h" +#include "ov.h" +#include "ov-file.h" +#include "ov-bool-mat.h" +#include "ov-typeinfo.h" + +// file by bool matrix ops. + +DEFBINOP (lshift, file, bool_matrix) +{ + CAST_BINOP_ARGS (const octave_file&, const octave_bool_matrix&); + + octave_stream *oct_stream = v1.stream_value (); + + if (oct_stream) + { + ostream *osp = oct_stream->output_stream (); + + if (osp) + { + ostream& os = *osp; + + v2.print_raw (os); + } + else + error ("invalid file specified for binary operator `<<'"); + } + + return octave_value (oct_stream, v1.stream_number ()); +} + +void +install_fil_bm_ops (void) +{ + INSTALL_BINOP (lshift, octave_file, octave_bool_matrix, lshift); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-fil-cm.cc @@ -0,0 +1,78 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <iostream.h> + +#include "mach-info.h" + +#include "error.h" +#include "oct-stream.h" +#include "ops.h" +#include "ov.h" +#include "ov-file.h" +#include "ov-cx-mat.h" +#include "ov-typeinfo.h" + +// file by complex matrix ops. + +DEFBINOP (lshift, file, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_file&, const octave_complex_matrix&); + + octave_stream *oct_stream = v1.stream_value (); + + if (oct_stream) + { + ostream *osp = oct_stream->output_stream (); + + if (osp) + { + ostream& os = *osp; + + v2.print_raw (os); + } + else + error ("invalid file specified for binary operator `<<'"); + } + + return octave_value (oct_stream, v1.stream_number ()); +} + +void +install_fil_cm_ops (void) +{ + INSTALL_BINOP (lshift, octave_file, octave_complex_matrix, lshift); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-fil-cs.cc @@ -0,0 +1,78 @@ +/* + +Copyright (C) 1996, 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <iostream.h>