Mercurial > hg > octave-nkf
changeset 515:e078f05f4aac
[project @ 1994-07-13 02:31:31 by jwe]
Initial revision
author | jwe |
---|---|
date | Wed, 13 Jul 1994 02:31:31 +0000 |
parents | 20d2061944ee |
children | 309fc59f66ee |
files | src/f-eval.cc src/f-fill.cc src/f-input.cc src/find.cc src/log.cc src/minmax.cc src/sort.cc |
diffstat | 7 files changed, 1670 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/src/f-eval.cc @@ -0,0 +1,124 @@ +// f-eval.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "unwind-prot.h" +#include "tree-const.h" +#include "variables.h" +#include "octave.h" +#include "input.h" +#include "parse.h" +#include "lex.h" +#include "f-eval.h" + +Octave_object +feval (const Octave_object& args, int nargout) +{ +// Assumes that we have been given the correct number of arguments. + + Octave_object retval; + + tree_fvc *fcn = is_valid_function (args(1), "feval", 1); + if (fcn != (tree_fvc *) NULL) + { + int nargin = args.length () - 1; + Octave_object tmp_args (nargin); + for (int i = 0; i < nargin; i++) + tmp_args(i) = args(i+1); + retval = fcn->eval (0, nargout, tmp_args); + } + + return retval; +} + +tree_constant +eval_string (const char *string, int print, int ans_assign, + int& parse_status) +{ + begin_unwind_frame ("eval_string"); + + unwind_protect_int (get_input_from_eval_string); + unwind_protect_ptr (global_command); + unwind_protect_ptr (current_eval_string); + + get_input_from_eval_string = 1; + current_eval_string = string; + + YY_BUFFER_STATE old_buf = current_buffer (); + YY_BUFFER_STATE new_buf = create_buffer ((FILE *) NULL); + + add_unwind_protect (restore_input_buffer, (void *) old_buf); + add_unwind_protect (delete_input_buffer, (void *) new_buf); + + switch_to_buffer (new_buf); + + unwind_protect_ptr (curr_sym_tab); + + reset_parser (); + + parse_status = yyparse (); + +// Important to reset the idea of where input is coming from before +// trying to eval the command we just parsed -- it might contain the +// name of an function file that still needs to be parsed! + + tree *command = global_command; + + run_unwind_frame ("eval_string"); + + tree_constant retval; + + if (parse_status == 0 && command != NULL_TREE) + { + retval = command->eval (print); + delete command; + } + + return retval; +} + +tree_constant +eval_string (const tree_constant& arg, int& parse_status) +{ + if (! arg.is_string_type ()) + { + error ("eval: expecting string argument"); + return -1; + } + + char *string = arg.string_value (); + +// Yes Virginia, we always print here... + + return eval_string (string, 1, 1, parse_status); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/f-fill.cc @@ -0,0 +1,175 @@ +// f-fill.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "tree-const.h" +#include "user-prefs.h" +#include "utils.h" +#include "error.h" +#include "f-fill.h" + +#ifndef MIN +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#endif + +static void +check_dimensions (int& nr, int& nc, const char *warnfor) +{ + if (nr < 0 || nc < 0) + { + if (user_pref.treat_neg_dim_as_zero) + nr = nc = 0; + else + error ("%s: can't create a matrix with negative dimensions", + warnfor); + } +} + +static void +get_dimensions (const tree_constant& a, const char *warn_for, + int& nr, int& nc) +{ + tree_constant tmpa = a.make_numeric (); + + if (tmpa.is_scalar_type ()) + { + double tmp = tmpa.double_value (); + nr = nc = NINT (tmp); + } + else + { + nr = tmpa.rows (); + nc = tmpa.columns (); + + if ((nr == 1 && nc == 2) || (nr == 2 && nc == 1)) + { + ColumnVector v = tmpa.to_vector (); + + nr = NINT (v.elem (0)); + nc = NINT (v.elem (1)); + } + else + warning ("%s (A): use %s (size (A)) instead", warn_for, warn_for); + } + + check_dimensions (nr, nc, warn_for); // May set error_state. +} + +static void +get_dimensions (const tree_constant& a, const tree_constant& b, + const char *warn_for, int& nr, int& nc) +{ + tree_constant tmpa = a.make_numeric (); + tree_constant tmpb = b.make_numeric (); + + if (tmpa.is_scalar_type () && tmpb.is_scalar_type ()) + { + nr = NINT (tmpa.double_value ()); + nc = NINT (tmpb.double_value ()); + + check_dimensions (nr, nc, warn_for); // May set error_state. + } + else + error ("%s: expecting two scalar arguments", warn_for); +} + +tree_constant +fill_matrix (const tree_constant& a, double val, const char *warn_for) +{ + int nr, nc; + get_dimensions (a, warn_for, nr, nc); + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, val); + + return m; +} + +tree_constant +fill_matrix (const tree_constant& a, const tree_constant& b, + double val, const char *warn_for) +{ + int nr, nc; + get_dimensions (a, b, warn_for, nr, nc); // May set error_state. + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, val); + + return m; +} + +tree_constant +identity_matrix (const tree_constant& a) +{ + int nr, nc; + get_dimensions (a, "eye", nr, nc); // May set error_state. + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, 0.0); + + if (nr > 0 && nc > 0) + { + int n = MIN (nr, nc); + for (int i = 0; i < n; i++) + m.elem (i, i) = 1.0; + } + + return m; +} + +tree_constant +identity_matrix (const tree_constant& a, const tree_constant& b) +{ + int nr, nc; + get_dimensions (a, b, "eye", nr, nc); // May set error_state. + + if (error_state) + return tree_constant (); + + Matrix m (nr, nc, 0.0); + + if (nr > 0 && nc > 0) + { + int n = MIN (nr, nc); + for (int i = 0; i < n; i++) + m.elem (i, i) = 1.0; + } + + return m; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/f-input.cc @@ -0,0 +1,135 @@ +// f-input.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include <string.h> + +#include "tree-const.h" +#include "octave-hist.h" +#include "pager.h" +#include "input.h" +#include "f-eval.h" +#include "f-input.h" + +static int +match_sans_spaces (const char *standard, const char *test) +{ + const char *tp = test; + while (*tp == ' ' || *tp == '\t') + tp++; + + const char *ep = test + strlen (test) - 1; + while (*ep == ' ' || *ep == '\t') + ep--; + + int len = ep - tp + 1; + + return (strncmp (standard, tp, len) == 0); +} + +tree_constant +get_user_input (const Octave_object& args, int nargout, int debug = 0) +{ + tree_constant retval; + + int nargin = args.length (); + + int read_as_string = 0; + + if (nargin == 3) + { + if (args(2).is_string_type () + && strcmp ("s", args(2).string_value ()) == 0) + read_as_string++; + else + { + error ("input: unrecognized second argument"); + return retval; + } + } + + char *prompt = "debug> "; + if (nargin > 1) + { + if (args(1).is_string_type ()) + prompt = args(1).string_value (); + else + { + error ("input: unrecognized argument"); + return retval; + } + } + + again: + + flush_output_to_pager (); + + char *input_buf = gnu_readline (prompt); + + if (input_buf != (char *) NULL) + { + if (input_buf) + maybe_save_history (input_buf); + + int len = strlen (input_buf); + + if (len < 1) + { + if (debug) + goto again; + else + return retval; + } + + if (match_sans_spaces ("exit", input_buf) + || match_sans_spaces ("quit", input_buf) + || match_sans_spaces ("return", input_buf)) + return tree_constant (); + else if (read_as_string) + retval = input_buf; + else + { + int parse_status = 0; + retval = eval_string (input_buf, 0, 0, parse_status); + if (debug && retval.is_defined ()) + retval.eval (1); + } + } + else + error ("input: reading user-input failed!"); + + if (debug) + goto again; + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/find.cc @@ -0,0 +1,212 @@ +// f-find.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "tree-const.h" +#include "error.h" +#include "f-find.h" + +static Octave_object +find_to_fortran_idx (const ColumnVector i_idx, const ColumnVector j_idx, + const tree_constant& val, int nr, int nc, int nargout) +{ + Octave_object retval (nargout); + + switch (nargout) + { + 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); + retval(0) = tree_constant (tmp, 1); +// If you want this to work more like Matlab, use the following line +// instead of the previous one. +// retval(0) = tree_constant (tmp, (nr != 1)); + } + break; + case 3: + retval(2) = val; + case 2: + retval(0) = tree_constant (i_idx, 1); +// If you want this to work more like Matlab, use the following line +// instead of the previous one. +// retval(0) = tree_constant (i_idx, (nr != 1)); + retval(1) = tree_constant (j_idx, 1); + break; + default: + panic_impossible (); + break; + } + + return retval; +} + +static Octave_object +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.elem (i, j) != 0.0) + count++; + + Matrix result; + Octave_object retval (nargout, result); + + 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.elem (i, j); + if (d != 0.0) + { + i_idx (count) = i + 1; + j_idx (count) = j + 1; + v (count) = d; + count++; + } + } + + tree_constant tmp (v, 1); + return find_to_fortran_idx (i_idx, j_idx, tmp, m_nr, m_nc, nargout); +} + +static Octave_object +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.elem (i, j) != 0.0) + count++; + + Matrix result; + Octave_object retval (nargout, result); + + 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.elem (i, j); + if (c != 0.0) + { + i_idx (count) = i; + j_idx (count) = j; + v (count) = c; + count++; + } + } + + tree_constant tmp (v, 1); + return find_to_fortran_idx (i_idx, j_idx, tmp, m_nr, m_nc, nargout); +} + +Octave_object +find_nonzero_elem_idx (const tree_constant& a, int nargout) +{ + Matrix result; + + nargout = (nargout == 0) ? 1 : nargout; + Octave_object retval (nargout, result); + + tree_constant tmp = a.make_numeric (); + + switch (tmp.const_type ()) + { + case tree_constant_rep::matrix_constant: + { + Matrix m = tmp.matrix_value (); + return find_nonzero_elem_idx (m, nargout); + } + break; + case tree_constant_rep::scalar_constant: + { + double d = tmp.double_value (); + if (d != 0.0) + { + retval(0) = 1.0; + if (nargout > 1) + retval(1) = 1.0; + if (nargout > 2) + retval(2) = d; + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = tmp.complex_matrix_value (); + return find_nonzero_elem_idx (m, nargout); + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex c = tmp.complex_value (); + if (c != 0.0) + { + retval(0) = 1.0; + if (nargout > 1) + retval(1) = 1.0; + if (nargout > 2) + retval(2) = c; + } + } + break; + default: + break; + } + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/log.cc @@ -0,0 +1,264 @@ +// f-log.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "EIG.h" + +#include "tree-const.h" +#include "user-prefs.h" +#include "error.h" +#include "gripes.h" +#include "f-log.h" + +// XXX FIXME XXX -- the next two functions (and expm) should really be just +// one... + +Octave_object +matrix_log (const tree_constant& a) +{ + Octave_object retval (1); + + tree_constant tmp = a.make_numeric ();; + + if (tmp.rows () == 0 || tmp.columns () == 0) + { + int flag = user_pref.propagate_empty_matrices; + if (flag != 0) + { + if (flag < 0) + gripe_empty_arg ("logm", 0); + Matrix m; + retval(0) = m; + return retval; + } + else + gripe_empty_arg ("logm", 1); + } + + switch (tmp.const_type ()) + { + case tree_constant_rep::matrix_constant: + { + Matrix m = tmp.matrix_value (); + + 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.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = log (real (elt)); + else + lambda.elem (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = tmp.complex_matrix_value (); + + 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.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = log (real (elt)); + else + lambda.elem (i) = log (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + break; + case tree_constant_rep::scalar_constant: + { + double d = tmp.double_value (); + if (d > 0.0) + retval(0) = log (d); + else + { + Complex dtmp (d); + retval(0) = log (dtmp); + } + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex c = tmp.complex_value (); + retval(0) = log (c); + } + break; + default: + break; + } + return retval; +} + +Octave_object +matrix_sqrt (const tree_constant& a) +{ + Octave_object retval (1); + + tree_constant tmp = a.make_numeric ();; + + if (tmp.rows () == 0 || tmp.columns () == 0) + { + int flag = user_pref.propagate_empty_matrices; + if (flag != 0) + { + if (flag < 0) + gripe_empty_arg ("sqrtm", 0); + Matrix m; + retval(0) = m; + return retval; + } + else + gripe_empty_arg ("sqrtm", 1); + } + + switch (tmp.const_type ()) + { + case tree_constant_rep::matrix_constant: + { + Matrix m = tmp.matrix_value (); + + 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.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = sqrt (real (elt)); + else + lambda.elem (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = tmp.complex_matrix_value (); + + 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.elem (i); + if (imag (elt) == 0.0 && real (elt) > 0.0) + lambda.elem (i) = sqrt (real (elt)); + else + lambda.elem (i) = sqrt (elt); + } + + ComplexDiagMatrix D (lambda); + ComplexMatrix result = Q * D * Q.inverse (); + + retval(0) = result; + } + } + break; + case tree_constant_rep::scalar_constant: + { + double d = tmp.double_value (); + if (d > 0.0) + retval(0) = sqrt (d); + else + { + Complex dtmp (d); + retval(0) = sqrt (dtmp); + } + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex c = tmp.complex_value (); + retval(0) = log (c); + } + break; + default: + break; + } + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/minmax.cc @@ -0,0 +1,495 @@ +// f-minmax.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include <math.h> + +#include "tree-const.h" +#include "error.h" +#include "f-minmax.h" + +#ifndef MAX +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#endif + +#ifndef MIN +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#endif + +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.elem (i, j); + double b_elem = b.elem (i, j); + result.elem (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++) + for (int i = 0; i < nr; i++) + { + double abs_a_elem = abs (a.elem (i, j)); + double abs_b_elem = abs (b.elem (i, j)); + if (abs_a_elem < abs_b_elem) + result.elem (i, j) = a.elem (i, j); + else + result.elem (i, j) = b.elem (i, j); + } + + 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.elem (i, j); + double b_elem = b.elem (i, j); + result.elem (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++) + for (int i = 0; i < nr; i++) + { + double abs_a_elem = abs (a.elem (i, j)); + double abs_b_elem = abs (b.elem (i, j)); + if (abs_a_elem > abs_b_elem) + result.elem (i, j) = a.elem (i, j); + else + result.elem (i, j) = b.elem (i, j); + } + + return result; +} + +Octave_object +column_min (const Octave_object& args, int nargout) +{ + Octave_object retval; + + tree_constant arg1; + tree_constant arg2; + tree_constant_rep::constant_type arg1_type = + tree_constant_rep::unknown_constant; + tree_constant_rep::constant_type arg2_type = + tree_constant_rep::unknown_constant; + + int nargin = args.length (); + + switch (nargin) + { + case 3: + arg2 = args(2).make_numeric (); + arg2_type = arg2.const_type (); +// Fall through... + case 2: + arg1 = args(1).make_numeric (); + arg1_type = arg1.const_type (); + break; + default: + panic_impossible (); + break; + } + + if (nargin == 2 && (nargout == 1 || nargout == 0)) + { + retval.resize (1); + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + retval(0) = arg1.double_value (); + break; + case tree_constant_rep::complex_scalar_constant: + retval(0) = arg1.complex_value (); + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + retval(0) = m.row_min (); + else + retval(0) = tree_constant (m.column_min (), 0); + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + retval(0) = m.row_min (); + else + retval(0) = tree_constant (m.column_min (), 0); + } + break; + default: + panic_impossible (); + break; + } + } + else if (nargin == 2 && nargout == 2) + { + retval.resize (2); + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + { + retval(0) = arg1.double_value (); + retval(1) = 1; + } + break; + case tree_constant_rep::complex_scalar_constant: + { + retval(0) = arg1.complex_value (); + retval(1) = 1; + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + { + retval(0) = m.row_min (); + retval(1) = m.row_min_loc (); + } + else + { + retval(0) = tree_constant (m.column_min (), 0); + retval(1) = tree_constant (m.column_min_loc (), 0); + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + { + retval(0) = m.row_min (); + retval(1) = m.row_min_loc (); + } + else + { + retval(0) = tree_constant (m.column_min (), 0); + retval(1) = tree_constant (m.column_min_loc (), 0); + } + } + break; + default: + panic_impossible (); + break; + } + } + else if (nargin == 3) + { + if (arg1.rows () == arg2.rows () + && arg1.columns () == arg2.columns ()) + { + retval.resize (1); + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + { + double result; + double a_elem = arg1.double_value (); + double b_elem = arg2.double_value (); + result = MIN (a_elem, b_elem); + retval(0) = result; + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex result; + Complex a_elem = arg1.complex_value (); + Complex b_elem = arg2.complex_value (); + if (abs (a_elem) < abs (b_elem)) + result = a_elem; + else + result = b_elem; + retval(0) = result; + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix result; + result = min (arg1.matrix_value (), arg2.matrix_value ()); + retval(0) = result; + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix result; + result = min (arg1.complex_matrix_value (), + arg2.complex_matrix_value ()); + retval(0) = result; + } + break; + default: + panic_impossible (); + break; + } + } + else + error ("min: nonconformant matrices"); + } + else + panic_impossible (); + + return retval; +} + +Octave_object +column_max (const Octave_object& args, int nargout) +{ + Octave_object retval; + + tree_constant arg1; + tree_constant arg2; + tree_constant_rep::constant_type arg1_type = + tree_constant_rep::unknown_constant; + tree_constant_rep::constant_type arg2_type = + tree_constant_rep::unknown_constant; + + int nargin = args.length (); + + switch (nargin) + { + case 3: + arg2 = args(2).make_numeric (); + arg2_type = arg2.const_type (); +// Fall through... + case 2: + arg1 = args(1).make_numeric (); + arg1_type = arg1.const_type (); + break; + default: + panic_impossible (); + break; + } + + if (nargin == 2 && (nargout == 1 || nargout == 0)) + { + retval.resize (1); + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + retval(0) = arg1.double_value (); + break; + case tree_constant_rep::complex_scalar_constant: + retval(0) = arg1.complex_value (); + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + retval(0) = m.row_max (); + else + retval(0) = tree_constant (m.column_max (), 0); + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + retval(0) = m.row_max (); + else + retval(0) = tree_constant (m.column_max (), 0); + } + break; + default: + panic_impossible (); + break; + } + } + else if (nargin == 2 && nargout == 2) + { + retval.resize (2); + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + { + retval(0) = arg1.double_value (); + retval(1) = 1; + } + break; + case tree_constant_rep::complex_scalar_constant: + { + retval(0) = arg1.complex_value (); + retval(1) = 1; + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix m = arg1.matrix_value (); + if (m.rows () == 1) + { + retval(0) = m.row_max (); + retval(1) = m.row_max_loc (); + } + else + { + retval(0) = tree_constant (m.column_max (), 0); + retval(1) = tree_constant (m.column_max_loc (), 0); + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix m = arg1.complex_matrix_value (); + if (m.rows () == 1) + { + retval(0) = m.row_max (); + retval(1) = m.row_max_loc (); + } + else + { + retval(0) = tree_constant (m.column_max (), 0); + retval(1) = tree_constant (m.column_max_loc (), 0); + } + } + break; + default: + panic_impossible (); + break; + } + } + else if (nargin == 3) + { + if (arg1.rows () == arg2.rows () + && arg1.columns () == arg2.columns ()) + { + retval.resize (1); + switch (arg1_type) + { + case tree_constant_rep::scalar_constant: + { + double result; + double a_elem = arg1.double_value (); + double b_elem = arg2.double_value (); + result = MAX (a_elem, b_elem); + retval(0) = result; + } + break; + case tree_constant_rep::complex_scalar_constant: + { + Complex result; + Complex a_elem = arg1.complex_value (); + Complex b_elem = arg2.complex_value (); + if (abs (a_elem) > abs (b_elem)) + result = a_elem; + else + result = b_elem; + retval(0) = result; + } + break; + case tree_constant_rep::matrix_constant: + { + Matrix result; + result = max (arg1.matrix_value (), arg2.matrix_value ()); + retval(0) = result; + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix result; + result = max (arg1.complex_matrix_value (), + arg2.complex_matrix_value ()); + retval(0) = result; + } + break; + default: + panic_impossible (); + break; + } + } + else + error ("max: nonconformant matrices"); + } + else + panic_impossible (); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/src/sort.cc @@ -0,0 +1,265 @@ +// f-sort.cc -*- C++ -*- +/* + +Copyright (C) 1994 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, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "tree-const.h" +#include "f-sort.h" + +static void +mx_sort (Matrix& m, Matrix& idx, int return_idx) +{ + int nr = m.rows (); + int nc = m.columns (); + idx.resize (nr, nc); + int i, j; + + if (return_idx) + { + for (j = 0; j < nc; j++) + for (i = 0; i < nr; i++) + idx.elem (i, j) = i+1; + } + + for (j = 0; j < nc; j++) + { + for (int gap = nr/2; gap > 0; gap /= 2) + for (i = gap; i < nr; i++) + for (int k = i - gap; + k >= 0 && m.elem (k, j) > m.elem (k+gap, j); + k -= gap) + { + double tmp = m.elem (k, j); + m.elem (k, j) = m.elem (k+gap, j); + m.elem (k+gap, j) = tmp; + + if (return_idx) + { + double tmp = idx.elem (k, j); + idx.elem (k, j) = idx.elem (k+gap, j); + idx.elem (k+gap, j) = tmp; + } + } + } +} + +static void +mx_sort (RowVector& v, RowVector& idx, int return_idx) +{ + int n = v.capacity (); + idx.resize (n); + int i; + + if (return_idx) + for (i = 0; i < n; i++) + idx.elem (i) = i+1; + + for (int gap = n/2; gap > 0; gap /= 2) + for (i = gap; i < n; i++) + for (int k = i - gap; + k >= 0 && v.elem (k) > v.elem (k+gap); + k -= gap) + { + double tmp = v.elem (k); + v.elem (k) = v.elem (k+gap); + v.elem (k+gap) = tmp; + + if (return_idx) + { + double tmp = idx.elem (k); + idx.elem (k) = idx.elem (k+gap); + idx.elem (k+gap) = tmp; + } + } +} + +static void +mx_sort (ComplexMatrix& cm, Matrix& idx, int return_idx) +{ + int nr = cm.rows (); + int nc = cm.columns (); + idx.resize (nr, nc); + int i, j; + + if (return_idx) + { + for (j = 0; j < nc; j++) + for (i = 0; i < nr; i++) + idx.elem (i, j) = i+1; + } + + for (j = 0; j < nc; j++) + { + for (int gap = nr/2; gap > 0; gap /= 2) + for (i = gap; i < nr; i++) + for (int k = i - gap; + k >= 0 && abs (cm.elem (k, j)) > abs (cm.elem (k+gap, j)); + k -= gap) + { + Complex ctmp = cm.elem (k, j); + cm.elem (k, j) = cm.elem (k+gap, j); + cm.elem (k+gap, j) = ctmp; + + if (return_idx) + { + double tmp = idx.elem (k, j); + idx.elem (k, j) = idx.elem (k+gap, j); + idx.elem (k+gap, j) = tmp; + } + } + } +} + +static void +mx_sort (ComplexRowVector& cv, RowVector& idx, int return_idx) +{ + int n = cv.capacity (); + idx.resize (n); + int i; + + if (return_idx) + for (i = 0; i < n; i++) + idx.elem (i) = i+1; + + for (int gap = n/2; gap > 0; gap /= 2) + for (i = gap; i < n; i++) + for (int k = i - gap; + k >= 0 && abs (cv.elem (k)) > abs (cv.elem (k+gap)); + k -= gap) + { + Complex tmp = cv.elem (k); + cv.elem (k) = cv.elem (k+gap); + cv.elem (k+gap) = tmp; + + if (return_idx) + { + double tmp = idx.elem (k); + idx.elem (k) = idx.elem (k+gap); + idx.elem (k+gap) = tmp; + } + } +} + +Octave_object +sort (const Octave_object& args, int nargout) +{ +// Assumes that we have been given the correct number of arguments. + + Octave_object retval; + + int return_idx = nargout > 1; + if (return_idx) + retval.resize (2); + else + retval.resize (1); + + switch (args(1).const_type ()) + { + case tree_constant_rep::scalar_constant: + { + retval(0) = args(1).double_value (); + if (return_idx) + retval(1) = 1.0; + } + break; + case tree_constant_rep::complex_scalar_constant: + { + retval(0) = args(1).complex_value (); + if (return_idx) + retval(1) = 1.0; + } + break; + case tree_constant_rep::string_constant: + case tree_constant_rep::range_constant: + case tree_constant_rep::matrix_constant: + { + Matrix m = args(1).to_matrix (); + if (m.rows () == 1) + { + int nc = m.columns (); + RowVector v (nc); + for (int i = 0; i < nc; i++) + v.elem (i) = m.elem (0, i); + RowVector idx; + mx_sort (v, idx, return_idx); + + retval(0) = tree_constant (v, 0); + if (return_idx) + retval(1) = tree_constant (idx, 0); + } + else + { +// Sorts m in place, optionally computes index Matrix. + Matrix idx; + mx_sort (m, idx, return_idx); + + retval(0) = m; + if (return_idx) + retval(1) = idx; + } + } + break; + case tree_constant_rep::complex_matrix_constant: + { + ComplexMatrix cm = args(1).complex_matrix_value (); + if (cm.rows () == 1) + { + int nc = cm.columns (); + ComplexRowVector cv (nc); + for (int i = 0; i < nc; i++) + cv.elem (i) = cm.elem (0, i); + RowVector idx; + mx_sort (cv, idx, return_idx); + + retval(0) = tree_constant (cv, 0); + if (return_idx) + retval(1) = tree_constant (idx, 0); + } + else + { +// Sorts cm in place, optionally computes index Matrix. + Matrix idx; + mx_sort (cm, idx, return_idx); + + retval(0) = cm; + if (return_idx) + retval(1) = idx; + } + } + break; + default: + panic_impossible (); + break; + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/