Mercurial > octave-nkf
view src/tc-rep.cc @ 636:fae2bd91c027
[project @ 1994-08-23 18:39:50 by jwe]
author | jwe |
---|---|
date | Tue, 23 Aug 1994 18:39:50 +0000 |
parents | 9aef0a53eee7 |
children | 36ba0576bd1b |
line wrap: on
line source
// tc-rep.cc -*- C++ -*- /* Copyright (C) 1992, 1993, 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 #if defined (__GNUG__) #pragma implementation #endif #include <ctype.h> #include <string.h> #include <fstream.h> #include <iostream.h> #include <strstream.h> #include "mx-base.h" #include "Range.h" #include "arith-ops.h" #include "variables.h" #include "sysdep.h" #include "error.h" #include "gripes.h" #include "user-prefs.h" #include "utils.h" #include "pager.h" #include "pr-output.h" #include "tree-const.h" #include "idx-vector.h" #include "tc-inlines.cc" // How about a few macros? #define TC_REP tree_constant::tree_constant_rep #ifndef MAX #define MAX(a,b) ((a) > (b) ? (a) : (b)) #endif #ifndef MIN #define MIN(a,b) ((a) < (b) ? (a) : (b)) #endif #ifndef ABS #define ABS(x) (((x) < 0) ? (-x) : (x)) #endif // The following are used by some of the functions in the // tree_constant_rep class that must deal with real and complex // matrices. This was not done with overloaded or virtual functions // from the Matrix class because there is no clean way to do that -- // the necessary functions (like elem) need to return values of // different types... // Given a tree_constant, and the names to be used for the real and // complex matrix and their dimensions, declare a real or complex // matrix, and initialize it from the tree_constant. Note that m, cm, // nr, and nc must not be previously declared, and they must not be // expressions. Since only one of the matrices will be defined after // this macro is used, only one set of dimesions is declared. // This macro only makes sense inside a friend or member function of // the tree_constant_rep class #define REP_RHS_MATRIX(tc,m,cm,nr,nc) \ int nr = 0; \ int nc = 0; \ Matrix m; \ ComplexMatrix cm; \ if ((tc).const_type () == TC_REP::complex_matrix_constant) \ { \ cm = (tc).complex_matrix_value (); \ nr = (cm).rows (); \ nc = (cm).columns (); \ } \ else if ((tc).const_type () == TC_REP::matrix_constant) \ { \ m = (tc).matrix_value (); \ nr = (m).rows (); \ nc = (m).columns (); \ } \ else \ abort (); // Assign a real or complex value to a tree_constant. // // This macro only makes sense inside a friend or member function of // the tree_constant_rep class. #define REP_ELEM_ASSIGN(i,j,rval,cval,real_type) \ do \ { \ if (type_tag == TC_REP::matrix_constant) \ { \ if (real_type) \ matrix->elem ((i), (j)) = (rval); \ else \ abort (); \ } \ else \ { \ if (real_type) \ complex_matrix->elem ((i), (j)) = (rval); \ else \ complex_matrix->elem ((i), (j)) = (cval); \ } \ } \ while (0) // Given a real and complex matrix and row and column dimensions, // declare both and size one of them. Only one of the matrices should // be used after this macro has been used. // This macro only makes sense inside a friend or member function of // the tree_constant_rep class. #define CRMATRIX(m,cm,nr,nc) \ Matrix m; \ ComplexMatrix cm; \ if (type_tag == TC_REP::matrix_constant) \ (m).resize ((nr), (nc)); \ else if (type_tag == complex_matrix_constant) \ (cm).resize ((nr), (nc)); \ else \ abort (); \ // Assign a real or complex matrix to a tree constant. // This macro only makes sense inside a friend or member function of // the tree_constant_rep class. #define ASSIGN_CRMATRIX_TO(tc,m,cm) \ do \ { \ if (type_tag == matrix_constant) \ tc = tree_constant (m); \ else \ tc = tree_constant (cm); \ } \ while (0) // Assign an element of this tree_constant_rep's real or complex // matrix to another real or complex matrix. // This macro only makes sense inside a friend or member function of // the tree_constant_rep class. #define CRMATRIX_ASSIGN_REP_ELEM(m,cm,i1,j1,i2,j2) \ do \ { \ if (type_tag == matrix_constant) \ (m).elem ((i1), (j1)) = matrix->elem ((i2), (j2)); \ else \ (cm).elem ((i1), (j1)) = complex_matrix->elem ((i2), (j2)); \ } \ while (0) // Assign a value to an element of a real or complex matrix. Assumes // that the lhs and rhs are either both real or both complex types. #define CRMATRIX_ASSIGN_ELEM(m,cm,i,j,rval,cval,real_type) \ do \ { \ if (real_type) \ (m).elem ((i), (j)) = (rval); \ else \ (cm).elem ((i), (j)) = (cval); \ } \ while (0) // A couple of handy helper functions. static int any_element_less_than (const Matrix& a, double val) { int nr = a.rows (); int nc = a.columns (); for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) if (a.elem (i, j) < val) return 1; return 0; } static int any_element_greater_than (const Matrix& a, double val) { int nr = a.rows (); int nc = a.columns (); for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) if (a.elem (i, j) > val) return 1; return 0; } static int any_element_is_complex (const ComplexMatrix& a) { int nr = a.rows (); int nc = a.columns (); for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) if (imag (a.elem (i, j)) != 0.0) return 1; return 0; } static int valid_scalar_indices (const Octave_object& args) { int nargin = args.length (); return ((nargin == 3 && args(2).valid_as_scalar_index () && args(1).valid_as_scalar_index ()) || (nargin == 2 && args(1).valid_as_scalar_index ())); } // Now, the classes. // The real representation of constants. TC_REP::tree_constant_rep (void) { type_tag = unknown_constant; orig_text = 0; } TC_REP::tree_constant_rep (double d) { scalar = d; type_tag = scalar_constant; orig_text = 0; } TC_REP::tree_constant_rep (const Matrix& m) { if (m.rows () == 1 && m.columns () == 1) { scalar = m.elem (0, 0); type_tag = scalar_constant; } else { matrix = new Matrix (m); type_tag = matrix_constant; } orig_text = 0; } TC_REP::tree_constant_rep (const DiagMatrix& d) { if (d.rows () == 1 && d.columns () == 1) { scalar = d.elem (0, 0); type_tag = scalar_constant; } else { matrix = new Matrix (d); type_tag = matrix_constant; } orig_text = 0; } TC_REP::tree_constant_rep (const RowVector& v, int prefer_column_vector) { int len = v.capacity (); if (len == 1) { scalar = v.elem (0); type_tag = scalar_constant; } else { int pcv = (prefer_column_vector < 0) ? user_pref.prefer_column_vectors : prefer_column_vector; if (pcv) { Matrix m (len, 1); for (int i = 0; i < len; i++) m.elem (i, 0) = v.elem (i); matrix = new Matrix (m); type_tag = matrix_constant; } else { Matrix m (1, len); for (int i = 0; i < len; i++) m.elem (0, i) = v.elem (i); matrix = new Matrix (m); type_tag = matrix_constant; } } orig_text = 0; } TC_REP::tree_constant_rep (const ColumnVector& v, int prefer_column_vector) { int len = v.capacity (); if (len == 1) { scalar = v.elem (0); type_tag = scalar_constant; } else { int pcv = (prefer_column_vector < 0) ? user_pref.prefer_column_vectors : prefer_column_vector; if (pcv) { Matrix m (len, 1); for (int i = 0; i < len; i++) m.elem (i, 0) = v.elem (i); matrix = new Matrix (m); type_tag = matrix_constant; } else { Matrix m (1, len); for (int i = 0; i < len; i++) m.elem (0, i) = v.elem (i); matrix = new Matrix (m); type_tag = matrix_constant; } } orig_text = 0; } TC_REP::tree_constant_rep (const Complex& c) { complex_scalar = new Complex (c); type_tag = complex_scalar_constant; orig_text = 0; } TC_REP::tree_constant_rep (const ComplexMatrix& m) { if (m.rows () == 1 && m.columns () == 1) { complex_scalar = new Complex (m.elem (0, 0)); type_tag = complex_scalar_constant; } else { complex_matrix = new ComplexMatrix (m); type_tag = complex_matrix_constant; } orig_text = 0; } TC_REP::tree_constant_rep (const ComplexDiagMatrix& d) { if (d.rows () == 1 && d.columns () == 1) { complex_scalar = new Complex (d.elem (0, 0)); type_tag = complex_scalar_constant; } else { complex_matrix = new ComplexMatrix (d); type_tag = complex_matrix_constant; } orig_text = 0; } TC_REP::tree_constant_rep (const ComplexRowVector& v, int prefer_column_vector) { int len = v.capacity (); if (len == 1) { complex_scalar = new Complex (v.elem (0)); type_tag = complex_scalar_constant; } else { int pcv = (prefer_column_vector < 0) ? user_pref.prefer_column_vectors : prefer_column_vector; if (pcv) { ComplexMatrix m (len, 1); for (int i = 0; i < len; i++) m.elem (i, 0) = v.elem (i); complex_matrix = new ComplexMatrix (m); type_tag = complex_matrix_constant; } else { ComplexMatrix m (1, len); for (int i = 0; i < len; i++) m.elem (0, i) = v.elem (i); complex_matrix = new ComplexMatrix (m); type_tag = complex_matrix_constant; } } orig_text = 0; } TC_REP::tree_constant_rep (const ComplexColumnVector& v, int prefer_column_vector) { int len = v.capacity (); if (len == 1) { complex_scalar = new Complex (v.elem (0)); type_tag = complex_scalar_constant; } else { int pcv = (prefer_column_vector < 0) ? user_pref.prefer_column_vectors : prefer_column_vector; if (pcv) { ComplexMatrix m (len, 1); for (int i = 0; i < len; i++) m.elem (i, 0) = v.elem (i); complex_matrix = new ComplexMatrix (m); type_tag = complex_matrix_constant; } else { ComplexMatrix m (1, len); for (int i = 0; i < len; i++) m.elem (0, i) = v.elem (i); complex_matrix = new ComplexMatrix (m); type_tag = complex_matrix_constant; } } orig_text = 0; } TC_REP::tree_constant_rep (const char *s) { string = strsave (s); type_tag = string_constant; orig_text = 0; } TC_REP::tree_constant_rep (double b, double l, double i) { range = new Range (b, l, i); int nel = range->nelem (); if (nel < 0) { delete range; type_tag = unknown_constant; if (nel == -1) ::error ("number of elements in range exceeds INT_MAX"); else ::error ("invalid range"); } else if (nel > 1) type_tag = range_constant; else { delete range; if (nel == 1) { scalar = b; type_tag = scalar_constant; } else if (nel == 0) { matrix = new Matrix (); type_tag = matrix_constant; } else panic_impossible (); } orig_text = 0; } TC_REP::tree_constant_rep (const Range& r) { if (r.nelem () > 1) { range = new Range (r); type_tag = range_constant; } else if (r.nelem () == 1) { scalar = r.base (); type_tag = scalar_constant; } else if (r.nelem () == 0) { matrix = new Matrix (); type_tag = matrix_constant; } else panic_impossible (); orig_text = 0; } TC_REP::tree_constant_rep (TC_REP::constant_type t) { assert (t == magic_colon); type_tag = magic_colon; orig_text = 0; } TC_REP::tree_constant_rep (const tree_constant_rep& t) { type_tag = t.type_tag; switch (t.type_tag) { case unknown_constant: break; case scalar_constant: scalar = t.scalar; break; case matrix_constant: matrix = new Matrix (*(t.matrix)); break; case string_constant: string = strsave (t.string); break; case complex_matrix_constant: complex_matrix = new ComplexMatrix (*(t.complex_matrix)); break; case complex_scalar_constant: complex_scalar = new Complex (*(t.complex_scalar)); break; case range_constant: range = new Range (*(t.range)); break; case magic_colon: break; default: panic_impossible (); break; } orig_text = strsave (t.orig_text); } TC_REP::~tree_constant_rep (void) { switch (type_tag) { case unknown_constant: case scalar_constant: case magic_colon: break; case matrix_constant: delete matrix; break; case complex_scalar_constant: delete complex_scalar; break; case complex_matrix_constant: delete complex_matrix; break; case string_constant: delete [] string; break; case range_constant: delete range; break; default: panic_impossible (); break; } delete [] orig_text; } #if defined (MDEBUG) void * TC_REP::operator new (size_t size) { tree_constant_rep *p = ::new tree_constant_rep; cerr << "TC_REP::new(): " << p << "\n"; return p; } void TC_REP::operator delete (void *p, size_t size) { cerr << "TC_REP::delete(): " << p << "\n"; ::delete p; } #endif int TC_REP::rows (void) const { int retval = -1; switch (type_tag) { case scalar_constant: case complex_scalar_constant: retval = 1; break; case string_constant: case range_constant: retval = (columns () > 0); break; case matrix_constant: retval = matrix->rows (); break; case complex_matrix_constant: retval = complex_matrix->rows (); break; case magic_colon: ::error ("invalid use of colon operator"); break; case unknown_constant: retval = 0; break; default: panic_impossible (); break; } return retval; } int TC_REP::columns (void) const { int retval = -1; switch (type_tag) { case scalar_constant: case complex_scalar_constant: retval = 1; break; case matrix_constant: retval = matrix->columns (); break; case complex_matrix_constant: retval = complex_matrix->columns (); break; case string_constant: retval = strlen (string); break; case range_constant: retval = range->nelem (); break; case magic_colon: ::error ("invalid use of colon operator"); break; case unknown_constant: retval = 0; break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::all (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.all (); } tree_constant retval; switch (type_tag) { case scalar_constant: { double status = (scalar != 0.0); retval = tree_constant (status); } break; case matrix_constant: { Matrix m = matrix->all (); retval = tree_constant (m); } break; case complex_scalar_constant: { double status = (*complex_scalar != 0.0); retval = tree_constant (status); } break; case complex_matrix_constant: { Matrix m = complex_matrix->all (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } tree_constant TC_REP::any (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.any (); } tree_constant retval; switch (type_tag) { case scalar_constant: { double status = (scalar != 0.0); retval = tree_constant (status); } break; case matrix_constant: { Matrix m = matrix->any (); retval = tree_constant (m); } break; case complex_scalar_constant: { double status = (*complex_scalar != 0.0); retval = tree_constant (status); } break; case complex_matrix_constant: { Matrix m = complex_matrix->any (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } int TC_REP::valid_as_scalar_index (void) const { return (type_tag == magic_colon || (type_tag == scalar_constant && NINT (scalar) == 1) || (type_tag == range_constant && range->nelem () == 1 && NINT (range->base ()) == 1)); } int TC_REP::is_true (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.is_true (); } int retval; switch (type_tag) { case scalar_constant: retval = (scalar != 0.0); break; case matrix_constant: { Matrix m = (matrix->all ()) . all (); retval = (m.rows () == 1 && m.columns () == 1 && m.elem (0, 0) != 0.0); } break; case complex_scalar_constant: retval = (*complex_scalar != 0.0); break; case complex_matrix_constant: { Matrix m = (complex_matrix->all ()) . all (); retval = (m.rows () == 1 && m.columns () == 1 && m.elem (0, 0) != 0.0); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } static void warn_implicit_conversion (const char *from, const char *to) { warning ("implicit conversion from %s to %s", from, to); } double TC_REP::double_value (int force_string_conversion) const { double retval = octave_NaN; switch (type_tag) { case scalar_constant: retval = scalar; break; case matrix_constant: { if (user_pref.do_fortran_indexing && rows () > 0 && columns () > 0) retval = matrix->elem (0, 0); else gripe_invalid_conversion ("real matrix", "real scalar"); } break; case complex_matrix_constant: case complex_scalar_constant: { int flag = user_pref.ok_to_lose_imaginary_part; if (flag < 0) warn_implicit_conversion ("complex scalar", "real scalar"); if (flag) { if (type_tag == complex_scalar_constant) retval = ::real (*complex_scalar); else if (type_tag == complex_matrix_constant) { if (user_pref.do_fortran_indexing && rows () > 0 && columns () > 0) retval = ::real (complex_matrix->elem (0, 0)); else gripe_invalid_conversion ("complex matrix", "real scalar"); } else panic_impossible (); } else gripe_invalid_conversion ("complex scalar", "real scalar"); } break; case string_constant: { int flag = force_string_conversion; if (! flag) flag = user_pref.implicit_str_to_num_ok; if (flag < 0) warn_implicit_conversion ("string", "real scalar"); int len = strlen (string); if (flag && (len == 1 || (len > 1 && user_pref.do_fortran_indexing))) retval = toascii ((int) string[0]); else gripe_invalid_conversion ("string", "real scalar"); } break; case range_constant: { int nel = range->nelem (); if (nel == 1 || (nel > 1 && user_pref.do_fortran_indexing)) retval = range->base (); else gripe_invalid_conversion ("range", "real scalar"); } break; default: gripe_invalid_conversion (type_as_string (), "real scalar"); break; } return retval; } Matrix TC_REP::matrix_value (int force_string_conversion) const { Matrix retval; switch (type_tag) { case scalar_constant: retval = Matrix (1, 1, scalar); break; case matrix_constant: retval = *matrix; break; case complex_scalar_constant: case complex_matrix_constant: { int flag = user_pref.ok_to_lose_imaginary_part; if (flag < 0) warn_implicit_conversion ("complex matrix", "real matrix"); if (flag) { if (type_tag == complex_scalar_constant) retval = Matrix (1, 1, ::real (*complex_scalar)); else if (type_tag == complex_matrix_constant) retval = ::real (*complex_matrix); else panic_impossible (); } else gripe_invalid_conversion ("complex matrix", "real matrix"); } break; case string_constant: { int flag = force_string_conversion; if (! flag) flag = user_pref.implicit_str_to_num_ok; if (flag < 0) warn_implicit_conversion ("string", "real matrix"); if (flag) { int len = strlen (string); retval.resize (1, len); if (len > 1) { for (int i = 0; i < len; i++) retval.elem (0, i) = toascii ((int) string[i]); } else if (len == 1) retval.elem (0, 0) = toascii ((int) string[0]); else panic_impossible (); } else gripe_invalid_conversion ("string", "real matrix"); } break; case range_constant: retval = range->matrix_value (); break; default: gripe_invalid_conversion (type_as_string (), "real matrix"); break; } return retval; } Complex TC_REP::complex_value (int force_string_conversion) const { Complex retval (octave_NaN, octave_NaN); switch (type_tag) { case complex_scalar_constant: retval = *complex_scalar; break; case scalar_constant: retval = scalar; break; case complex_matrix_constant: case matrix_constant: { if (user_pref.do_fortran_indexing && rows () > 0 && columns () > 0) { if (type_tag == complex_matrix_constant) retval = complex_matrix->elem (0, 0); else retval = matrix->elem (0, 0); } else gripe_invalid_conversion ("real matrix", "real scalar"); } break; case string_constant: { int flag = force_string_conversion; if (! flag) flag = user_pref.implicit_str_to_num_ok; if (flag < 0) warn_implicit_conversion ("string", "complex scalar"); int len = strlen (string); if (flag && (len == 1 || (len > 1 && user_pref.do_fortran_indexing))) retval = toascii ((int) string[0]); else gripe_invalid_conversion ("string", "complex scalar"); } break; case range_constant: { int nel = range->nelem (); if (nel == 1 || (nel > 1 && user_pref.do_fortran_indexing)) retval = range->base (); else gripe_invalid_conversion ("range", "complex scalar"); } break; default: gripe_invalid_conversion (type_as_string (), "complex scalar"); break; } return retval; } ComplexMatrix TC_REP::complex_matrix_value (int force_string_conversion) const { ComplexMatrix retval; switch (type_tag) { case scalar_constant: retval = ComplexMatrix (1, 1, Complex (scalar)); break; case complex_scalar_constant: retval = ComplexMatrix (1, 1, *complex_scalar); break; case matrix_constant: retval = ComplexMatrix (*matrix); break; case complex_matrix_constant: retval = *complex_matrix; break; case string_constant: { int flag = force_string_conversion; if (! flag) flag = user_pref.implicit_str_to_num_ok; if (flag < 0) warn_implicit_conversion ("string", "complex matrix"); if (flag) { int len = strlen (string); retval.resize (1, len); if (len > 1) { for (int i = 0; i < len; i++) retval.elem (0, i) = toascii ((int) string[i]); } else if (len == 1) retval.elem (0, 0) = toascii ((int) string[0]); else panic_impossible (); } else gripe_invalid_conversion ("string", "real matrix"); } break; case range_constant: retval = range->matrix_value (); break; default: gripe_invalid_conversion (type_as_string (), "complex matrix"); break; } return retval; } char * TC_REP::string_value (void) const { if (type_tag == string_constant) return string; else { gripe_invalid_conversion (type_as_string (), "string"); return 0; } } Range TC_REP::range_value (void) const { assert (type_tag == range_constant); return *range; } // This could be made more efficient by doing all the work here rather // than relying on matrix_value() to do any possible type conversions. ColumnVector TC_REP::vector_value (int force_string_conversion, int force_vector_conversion) const { ColumnVector retval; Matrix m = matrix_value (force_string_conversion); if (error_state) return retval; int nr = m.rows (); int nc = m.columns (); if (nr == 1) { retval.resize (nc); for (int i = 0; i < nc; i++) retval.elem (i) = m (0, i); } else if (nc == 1) { retval.resize (nr); for (int i = 0; i < nr; i++) retval.elem (i) = m.elem (i, 0); } else if (nr > 0 && nc > 0 && (user_pref.do_fortran_indexing || force_vector_conversion)) { retval.resize (nr * nc); int k = 0; for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) retval.elem (k++) = m.elem (i, j); } else gripe_invalid_conversion ("real matrix", "real vector"); return retval; } // This could be made more efficient by doing all the work here rather // than relying on complex_matrix_value() to do any possible type // conversions. ComplexColumnVector TC_REP::complex_vector_value (int force_string_conversion, int force_vector_conversion) const { ComplexColumnVector retval; ComplexMatrix m = complex_matrix_value (force_string_conversion); if (error_state) return retval; int nr = m.rows (); int nc = m.columns (); if (nr == 1) { retval.resize (nc); for (int i = 0; i < nc; i++) retval.elem (i) = m (0, i); } else if (nc == 1) { retval.resize (nr); for (int i = 0; i < nr; i++) retval.elem (i) = m.elem (i, 0); } else if (nr > 0 && nc > 0 && (user_pref.do_fortran_indexing || force_vector_conversion)) { retval.resize (nr * nc); int k = 0; for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) retval.elem (k++) = m.elem (i, j); } else gripe_invalid_conversion ("complex matrix", "complex vector"); return retval; } tree_constant TC_REP::convert_to_str (void) { tree_constant retval; switch (type_tag) { case complex_scalar_constant: case scalar_constant: { double d = double_value (); int i = NINT (d); // Warn about out of range conversions? char s[2]; s[0] = (char) i; s[1] = '\0'; retval = tree_constant (s); } break; case complex_matrix_constant: case matrix_constant: { ColumnVector v = vector_value (); int len = v.length (); if (len == 0) ::error ("can only convert vectors and scalars to strings"); else { char *s = new char [len+1]; s[len] = '\0'; for (int i = 0; i < len; i++) { double d = v.elem (i); int ival = NINT (d); // Warn about out of range conversions? s[i] = (char) ival; } retval = tree_constant (s); delete [] s; } } break; case range_constant: { Range r = range_value (); double b = r.base (); double incr = r.inc (); int nel = r.nelem (); char *s = new char [nel+1]; s[nel] = '\0'; for (int i = 0; i < nel; i++) { double d = b + i * incr; int ival = NINT (d); // Warn about out of range conversions? s[i] = (char) ival; } retval = tree_constant (s); delete [] s; } break; case string_constant: retval = string; break; case magic_colon: default: panic_impossible (); break; } return retval; } void TC_REP::convert_to_row_or_column_vector (void) { assert (type_tag == matrix_constant || type_tag == complex_matrix_constant); int nr = rows (); int nc = columns (); if (nr == 1 || nc == 1) return; int len = nr * nc; assert (len > 0); int new_nr = 1; int new_nc = 1; if (user_pref.prefer_column_vectors) new_nr = len; else new_nc = len; if (type_tag == matrix_constant) { Matrix *m = new Matrix (new_nr, new_nc); double *cop_out = matrix->fortran_vec (); for (int i = 0; i < len; i++) { if (new_nr == 1) m->elem (0, i) = *cop_out++; else m->elem (i, 0) = *cop_out++; } delete matrix; matrix = m; } else { ComplexMatrix *cm = new ComplexMatrix (new_nr, new_nc); Complex *cop_out = complex_matrix->fortran_vec (); for (int i = 0; i < len; i++) { if (new_nr == 1) cm->elem (0, i) = *cop_out++; else cm->elem (i, 0) = *cop_out++; } delete complex_matrix; complex_matrix = cm; } } void TC_REP::force_numeric (int force_str_conv) { switch (type_tag) { case scalar_constant: case matrix_constant: case complex_scalar_constant: case complex_matrix_constant: break; case string_constant: { if (! force_str_conv && ! user_pref.implicit_str_to_num_ok) { ::error ("failed to convert `%s' to a numeric type --", string); ::error ("default conversion turned off"); // Abort! jump_to_top_level (); } int len = strlen (string); if (len > 1) { type_tag = matrix_constant; Matrix *tm = new Matrix (1, len); for (int i = 0; i < len; i++) tm->elem (0, i) = toascii ((int) string[i]); matrix = tm; } else if (len == 1) { type_tag = scalar_constant; scalar = toascii ((int) string[0]); } else if (len == 0) { type_tag = matrix_constant; matrix = new Matrix (0, 0); } else panic_impossible (); } break; case range_constant: { int len = range->nelem (); if (len > 1) { type_tag = matrix_constant; Matrix *tm = new Matrix (1, len); double b = range->base (); double increment = range->inc (); for (int i = 0; i < len; i++) tm->elem (0, i) = b + i * increment; matrix = tm; } else if (len == 1) { type_tag = scalar_constant; scalar = range->base (); } } break; case magic_colon: default: panic_impossible (); break; } } tree_constant TC_REP::make_numeric (int force_str_conv) const { tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar); break; case matrix_constant: retval = tree_constant (*matrix); break; case complex_scalar_constant: retval = tree_constant (*complex_scalar); break; case complex_matrix_constant: retval = tree_constant (*complex_matrix); break; case string_constant: retval = tree_constant (string); retval.force_numeric (force_str_conv); break; case range_constant: retval = tree_constant (*range); retval.force_numeric (force_str_conv); break; case magic_colon: default: panic_impossible (); break; } return retval; } void TC_REP::bump_value (tree_expression::type etype) { switch (etype) { case tree_expression::increment: switch (type_tag) { case scalar_constant: scalar++; break; case matrix_constant: *matrix = *matrix + 1.0; break; case complex_scalar_constant: *complex_scalar = *complex_scalar + 1.0; break; case complex_matrix_constant: *complex_matrix = *complex_matrix + 1.0; break; case string_constant: ::error ("string++ and ++string not implemented yet, ok?"); break; case range_constant: range->set_base (range->base () + 1.0); range->set_limit (range->limit () + 1.0); break; case magic_colon: default: panic_impossible (); break; } break; case tree_expression::decrement: switch (type_tag) { case scalar_constant: scalar--; break; case matrix_constant: *matrix = *matrix - 1.0; break; case string_constant: ::error ("string-- and -- string not implemented yet, ok?"); break; case range_constant: range->set_base (range->base () - 1.0); range->set_limit (range->limit () - 1.0); break; case magic_colon: default: panic_impossible (); break; } break; default: panic_impossible (); break; } } void TC_REP::resize (int i, int j) { switch (type_tag) { case matrix_constant: matrix->resize (i, j); break; case complex_matrix_constant: complex_matrix->resize (i, j); break; default: panic_impossible (); break; } } void TC_REP::resize (int i, int j, double val) { switch (type_tag) { case matrix_constant: matrix->resize (i, j, val); break; case complex_matrix_constant: complex_matrix->resize (i, j, val); break; default: panic_impossible (); break; } } void TC_REP::maybe_resize (int i, int j) { int nr = rows (); int nc = columns (); i++; j++; assert (i > 0 && j > 0); if (i > nr || j > nc) { if (user_pref.resize_on_range_error) resize (MAX (i, nr), MAX (j, nc), 0.0); else { if (i > nr) ::error ("row index = %d exceeds max row dimension = %d", i, nr); if (j > nc) ::error ("column index = %d exceeds max column dimension = %d", j, nc); } } } void TC_REP::maybe_resize (int i, force_orient f_orient) { int nr = rows (); int nc = columns (); i++; assert (i >= 0 && (nr <= 1 || nc <= 1)); // This function never reduces the size of a vector, and all vectors // have dimensions of at least 0x0. If i is 0, it is either because // a vector has been indexed with a vector of all zeros (in which case // the index vector is empty and nothing will happen) or a vector has // been indexed with 0 (an error which will be caught elsewhere). if (i == 0) return; if (nr <= 1 && nc <= 1 && i >= 1) { if (user_pref.resize_on_range_error) { if (f_orient == row_orient) resize (1, i, 0.0); else if (f_orient == column_orient) resize (i, 1, 0.0); else if (user_pref.prefer_column_vectors) resize (i, 1, 0.0); else resize (1, i, 0.0); } else ::error ("matrix index = %d exceeds max dimension = %d", i, nc); } else if (nr == 1 && i > nc) { if (user_pref.resize_on_range_error) resize (1, i, 0.0); else ::error ("matrix index = %d exceeds max dimension = %d", i, nc); } else if (nc == 1 && i > nr) { if (user_pref.resize_on_range_error) resize (i, 1, 0.0); else ::error ("matrix index = %d exceeds max dimension = ", i, nc); } } void TC_REP::stash_original_text (char *s) { orig_text = strsave (s); } // Indexing functions. tree_constant TC_REP::do_index (const Octave_object& args) { tree_constant retval; if (error_state) return retval; if (rows () == 0 || columns () == 0) { ::error ("attempt to index empty matrix"); return retval; } switch (type_tag) { case complex_scalar_constant: case scalar_constant: retval = do_scalar_index (args); break; case complex_matrix_constant: case matrix_constant: retval = do_matrix_index (args); break; case string_constant: gripe_string_invalid (); // retval = do_string_index (args); break; case magic_colon: case range_constant: // This isn\'t great, but it\'s easier than implementing a lot of // range indexing functions. force_numeric (); assert (type_tag != magic_colon && type_tag != range_constant); retval = do_index (args); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_scalar_index (const Octave_object& args) const { tree_constant retval; if (valid_scalar_indices (args)) { if (type_tag == scalar_constant) retval = scalar; else if (type_tag == complex_scalar_constant) retval = *complex_scalar; else panic_impossible (); return retval; } else { int rows = 0; int cols = 0; int nargin = args.length (); switch (nargin) { case 3: { if (args(2).is_matrix_type ()) { Matrix mj = args(2).matrix_value (); idx_vector j (mj, user_pref.do_fortran_indexing, ""); if (! j) return retval; int len = j.length (); if (len == j.ones_count ()) cols = len; } else if (args(2).const_type () == magic_colon || (args(2).is_scalar_type () && NINT (args(2).double_value ()) == 1)) { cols = 1; } else break; } // Fall through... case 2: { if (args(1).is_matrix_type ()) { Matrix mi = args(1).matrix_value (); idx_vector i (mi, user_pref.do_fortran_indexing, ""); if (! i) return retval; int len = i.length (); if (len == i.ones_count ()) rows = len; } else if (args(1).const_type () == magic_colon || (args(1).is_scalar_type () && NINT (args(1).double_value ()) == 1)) { rows = 1; } else if (args(1).is_scalar_type () && NINT (args(1).double_value ()) == 0) { return Matrix (); } else break; if (cols == 0) { if (user_pref.prefer_column_vectors) cols = 1; else { cols = rows; rows = 1; } } if (type_tag == scalar_constant) { return Matrix (rows, cols, scalar); } else if (type_tag == complex_scalar_constant) { return ComplexMatrix (rows, cols, *complex_scalar); } else panic_impossible (); } break; default: ::error ("invalid number of arguments for scalar type"); return tree_constant (); break; } } ::error ("index invalid or out of range for scalar type"); return tree_constant (); } tree_constant TC_REP::do_matrix_index (const Octave_object& args) const { tree_constant retval; int nargin = args.length (); switch (nargin) { case 2: if (args.length () <= 0) ::error ("matrix index is null"); else if (args(1).is_undefined ()) ::error ("matrix index is a null expression"); else retval = do_matrix_index (args(1)); break; case 3: if (args.length () <= 0) ::error ("matrix indices are null"); else if (args(1).is_undefined ()) ::error ("first matrix index is a null expression"); else if (args(2).is_undefined ()) ::error ("second matrix index is a null expression"); else retval = do_matrix_index (args(1), args(2)); break; default: ::error ("too many indices for matrix expression"); break; } return retval; } tree_constant TC_REP::do_matrix_index (const tree_constant& i_arg) const { tree_constant retval; int nr = rows (); int nc = columns (); if (user_pref.do_fortran_indexing) retval = fortran_style_matrix_index (i_arg); else if (nr <= 1 || nc <= 1) retval = do_vector_index (i_arg); else ::error ("single index only valid for row or column vector"); return retval; } tree_constant TC_REP::do_matrix_index (const tree_constant& i_arg, const tree_constant& j_arg) const { tree_constant retval; tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type itype = tmp_i.const_type (); switch (itype) { case complex_scalar_constant: case scalar_constant: { int i = tree_to_mat_idx (tmp_i.double_value ()); if (index_check (i, "row") < 0) return tree_constant (); retval = do_matrix_index (i, j_arg); } break; case complex_matrix_constant: case matrix_constant: { Matrix mi = tmp_i.matrix_value (); idx_vector iv (mi, user_pref.do_fortran_indexing, "row", rows ()); if (! iv) return tree_constant (); if (iv.length () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else retval = do_matrix_index (iv, j_arg); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range ri = tmp_i.range_value (); int nr = rows (); if (nr == 2 && is_zero_one (ri)) { retval = do_matrix_index (1, j_arg); } else if (nr == 2 && is_one_zero (ri)) { retval = do_matrix_index (0, j_arg); } else { if (index_check (ri, "row") < 0) return tree_constant (); retval = do_matrix_index (ri, j_arg); } } break; case magic_colon: retval = do_matrix_index (magic_colon, j_arg); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_matrix_index (TC_REP::constant_type mci) const { assert (mci == magic_colon); tree_constant retval; int nr = rows (); int nc = columns (); int size = nr * nc; if (size > 0) { CRMATRIX (m, cm, size, 1); int idx = 0; for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) { CRMATRIX_ASSIGN_REP_ELEM (m, cm, idx, 0, i, j); idx++; } ASSIGN_CRMATRIX_TO (retval, m, cm); } return retval; } tree_constant TC_REP::fortran_style_matrix_index (const tree_constant& i_arg) const { tree_constant retval; tree_constant tmp_i = i_arg.make_numeric_or_magic (); TC_REP::constant_type itype = tmp_i.const_type (); int nr = rows (); int nc = columns (); switch (itype) { case complex_scalar_constant: case scalar_constant: { int i = NINT (tmp_i.double_value ()); int ii = fortran_row (i, nr) - 1; int jj = fortran_column (i, nr) - 1; if (index_check (i-1, "") < 0) return tree_constant (); if (range_max_check (i-1, nr * nc) < 0) return tree_constant (); retval = do_matrix_index (ii, jj); } break; case complex_matrix_constant: case matrix_constant: { Matrix mi = tmp_i.matrix_value (); if (mi.rows () == 0 || mi.columns () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else { // Yes, we really do want to call this with mi. retval = fortran_style_matrix_index (mi); } } break; case string_constant: gripe_string_invalid (); break; case range_constant: gripe_range_invalid (); break; case magic_colon: retval = do_matrix_index (magic_colon); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::fortran_style_matrix_index (const Matrix& mi) const { assert (is_matrix_type ()); tree_constant retval; int nr = rows (); int nc = columns (); int len = nr * nc; int index_nr = mi.rows (); int index_nc = mi.columns (); if (index_nr >= 1 && index_nc >= 1) { const double *cop_out = 0; const Complex *c_cop_out = 0; int real_type = type_tag == matrix_constant; if (real_type) cop_out = matrix->data (); else c_cop_out = complex_matrix->data (); const double *cop_out_index = mi.data (); idx_vector iv (mi, 1, "", len); if (! iv) return tree_constant (); int result_size = iv.length (); if (nc == 1 || (nr != 1 && iv.one_zero_only ())) { CRMATRIX (m, cm, result_size, 1); for (int i = 0; i < result_size; i++) { int idx = iv.elem (i); CRMATRIX_ASSIGN_ELEM (m, cm, i, 0, cop_out [idx], c_cop_out [idx], real_type); } ASSIGN_CRMATRIX_TO (retval, m, cm); } else if (nr == 1) { CRMATRIX (m, cm, 1, result_size); for (int i = 0; i < result_size; i++) { int idx = iv.elem (i); CRMATRIX_ASSIGN_ELEM (m, cm, 0, i, cop_out [idx], c_cop_out [idx], real_type); } ASSIGN_CRMATRIX_TO (retval, m, cm); } else { CRMATRIX (m, cm, index_nr, index_nc); for (int j = 0; j < index_nc; j++) for (int i = 0; i < index_nr; i++) { double tmp = *cop_out_index++; int idx = tree_to_mat_idx (tmp); CRMATRIX_ASSIGN_ELEM (m, cm, i, j, cop_out [idx], c_cop_out [idx], real_type); } ASSIGN_CRMATRIX_TO (retval, m, cm); } } else { if (index_nr == 0 || index_nc == 0) ::error ("empty matrix invalid as index"); else ::error ("invalid matrix index"); return tree_constant (); } return retval; } tree_constant TC_REP::do_vector_index (const tree_constant& i_arg) const { tree_constant retval; tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type itype = tmp_i.const_type (); int nr = rows (); int nc = columns (); int len = MAX (nr, nc); assert ((nr == 1 || nc == 1) && ! user_pref.do_fortran_indexing); int swap_indices = (nr == 1); switch (itype) { case complex_scalar_constant: case scalar_constant: { int i = tree_to_mat_idx (tmp_i.double_value ()); if (index_check (i, "") < 0) return tree_constant (); if (swap_indices) { if (range_max_check (i, nc) < 0) return tree_constant (); retval = do_matrix_index (0, i); } else { if (range_max_check (i, nr) < 0) return tree_constant (); retval = do_matrix_index (i, 0); } } break; case complex_matrix_constant: case matrix_constant: { Matrix mi = tmp_i.matrix_value (); if (mi.rows () == 0 || mi.columns () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else { idx_vector iv (mi, user_pref.do_fortran_indexing, "", len); if (! iv) return tree_constant (); if (swap_indices) { if (range_max_check (iv.max (), nc) < 0) return tree_constant (); retval = do_matrix_index (0, iv); } else { if (range_max_check (iv.max (), nr) < 0) return tree_constant (); retval = do_matrix_index (iv, 0); } } } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range ri = tmp_i.range_value (); if (len == 2 && is_zero_one (ri)) { if (swap_indices) retval = do_matrix_index (0, 1); else retval = do_matrix_index (1, 0); } else if (len == 2 && is_one_zero (ri)) { retval = do_matrix_index (0, 0); } else { if (index_check (ri, "") < 0) return tree_constant (); if (swap_indices) { if (range_max_check (tree_to_mat_idx (ri.max ()), nc) < 0) return tree_constant (); retval = do_matrix_index (0, ri); } else { if (range_max_check (tree_to_mat_idx (ri.max ()), nr) < 0) return tree_constant (); retval = do_matrix_index (ri, 0); } } } break; case magic_colon: if (swap_indices) retval = do_matrix_index (0, magic_colon); else retval = do_matrix_index (magic_colon, 0); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_matrix_index (int i, const tree_constant& j_arg) const { tree_constant retval; tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int nr = rows (); int nc = columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return tree_constant (); if (range_max_check (i, j, nr, nc) < 0) return tree_constant (); retval = do_matrix_index (i, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); if (! jv) return tree_constant (); if (jv.length () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else { if (range_max_check (i, jv.max (), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (i, jv); } } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (nc == 2 && is_zero_one (rj)) { retval = do_matrix_index (i, 1); } else if (nc == 2 && is_one_zero (rj)) { retval = do_matrix_index (i, 0); } else { if (index_check (rj, "column") < 0) return tree_constant (); if (range_max_check (i, tree_to_mat_idx (rj.max ()), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (i, rj); } } break; case magic_colon: if (range_max_check (i, 0, nr, nc) < 0) return tree_constant (); retval = do_matrix_index (i, magic_colon); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_matrix_index (const idx_vector& iv, const tree_constant& j_arg) const { tree_constant retval; tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int nr = rows (); int nc = columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return tree_constant (); if (range_max_check (iv.max (), j, nr, nc) < 0) return tree_constant (); retval = do_matrix_index (iv, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); if (! jv) return tree_constant (); if (jv.length () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else { if (range_max_check (iv.max (), jv.max (), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (iv, jv); } } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (nc == 2 && is_zero_one (rj)) { retval = do_matrix_index (iv, 1); } else if (nc == 2 && is_one_zero (rj)) { retval = do_matrix_index (iv, 0); } else { if (index_check (rj, "column") < 0) return tree_constant (); if (range_max_check (iv.max (), tree_to_mat_idx (rj.max ()), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (iv, rj); } } break; case magic_colon: if (range_max_check (iv.max (), 0, nr, nc) < 0) return tree_constant (); retval = do_matrix_index (iv, magic_colon); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_matrix_index (const Range& ri, const tree_constant& j_arg) const { tree_constant retval; tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int nr = rows (); int nc = columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return tree_constant (); if (range_max_check (tree_to_mat_idx (ri.max ()), j, nr, nc) < 0) return tree_constant (); retval = do_matrix_index (ri, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); if (! jv) return tree_constant (); if (jv.length () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else { if (range_max_check (tree_to_mat_idx (ri.max ()), jv.max (), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (ri, jv); } } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (nc == 2 && is_zero_one (rj)) { retval = do_matrix_index (ri, 1); } else if (nc == 2 && is_one_zero (rj)) { retval = do_matrix_index (ri, 0); } else { if (index_check (rj, "column") < 0) return tree_constant (); if (range_max_check (tree_to_mat_idx (ri.max ()), tree_to_mat_idx (rj.max ()), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (ri, rj); } } break; case magic_colon: retval = do_matrix_index (ri, magic_colon); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_matrix_index (TC_REP::constant_type mci, const tree_constant& j_arg) const { tree_constant retval; tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int nr = rows (); int nc = columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return tree_constant (); if (range_max_check (0, j, nr, nc) < 0) return tree_constant (); retval = do_matrix_index (magic_colon, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", nc); if (! jv) return tree_constant (); if (jv.length () == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else { if (range_max_check (0, jv.max (), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (magic_colon, jv); } } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (nc == 2 && is_zero_one (rj)) { retval = do_matrix_index (magic_colon, 1); } else if (nc == 2 && is_one_zero (rj)) { retval = do_matrix_index (magic_colon, 0); } else { if (index_check (rj, "column") < 0) return tree_constant (); if (range_max_check (0, tree_to_mat_idx (rj.max ()), nr, nc) < 0) return tree_constant (); retval = do_matrix_index (magic_colon, rj); } } break; case magic_colon: retval = do_matrix_index (magic_colon, magic_colon); break; default: panic_impossible (); break; } return retval; } tree_constant TC_REP::do_matrix_index (int i, int j) const { tree_constant retval; if (type_tag == matrix_constant) retval = tree_constant (matrix->elem (i, j)); else retval = tree_constant (complex_matrix->elem (i, j)); return retval; } tree_constant TC_REP::do_matrix_index (int i, const idx_vector& jv) const { tree_constant retval; int jlen = jv.capacity (); CRMATRIX (m, cm, 1, jlen); for (int j = 0; j < jlen; j++) { int col = jv.elem (j); CRMATRIX_ASSIGN_REP_ELEM (m, cm, 0, j, i, col); } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (int i, const Range& rj) const { tree_constant retval; int jlen = rj.nelem (); CRMATRIX (m, cm, 1, jlen); double b = rj.base (); double increment = rj.inc (); for (int j = 0; j < jlen; j++) { double tmp = b + j * increment; int col = tree_to_mat_idx (tmp); CRMATRIX_ASSIGN_REP_ELEM (m, cm, 0, j, i, col); } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (int i, TC_REP::constant_type mcj) const { assert (mcj == magic_colon); tree_constant retval; int nc = columns (); CRMATRIX (m, cm, 1, nc); for (int j = 0; j < nc; j++) { CRMATRIX_ASSIGN_REP_ELEM (m, cm, 0, j, i, j); } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const idx_vector& iv, int j) const { tree_constant retval; int ilen = iv.capacity (); CRMATRIX (m, cm, ilen, 1); for (int i = 0; i < ilen; i++) { int row = iv.elem (i); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, 0, row, j); } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const idx_vector& iv, const idx_vector& jv) const { tree_constant retval; int ilen = iv.capacity (); int jlen = jv.capacity (); CRMATRIX (m, cm, ilen, jlen); for (int i = 0; i < ilen; i++) { int row = iv.elem (i); for (int j = 0; j < jlen; j++) { int col = jv.elem (j); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const idx_vector& iv, const Range& rj) const { tree_constant retval; int ilen = iv.capacity (); int jlen = rj.nelem (); CRMATRIX (m, cm, ilen, jlen); double b = rj.base (); double increment = rj.inc (); for (int i = 0; i < ilen; i++) { int row = iv.elem (i); for (int j = 0; j < jlen; j++) { double tmp = b + j * increment; int col = tree_to_mat_idx (tmp); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const idx_vector& iv, TC_REP::constant_type mcj) const { assert (mcj == magic_colon); tree_constant retval; int nc = columns (); int ilen = iv.capacity (); CRMATRIX (m, cm, ilen, nc); for (int j = 0; j < nc; j++) { for (int i = 0; i < ilen; i++) { int row = iv.elem (i); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, j); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const Range& ri, int j) const { tree_constant retval; int ilen = ri.nelem (); CRMATRIX (m, cm, ilen, 1); double b = ri.base (); double increment = ri.inc (); for (int i = 0; i < ilen; i++) { double tmp = b + i * increment; int row = tree_to_mat_idx (tmp); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, 0, row, j); } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const Range& ri, const idx_vector& jv) const { tree_constant retval; int ilen = ri.nelem (); int jlen = jv.capacity (); CRMATRIX (m, cm, ilen, jlen); double b = ri.base (); double increment = ri.inc (); for (int i = 0; i < ilen; i++) { double tmp = b + i * increment; int row = tree_to_mat_idx (tmp); for (int j = 0; j < jlen; j++) { int col = jv.elem (j); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const Range& ri, const Range& rj) const { tree_constant retval; int ilen = ri.nelem (); int jlen = rj.nelem (); CRMATRIX (m, cm, ilen, jlen); double ib = ri.base (); double iinc = ri.inc (); double jb = rj.base (); double jinc = rj.inc (); for (int i = 0; i < ilen; i++) { double itmp = ib + i * iinc; int row = tree_to_mat_idx (itmp); for (int j = 0; j < jlen; j++) { double jtmp = jb + j * jinc; int col = tree_to_mat_idx (jtmp); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, col); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (const Range& ri, TC_REP::constant_type mcj) const { assert (mcj == magic_colon); tree_constant retval; int nc = columns (); int ilen = ri.nelem (); CRMATRIX (m, cm, ilen, nc); double ib = ri.base (); double iinc = ri.inc (); for (int i = 0; i < ilen; i++) { double itmp = ib + i * iinc; int row = tree_to_mat_idx (itmp); for (int j = 0; j < nc; j++) { CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, row, j); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (TC_REP::constant_type mci, int j) const { assert (mci == magic_colon); tree_constant retval; int nr = rows (); CRMATRIX (m, cm, nr, 1); for (int i = 0; i < nr; i++) { CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, 0, i, j); } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (TC_REP::constant_type mci, const idx_vector& jv) const { assert (mci == magic_colon); tree_constant retval; int nr = rows (); int jlen = jv.capacity (); CRMATRIX (m, cm, nr, jlen); for (int i = 0; i < nr; i++) { for (int j = 0; j < jlen; j++) { int col = jv.elem (j); CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, i, col); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (TC_REP::constant_type mci, const Range& rj) const { assert (mci == magic_colon); tree_constant retval; int nr = rows (); int jlen = rj.nelem (); CRMATRIX (m, cm, nr, jlen); double jb = rj.base (); double jinc = rj.inc (); for (int j = 0; j < jlen; j++) { double jtmp = jb + j * jinc; int col = tree_to_mat_idx (jtmp); for (int i = 0; i < nr; i++) { CRMATRIX_ASSIGN_REP_ELEM (m, cm, i, j, i, col); } } ASSIGN_CRMATRIX_TO (retval, m, cm); return retval; } tree_constant TC_REP::do_matrix_index (TC_REP::constant_type mci, TC_REP::constant_type mcj) const { tree_constant retval; assert (mci == magic_colon && mcj == magic_colon); switch (type_tag) { case complex_scalar_constant: retval = *complex_scalar; break; case scalar_constant: retval = scalar; break; case complex_matrix_constant: retval = *complex_matrix; break; case matrix_constant: retval = *matrix; break; case range_constant: retval = *range; break; case string_constant: retval = string; break; case magic_colon: default: panic_impossible (); break; } return retval; } // Top-level tree-constant function that handles assignments. Only // decide if the left-hand side is currently a scalar or a matrix and // hand off to other functions to do the real work. void TC_REP::assign (const tree_constant& rhs, const Octave_object& args) { tree_constant rhs_tmp = rhs.make_numeric (); // This is easier than actually handling assignments to strings. // An assignment to a range will normally require a conversion to a // vector since it will normally destroy the equally-spaced property // of the range elements. if (type_tag == string_constant || type_tag == range_constant) force_numeric (); switch (type_tag) { case complex_scalar_constant: case scalar_constant: case unknown_constant: do_scalar_assignment (rhs_tmp, args); break; case complex_matrix_constant: case matrix_constant: do_matrix_assignment (rhs_tmp, args); break; case string_constant: ::error ("invalid assignment to string type"); break; case range_constant: case magic_colon: default: panic_impossible (); break; } } // Assignments to scalars. If resize_on_range_error is true, // this can convert the left-hand side to a matrix. void TC_REP::do_scalar_assignment (const tree_constant& rhs, const Octave_object& args) { assert (type_tag == unknown_constant || type_tag == scalar_constant || type_tag == complex_scalar_constant); int nargin = args.length (); if ((rhs.is_scalar_type () || rhs.is_zero_by_zero ()) && valid_scalar_indices (args)) { if (rhs.is_zero_by_zero ()) { if (type_tag == complex_scalar_constant) delete complex_scalar; matrix = new Matrix (0, 0); type_tag = matrix_constant; } else if (type_tag == unknown_constant || type_tag == scalar_constant) { if (rhs.const_type () == scalar_constant) { scalar = rhs.double_value (); type_tag = scalar_constant; } else if (rhs.const_type () == complex_scalar_constant) { complex_scalar = new Complex (rhs.complex_value ()); type_tag = complex_scalar_constant; } else { ::error ("invalid assignment to scalar"); return; } } else { if (rhs.const_type () == scalar_constant) { delete complex_scalar; scalar = rhs.double_value (); type_tag = scalar_constant; } else if (rhs.const_type () == complex_scalar_constant) { *complex_scalar = rhs.complex_value (); type_tag = complex_scalar_constant; } else { ::error ("invalid assignment to scalar"); return; } } } else if (user_pref.resize_on_range_error) { TC_REP::constant_type old_type_tag = type_tag; if (type_tag == complex_scalar_constant) { Complex *old_complex = complex_scalar; complex_matrix = new ComplexMatrix (1, 1, *complex_scalar); type_tag = complex_matrix_constant; delete old_complex; } else if (type_tag == scalar_constant) { matrix = new Matrix (1, 1, scalar); type_tag = matrix_constant; } // If there is an error, the call to do_matrix_assignment should not // destroy the current value. // TC_REP::eval(int) will take // care of converting single element matrices back to scalars. do_matrix_assignment (rhs, args); // I don't think there's any other way to revert back to unknown // constant types, so here it is. if (old_type_tag == unknown_constant && error_state) { if (type_tag == matrix_constant) delete matrix; else if (type_tag == complex_matrix_constant) delete complex_matrix; type_tag = unknown_constant; } } else if (nargin > 3 || nargin < 2) ::error ("invalid index expression for scalar type"); else ::error ("index invalid or out of range for scalar type"); } // Assignments to matrices (and vectors). // // For compatibility with Matlab, we allow assignment of an empty // matrix to an expression with empty indices to do nothing. void TC_REP::do_matrix_assignment (const tree_constant& rhs, const Octave_object& args) { assert (type_tag == unknown_constant || type_tag == matrix_constant || type_tag == complex_matrix_constant); if (type_tag == matrix_constant && rhs.is_complex_type ()) { Matrix *old_matrix = matrix; complex_matrix = new ComplexMatrix (*matrix); type_tag = complex_matrix_constant; delete old_matrix; } else if (type_tag == unknown_constant) { if (rhs.is_complex_type ()) { complex_matrix = new ComplexMatrix (); type_tag = complex_matrix_constant; } else { matrix = new Matrix (); type_tag = matrix_constant; } } int nargin = args.length (); // The do_matrix_assignment functions can't handle empty matrices, so // don't let any pass through here. switch (nargin) { case 2: if (args.length () <= 0) ::error ("matrix index is null"); else if (args(1).is_undefined ()) ::error ("matrix index is undefined"); else do_matrix_assignment (rhs, args(1)); break; case 3: if (args.length () <= 0) ::error ("matrix indices are null"); else if (args(1).is_undefined ()) ::error ("first matrix index is undefined"); else if (args(2).is_undefined ()) ::error ("second matrix index is undefined"); else if (args(1).is_empty () || args(2).is_empty ()) { if (! rhs.is_empty ()) { ::error ("in assignment expression, a matrix index is empty"); ::error ("but hte right hand side is not an empty matrix"); } // XXX FIXME XXX -- to really be correct here, we should probably // check to see if the assignment conforms, but that seems like more // work than it's worth right now... } else do_matrix_assignment (rhs, args(1), args(2)); break; default: ::error ("too many indices for matrix expression"); break; } } // Matrix assignments indexed by a single value. void TC_REP::do_matrix_assignment (const tree_constant& rhs, const tree_constant& i_arg) { int nr = rows (); int nc = columns (); if (user_pref.do_fortran_indexing || nr <= 1 || nc <= 1) { if (i_arg.is_empty ()) { if (! rhs.is_empty ()) { ::error ("in assignment expression, matrix index is empty but"); ::error ("right hand side is not an empty matrix"); } // XXX FIXME XXX -- to really be correct here, we should probably // check to see if the assignment conforms, but that seems like more // work than it's worth right now... // The assignment functions can't handle empty matrices, so don't let // any pass through here. return; } // We can't handle the case of assigning to a vector first, since even // then, the two operations are not equivalent. For example, the // expression V(:) = M is handled differently depending on whether the // user specified do_fortran_indexing = "true". if (user_pref.do_fortran_indexing) fortran_style_matrix_assignment (rhs, i_arg); else if (nr <= 1 || nc <= 1) vector_assignment (rhs, i_arg); else panic_impossible (); } else ::error ("single index only valid for row or column vector"); } // Fortran-style assignments. Matrices are assumed to be stored in // column-major order and it is ok to use a single index for // multi-dimensional matrices. void TC_REP::fortran_style_matrix_assignment (const tree_constant& rhs, const tree_constant& i_arg) { tree_constant tmp_i = i_arg.make_numeric_or_magic (); TC_REP::constant_type itype = tmp_i.const_type (); int nr = rows (); int nc = columns (); int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); switch (itype) { case complex_scalar_constant: case scalar_constant: { int i = NINT (tmp_i.double_value ()); int idx = i - 1; if (rhs_nr == 0 && rhs_nc == 0) { if (idx < nr * nc) { convert_to_row_or_column_vector (); nr = rows (); nc = columns (); if (nr == 1) delete_column (idx); else if (nc == 1) delete_row (idx); else panic_impossible (); } return; } if (index_check (idx, "") < 0) return; if (nr <= 1 || nc <= 1) { maybe_resize (idx); if (error_state) return; } else if (range_max_check (idx, nr * nc) < 0) return; nr = rows (); nc = columns (); if (! indexed_assign_conforms (1, 1, rhs_nr, rhs_nc)) { ::error ("for A(int) = X: X must be a scalar"); return; } int ii = fortran_row (i, nr) - 1; int jj = fortran_column (i, nr) - 1; do_matrix_assignment (rhs, ii, jj); } break; case complex_matrix_constant: case matrix_constant: { Matrix mi = tmp_i.matrix_value (); int len = nr * nc; idx_vector ii (mi, 1, "", len); // Always do fortran indexing here... if (! ii) return; if (rhs_nr == 0 && rhs_nc == 0) { ii.sort_uniq (); int num_to_delete = 0; for (int i = 0; i < ii.length (); i++) { if (ii.elem (i) < len) num_to_delete++; else break; } if (num_to_delete > 0) { if (num_to_delete != ii.length ()) ii.shorten (num_to_delete); convert_to_row_or_column_vector (); nr = rows (); nc = columns (); if (nr == 1) delete_columns (ii); else if (nc == 1) delete_rows (ii); else panic_impossible (); } return; } if (nr <= 1 || nc <= 1) { maybe_resize (ii.max ()); if (error_state) return; } else if (range_max_check (ii.max (), len) < 0) return; int ilen = ii.capacity (); if (ilen != rhs_nr * rhs_nc) { ::error ("A(matrix) = X: X and matrix must have the same number"); ::error ("of elements"); } else if (ilen == 1 && rhs.is_scalar_type ()) { int nr = rows (); int idx = ii.elem (0); int ii = fortran_row (idx + 1, nr) - 1; int jj = fortran_column (idx + 1, nr) - 1; if (rhs.const_type () == scalar_constant) matrix->elem (ii, jj) = rhs.double_value (); else if (rhs.const_type () == complex_scalar_constant) complex_matrix->elem (ii, jj) = rhs.complex_value (); else panic_impossible (); } else fortran_style_matrix_assignment (rhs, ii); } break; case string_constant: gripe_string_invalid (); break; case range_constant: gripe_range_invalid (); break; case magic_colon: // a(:) = [] is equivalent to a(:,:) = []. if (rhs_nr == 0 && rhs_nc == 0) do_matrix_assignment (rhs, magic_colon, magic_colon); else fortran_style_matrix_assignment (rhs, magic_colon); break; default: panic_impossible (); break; } } // Fortran-style assignment for vector index. void TC_REP::fortran_style_matrix_assignment (const tree_constant& rhs, idx_vector& i) { assert (rhs.is_matrix_type ()); int ilen = i.capacity (); REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int len = rhs_nr * rhs_nc; if (len == ilen) { int nr = rows (); if (rhs.const_type () == matrix_constant) { double *cop_out = rhs_m.fortran_vec (); for (int k = 0; k < len; k++) { int ii = fortran_row (i.elem (k) + 1, nr) - 1; int jj = fortran_column (i.elem (k) + 1, nr) - 1; matrix->elem (ii, jj) = *cop_out++; } } else { Complex *cop_out = rhs_cm.fortran_vec (); for (int k = 0; k < len; k++) { int ii = fortran_row (i.elem (k) + 1, nr) - 1; int jj = fortran_column (i.elem (k) + 1, nr) - 1; complex_matrix->elem (ii, jj) = *cop_out++; } } } else ::error ("number of rows and columns must match for indexed assignment"); } // Fortran-style assignment for colon index. void TC_REP::fortran_style_matrix_assignment (const tree_constant& rhs, TC_REP::constant_type mci) { assert (rhs.is_matrix_type () && mci == TC_REP::magic_colon); int nr = rows (); int nc = columns (); REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int rhs_size = rhs_nr * rhs_nc; if (rhs_size == 0) { if (rhs.const_type () == matrix_constant) { delete matrix; matrix = new Matrix (0, 0); return; } else panic_impossible (); } else if (nr*nc != rhs_size) { ::error ("A(:) = X: X and A must have the same number of elements"); return; } if (rhs.const_type () == matrix_constant) { double *cop_out = rhs_m.fortran_vec (); for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) matrix->elem (i, j) = *cop_out++; } else { Complex *cop_out = rhs_cm.fortran_vec (); for (int j = 0; j < nc; j++) for (int i = 0; i < nr; i++) complex_matrix->elem (i, j) = *cop_out++; } } // Assignments to vectors. Hand off to other functions once we know // what kind of index we have. For a colon, it is the same as // assignment to a matrix indexed by two colons. void TC_REP::vector_assignment (const tree_constant& rhs, const tree_constant& i_arg) { int nr = rows (); int nc = columns (); assert ((nr == 1 || nc == 1 || (nr == 0 && nc == 0)) && ! user_pref.do_fortran_indexing); tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type itype = tmp_i.const_type (); switch (itype) { case complex_scalar_constant: case scalar_constant: { int i = tree_to_mat_idx (tmp_i.double_value ()); if (index_check (i, "") < 0) return; do_vector_assign (rhs, i); } break; case complex_matrix_constant: case matrix_constant: { Matrix mi = tmp_i.matrix_value (); int len = nr * nc; idx_vector iv (mi, user_pref.do_fortran_indexing, "", len); if (! iv) return; do_vector_assign (rhs, iv); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range ri = tmp_i.range_value (); int len = nr * nc; if (len == 2 && is_zero_one (ri)) { do_vector_assign (rhs, 1); } else if (len == 2 && is_one_zero (ri)) { do_vector_assign (rhs, 0); } else { if (index_check (ri, "") < 0) return; do_vector_assign (rhs, ri); } } break; case magic_colon: { int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); if (! indexed_assign_conforms (nr, nc, rhs_nr, rhs_nc)) { ::error ("A(:) = X: X and A must have the same dimensions"); return; } do_matrix_assignment (rhs, magic_colon, magic_colon); } break; default: panic_impossible (); break; } } // Check whether an indexed assignment to a vector is valid. void TC_REP::check_vector_assign (int rhs_nr, int rhs_nc, int ilen, const char *rm) { int nr = rows (); int nc = columns (); if ((nr == 1 && nc == 1) || nr == 0 || nc == 0) // No orientation. { if (! (ilen == rhs_nr || ilen == rhs_nc)) { ::error ("A(%s) = X: X and %s must have the same number of elements", rm, rm); } } else if (nr == 1) // Preserve current row orientation. { if (! (rhs_nr == 1 && rhs_nc == ilen)) { ::error ("A(%s) = X: where A is a row vector, X must also be a", rm); ::error ("row vector with the same number of elements as %s", rm); } } else if (nc == 1) // Preserve current column orientation. { if (! (rhs_nc == 1 && rhs_nr == ilen)) { ::error ("A(%s) = X: where A is a column vector, X must also be", rm); ::error ("a column vector with the same number of elements as %s", rm); } } else panic_impossible (); } // Assignment to a vector with an integer index. void TC_REP::do_vector_assign (const tree_constant& rhs, int i) { int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); if (indexed_assign_conforms (1, 1, rhs_nr, rhs_nc)) { maybe_resize (i); if (error_state) return; int nr = rows (); int nc = columns (); if (nr == 1) { REP_ELEM_ASSIGN (0, i, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } else if (nc == 1) { REP_ELEM_ASSIGN (i, 0, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } else panic_impossible (); } else if (rhs_nr == 0 && rhs_nc == 0) { int nr = rows (); int nc = columns (); int len = MAX (nr, nc); if (i < 0 || i >= len) { ::error ("A(int) = []: index out of range"); return; } if (nr == 1) delete_column (i); else if (nc == 1) delete_row (i); else panic_impossible (); } else { ::error ("for A(int) = X: X must be a scalar"); return; } } // Assignment to a vector with a vector index. void TC_REP::do_vector_assign (const tree_constant& rhs, idx_vector& iv) { if (rhs.is_zero_by_zero ()) { int nr = rows (); int nc = columns (); int len = MAX (nr, nc); if (iv.max () >= len) { ::error ("A(matrix) = []: index out of range"); return; } if (nr == 1) delete_columns (iv); else if (nc == 1) delete_rows (iv); else panic_impossible (); } else if (rhs.is_scalar_type ()) { int nr = rows (); int nc = columns (); if (iv.capacity () == 1) { int idx = iv.elem (0); if (nr == 1) { REP_ELEM_ASSIGN (0, idx, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } else if (nc == 1) { REP_ELEM_ASSIGN (idx, 0, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } else panic_impossible (); } else { if (nr == 1) { ::error ("A(matrix) = X: where A is a row vector, X must also be a"); ::error ("row vector with the same number of elements as matrix"); } else if (nc == 1) { ::error ("A(matrix) = X: where A is a column vector, X must also be a"); ::error ("column vector with the same number of elements as matrix"); } else panic_impossible (); } } else if (rhs.is_matrix_type ()) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int ilen = iv.capacity (); check_vector_assign (rhs_nr, rhs_nc, ilen, "matrix"); if (error_state) return; force_orient f_orient = no_orient; if (rhs_nr == 1 && rhs_nc != 1) f_orient = row_orient; else if (rhs_nc == 1 && rhs_nr != 1) f_orient = column_orient; maybe_resize (iv.max (), f_orient); if (error_state) return; int nr = rows (); int nc = columns (); if (nr == 1) { for (int i = 0; i < iv.capacity (); i++) REP_ELEM_ASSIGN (0, iv.elem (i), rhs_m.elem (0, i), rhs_cm.elem (0, i), rhs.is_real_type ()); } else if (nc == 1) { for (int i = 0; i < iv.capacity (); i++) REP_ELEM_ASSIGN (iv.elem (i), 0, rhs_m.elem (i, 0), rhs_cm.elem (i, 0), rhs.is_real_type ()); } else panic_impossible (); } else panic_impossible (); } // Assignment to a vector with a range index. void TC_REP::do_vector_assign (const tree_constant& rhs, Range& ri) { if (rhs.is_zero_by_zero ()) { int nr = rows (); int nc = columns (); int len = MAX (nr, nc); int b = tree_to_mat_idx (ri.min ()); int l = tree_to_mat_idx (ri.max ()); if (b < 0 || l >= len) { ::error ("A(range) = []: index out of range"); return; } if (nr == 1) delete_columns (ri); else if (nc == 1) delete_rows (ri); else panic_impossible (); } else if (rhs.is_scalar_type ()) { int nr = rows (); int nc = columns (); if (nr == 1) { ::error ("A(range) = X: where A is a row vector, X must also be a"); ::error ("row vector with the same number of elements as range"); } else if (nc == 1) { ::error ("A(range) = X: where A is a column vector, X must also be a"); ::error ("column vector with the same number of elements as range"); } else panic_impossible (); } else if (rhs.is_matrix_type ()) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int ilen = ri.nelem (); check_vector_assign (rhs_nr, rhs_nc, ilen, "range"); if (error_state) return; force_orient f_orient = no_orient; if (rhs_nr == 1 && rhs_nc != 1) f_orient = row_orient; else if (rhs_nc == 1 && rhs_nr != 1) f_orient = column_orient; maybe_resize (tree_to_mat_idx (ri.max ()), f_orient); if (error_state) return; int nr = rows (); int nc = columns (); double b = ri.base (); double increment = ri.inc (); if (nr == 1) { for (int i = 0; i < ri.nelem (); i++) { double tmp = b + i * increment; int col = tree_to_mat_idx (tmp); REP_ELEM_ASSIGN (0, col, rhs_m.elem (0, i), rhs_cm.elem (0, i), rhs.is_real_type ()); } } else if (nc == 1) { for (int i = 0; i < ri.nelem (); i++) { double tmp = b + i * increment; int row = tree_to_mat_idx (tmp); REP_ELEM_ASSIGN (row, 0, rhs_m.elem (i, 0), rhs_cm.elem (i, 0), rhs.is_real_type ()); } } else panic_impossible (); } else panic_impossible (); } // Matrix assignment indexed by two values. This function determines // the type of the first arugment, checks as much as possible, and // then calls one of a set of functions to handle the specific cases: // // M (integer, arg2) = RHS (MA1) // M (vector, arg2) = RHS (MA2) // M (range, arg2) = RHS (MA3) // M (colon, arg2) = RHS (MA4) // // Each of those functions determines the type of the second argument // and calls another function to handle the real work of doing the // assignment. void TC_REP::do_matrix_assignment (const tree_constant& rhs, const tree_constant& i_arg, const tree_constant& j_arg) { tree_constant tmp_i = i_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type itype = tmp_i.const_type (); switch (itype) { case complex_scalar_constant: case scalar_constant: { int i = tree_to_mat_idx (tmp_i.double_value ()); if (index_check (i, "row") < 0) return; do_matrix_assignment (rhs, i, j_arg); } break; case complex_matrix_constant: case matrix_constant: { Matrix mi = tmp_i.matrix_value (); idx_vector iv (mi, user_pref.do_fortran_indexing, "row", rows ()); if (! iv) return; do_matrix_assignment (rhs, iv, j_arg); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range ri = tmp_i.range_value (); int nr = rows (); if (nr == 2 && is_zero_one (ri)) { do_matrix_assignment (rhs, 1, j_arg); } else if (nr == 2 && is_one_zero (ri)) { do_matrix_assignment (rhs, 0, j_arg); } else { if (index_check (ri, "row") < 0) return; do_matrix_assignment (rhs, ri, j_arg); } } break; case magic_colon: do_matrix_assignment (rhs, magic_colon, j_arg); break; default: panic_impossible (); break; } } /* MA1 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, int i, const tree_constant& j_arg) { tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return; if (! indexed_assign_conforms (1, 1, rhs_nr, rhs_nc)) { ::error ("A(int,int) = X, X must be a scalar"); return; } maybe_resize (i, j); if (error_state) return; do_matrix_assignment (rhs, i, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", columns ()); if (! jv) return; if (! indexed_assign_conforms (1, jv.capacity (), rhs_nr, rhs_nc)) { ::error ("A(int,matrix) = X: X must be a row vector with the same"); ::error ("number of elements as matrix"); return; } maybe_resize (i, jv.max ()); if (error_state) return; do_matrix_assignment (rhs, i, jv); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (! indexed_assign_conforms (1, rj.nelem (), rhs_nr, rhs_nc)) { ::error ("A(int,range) = X: X must be a row vector with the same"); ::error ("number of elements as range"); return; } int nc = columns (); if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, i, 1); } else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, i, 0); } else { if (index_check (rj, "column") < 0) return; maybe_resize (i, tree_to_mat_idx (rj.max ())); if (error_state) return; do_matrix_assignment (rhs, i, rj); } } break; case magic_colon: { int nc = columns (); int nr = rows (); if (nc == 0 && nr == 0 && rhs_nr == 1) { if (rhs.is_complex_type ()) { complex_matrix = new ComplexMatrix (); type_tag = complex_matrix_constant; } else { matrix = new Matrix (); type_tag = matrix_constant; } maybe_resize (i, rhs_nc-1); if (error_state) return; } else if (indexed_assign_conforms (1, nc, rhs_nr, rhs_nc)) { maybe_resize (i, nc-1); if (error_state) return; } else if (rhs_nr == 0 && rhs_nc == 0) { if (i < 0 || i >= nr) { ::error ("A(int,:) = []: row index out of range"); return; } } else { ::error ("A(int,:) = X: X must be a row vector with the same"); ::error ("number of columns as A"); return; } do_matrix_assignment (rhs, i, magic_colon); } break; default: panic_impossible (); break; } } /* MA2 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, idx_vector& iv, const tree_constant& j_arg) { tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return; if (! indexed_assign_conforms (iv.capacity (), 1, rhs_nr, rhs_nc)) { ::error ("A(matrix,int) = X: X must be a column vector with the"); ::error ("same number of elements as matrix"); return; } maybe_resize (iv.max (), j); if (error_state) return; do_matrix_assignment (rhs, iv, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", columns ()); if (! jv) return; if (! indexed_assign_conforms (iv.capacity (), jv.capacity (), rhs_nr, rhs_nc)) { ::error ("A(r_mat,c_mat) = X: the number of rows in X must match"); ::error ("the number of elements in r_mat and the number of"); ::error ("columns in X must match the number of elements in c_mat"); return; } maybe_resize (iv.max (), jv.max ()); if (error_state) return; do_matrix_assignment (rhs, iv, jv); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (! indexed_assign_conforms (iv.capacity (), rj.nelem (), rhs_nr, rhs_nc)) { ::error ("A(matrix,range) = X: the number of rows in X must match"); ::error ("the number of elements in matrix and the number of"); ::error ("columns in X must match the number of elements in range"); return; } int nc = columns (); if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, iv, 1); } else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, iv, 0); } else { if (index_check (rj, "column") < 0) return; maybe_resize (iv.max (), tree_to_mat_idx (rj.max ())); if (error_state) return; do_matrix_assignment (rhs, iv, rj); } } break; case magic_colon: { int nc = columns (); int new_nc = nc; if (nc == 0) new_nc = rhs_nc; if (indexed_assign_conforms (iv.capacity (), new_nc, rhs_nr, rhs_nc)) { maybe_resize (iv.max (), new_nc-1); if (error_state) return; } else if (rhs_nr == 0 && rhs_nc == 0) { if (iv.max () >= rows ()) { ::error ("A(matrix,:) = []: row index out of range"); return; } } else { ::error ("A(matrix,:) = X: the number of rows in X must match the"); ::error ("number of elements in matrix, and the number of columns"); ::error ("in X must match the number of columns in A"); return; } do_matrix_assignment (rhs, iv, magic_colon); } break; default: panic_impossible (); break; } } /* MA3 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, Range& ri, const tree_constant& j_arg) { tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return; if (! indexed_assign_conforms (ri.nelem (), 1, rhs_nr, rhs_nc)) { ::error ("A(range,int) = X: X must be a column vector with the"); ::error ("same number of elements as range"); return; } maybe_resize (tree_to_mat_idx (ri.max ()), j); if (error_state) return; do_matrix_assignment (rhs, ri, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", columns ()); if (! jv) return; if (! indexed_assign_conforms (ri.nelem (), jv.capacity (), rhs_nr, rhs_nc)) { ::error ("A(range,matrix) = X: the number of rows in X must match"); ::error ("the number of elements in range and the number of"); ::error ("columns in X must match the number of elements in matrix"); return; } maybe_resize (tree_to_mat_idx (ri.max ()), jv.max ()); if (error_state) return; do_matrix_assignment (rhs, ri, jv); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); if (! indexed_assign_conforms (ri.nelem (), rj.nelem (), rhs_nr, rhs_nc)) { ::error ("A(r_range,c_range) = X: the number of rows in X must"); ::error ("match the number of elements in r_range and the number"); ::error ("of columns in X must match the number of elements in"); ::error ("c_range"); return; } int nc = columns (); if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, ri, 1); } else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, ri, 0); } else { if (index_check (rj, "column") < 0) return; maybe_resize (tree_to_mat_idx (ri.max ()), tree_to_mat_idx (rj.max ())); if (error_state) return; do_matrix_assignment (rhs, ri, rj); } } break; case magic_colon: { int nc = columns (); int new_nc = nc; if (nc == 0) new_nc = rhs_nc; if (indexed_assign_conforms (ri.nelem (), new_nc, rhs_nr, rhs_nc)) { maybe_resize (tree_to_mat_idx (ri.max ()), new_nc-1); if (error_state) return; } else if (rhs_nr == 0 && rhs_nc == 0) { int b = tree_to_mat_idx (ri.min ()); int l = tree_to_mat_idx (ri.max ()); if (b < 0 || l >= rows ()) { ::error ("A(range,:) = []: row index out of range"); return; } } else { ::error ("A(range,:) = X: the number of rows in X must match the"); ::error ("number of elements in range, and the number of columns"); ::error ("in X must match the number of columns in A"); return; } do_matrix_assignment (rhs, ri, magic_colon); } break; default: panic_impossible (); break; } } /* MA4 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, TC_REP::constant_type i, const tree_constant& j_arg) { tree_constant tmp_j = j_arg.make_numeric_or_range_or_magic (); TC_REP::constant_type jtype = tmp_j.const_type (); int rhs_nr = rhs.rows (); int rhs_nc = rhs.columns (); switch (jtype) { case complex_scalar_constant: case scalar_constant: { int j = tree_to_mat_idx (tmp_j.double_value ()); if (index_check (j, "column") < 0) return; int nr = rows (); int nc = columns (); if (nr == 0 && nc == 0 && rhs_nc == 1) { if (rhs.is_complex_type ()) { complex_matrix = new ComplexMatrix (); type_tag = complex_matrix_constant; } else { matrix = new Matrix (); type_tag = matrix_constant; } maybe_resize (rhs_nr-1, j); if (error_state) return; } else if (indexed_assign_conforms (nr, 1, rhs_nr, rhs_nc)) { maybe_resize (nr-1, j); if (error_state) return; } else if (rhs_nr == 0 && rhs_nc == 0) { if (j < 0 || j >= nc) { ::error ("A(:,int) = []: column index out of range"); return; } } else { ::error ("A(:,int) = X: X must be a column vector with the same"); ::error ("number of rows as A"); return; } do_matrix_assignment (rhs, magic_colon, j); } break; case complex_matrix_constant: case matrix_constant: { Matrix mj = tmp_j.matrix_value (); idx_vector jv (mj, user_pref.do_fortran_indexing, "column", columns ()); if (! jv) return; int nr = rows (); int new_nr = nr; if (nr == 0) new_nr = rhs_nr; if (indexed_assign_conforms (new_nr, jv.capacity (), rhs_nr, rhs_nc)) { maybe_resize (new_nr-1, jv.max ()); if (error_state) return; } else if (rhs_nr == 0 && rhs_nc == 0) { if (jv.max () >= columns ()) { ::error ("A(:,matrix) = []: column index out of range"); return; } } else { ::error ("A(:,matrix) = X: the number of rows in X must match the"); ::error ("number of rows in A, and the number of columns in X must"); ::error ("match the number of elements in matrix"); return; } do_matrix_assignment (rhs, magic_colon, jv); } break; case string_constant: gripe_string_invalid (); break; case range_constant: { Range rj = tmp_j.range_value (); int nr = rows (); int new_nr = nr; if (nr == 0) new_nr = rhs_nr; if (indexed_assign_conforms (new_nr, rj.nelem (), rhs_nr, rhs_nc)) { int nc = columns (); if (nc == 2 && is_zero_one (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, magic_colon, 1); } else if (nc == 2 && is_one_zero (rj) && rhs_nc == 1) { do_matrix_assignment (rhs, magic_colon, 0); } else { if (index_check (rj, "column") < 0) return; maybe_resize (new_nr-1, tree_to_mat_idx (rj.max ())); if (error_state) return; } } else if (rhs_nr == 0 && rhs_nc == 0) { int b = tree_to_mat_idx (rj.min ()); int l = tree_to_mat_idx (rj.max ()); if (b < 0 || l >= columns ()) { ::error ("A(:,range) = []: column index out of range"); return; } } else { ::error ("A(:,range) = X: the number of rows in X must match the"); ::error ("number of rows in A, and the number of columns in X"); ::error ("must match the number of elements in range"); return; } do_matrix_assignment (rhs, magic_colon, rj); } break; case magic_colon: // a(:,:) = foo is equivalent to a = foo. do_matrix_assignment (rhs, magic_colon, magic_colon); break; default: panic_impossible (); break; } } // Functions that actually handle assignment to a matrix using two // index values. // // idx2 // +---+---+----+----+ // idx1 | i | v | r | c | // ---------+---+---+----+----+ // integer | 1 | 5 | 9 | 13 | // ---------+---+---+----+----+ // vector | 2 | 6 | 10 | 14 | // ---------+---+---+----+----+ // range | 3 | 7 | 11 | 15 | // ---------+---+---+----+----+ // colon | 4 | 8 | 12 | 16 | // ---------+---+---+----+----+ /* 1 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, int i, int j) { REP_ELEM_ASSIGN (i, j, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } /* 2 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, int i, idx_vector& jv) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); for (int j = 0; j < jv.capacity (); j++) REP_ELEM_ASSIGN (i, jv.elem (j), rhs_m.elem (0, j), rhs_cm.elem (0, j), rhs.is_real_type ()); } /* 3 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, int i, Range& rj) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); double b = rj.base (); double increment = rj.inc (); for (int j = 0; j < rj.nelem (); j++) { double tmp = b + j * increment; int col = tree_to_mat_idx (tmp); REP_ELEM_ASSIGN (i, col, rhs_m.elem (0, j), rhs_cm.elem (0, j), rhs.is_real_type ()); } } /* 4 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, int i, TC_REP::constant_type mcj) { assert (mcj == magic_colon); int nc = columns (); if (rhs.is_zero_by_zero ()) { delete_row (i); } else if (rhs.is_matrix_type ()) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); for (int j = 0; j < nc; j++) REP_ELEM_ASSIGN (i, j, rhs_m.elem (0, j), rhs_cm.elem (0, j), rhs.is_real_type ()); } else if (rhs.is_scalar_type () && nc == 1) { REP_ELEM_ASSIGN (i, 0, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } else panic_impossible (); } /* 5 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, idx_vector& iv, int j) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); for (int i = 0; i < iv.capacity (); i++) { int row = iv.elem (i); REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, 0), rhs_cm.elem (i, 0), rhs.is_real_type ()); } } /* 6 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, idx_vector& iv, idx_vector& jv) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); for (int i = 0; i < iv.capacity (); i++) { int row = iv.elem (i); for (int j = 0; j < jv.capacity (); j++) { int col = jv.elem (j); REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } /* 7 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, idx_vector& iv, Range& rj) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); double b = rj.base (); double increment = rj.inc (); for (int i = 0; i < iv.capacity (); i++) { int row = iv.elem (i); for (int j = 0; j < rj.nelem (); j++) { double tmp = b + j * increment; int col = tree_to_mat_idx (tmp); REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } /* 8 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, idx_vector& iv, TC_REP::constant_type mcj) { assert (mcj == magic_colon); if (rhs.is_zero_by_zero ()) { delete_rows (iv); } else { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int nc = columns (); for (int j = 0; j < nc; j++) { for (int i = 0; i < iv.capacity (); i++) { int row = iv.elem (i); REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } } /* 9 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, Range& ri, int j) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); double b = ri.base (); double increment = ri.inc (); for (int i = 0; i < ri.nelem (); i++) { double tmp = b + i * increment; int row = tree_to_mat_idx (tmp); REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, 0), rhs_cm.elem (i, 0), rhs.is_real_type ()); } } /* 10 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, Range& ri, idx_vector& jv) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); double b = ri.base (); double increment = ri.inc (); for (int j = 0; j < jv.capacity (); j++) { int col = jv.elem (j); for (int i = 0; i < ri.nelem (); i++) { double tmp = b + i * increment; int row = tree_to_mat_idx (tmp); REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), rhs_m.elem (i, j), rhs.is_real_type ()); } } } /* 11 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, Range& ri, Range& rj) { double ib = ri.base (); double iinc = ri.inc (); double jb = rj.base (); double jinc = rj.inc (); REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); for (int i = 0; i < ri.nelem (); i++) { double itmp = ib + i * iinc; int row = tree_to_mat_idx (itmp); for (int j = 0; j < rj.nelem (); j++) { double jtmp = jb + j * jinc; int col = tree_to_mat_idx (jtmp); REP_ELEM_ASSIGN (row, col, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } /* 12 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, Range& ri, TC_REP::constant_type mcj) { assert (mcj == magic_colon); if (rhs.is_zero_by_zero ()) { delete_rows (ri); } else { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); double ib = ri.base (); double iinc = ri.inc (); int nc = columns (); for (int i = 0; i < ri.nelem (); i++) { double itmp = ib + i * iinc; int row = tree_to_mat_idx (itmp); for (int j = 0; j < nc; j++) REP_ELEM_ASSIGN (row, j, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } /* 13 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, TC_REP::constant_type mci, int j) { assert (mci == magic_colon); int nr = rows (); if (rhs.is_zero_by_zero ()) { delete_column (j); } else if (rhs.is_matrix_type ()) { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); for (int i = 0; i < nr; i++) REP_ELEM_ASSIGN (i, j, rhs_m.elem (i, 0), rhs_cm.elem (i, 0), rhs.is_real_type ()); } else if (rhs.is_scalar_type () && nr == 1) { REP_ELEM_ASSIGN (0, j, rhs.double_value (), rhs.complex_value (), rhs.is_real_type ()); } else panic_impossible (); } /* 14 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, TC_REP::constant_type mci, idx_vector& jv) { assert (mci == magic_colon); if (rhs.is_zero_by_zero ()) { delete_columns (jv); } else { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int nr = rows (); for (int i = 0; i < nr; i++) { for (int j = 0; j < jv.capacity (); j++) { int col = jv.elem (j); REP_ELEM_ASSIGN (i, col, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } } /* 15 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, TC_REP::constant_type mci, Range& rj) { assert (mci == magic_colon); if (rhs.is_zero_by_zero ()) { delete_columns (rj); } else { REP_RHS_MATRIX (rhs, rhs_m, rhs_cm, rhs_nr, rhs_nc); int nr = rows (); double jb = rj.base (); double jinc = rj.inc (); for (int j = 0; j < rj.nelem (); j++) { double jtmp = jb + j * jinc; int col = tree_to_mat_idx (jtmp); for (int i = 0; i < nr; i++) { REP_ELEM_ASSIGN (i, col, rhs_m.elem (i, j), rhs_cm.elem (i, j), rhs.is_real_type ()); } } } } /* 16 */ void TC_REP::do_matrix_assignment (const tree_constant& rhs, TC_REP::constant_type mci, TC_REP::constant_type mcj) { assert (mci == magic_colon && mcj == magic_colon); switch (type_tag) { case scalar_constant: break; case matrix_constant: delete matrix; break; case complex_scalar_constant: delete complex_scalar; break; case complex_matrix_constant: delete complex_matrix; break; case string_constant: delete [] string; break; case range_constant: delete range; break; case magic_colon: default: panic_impossible (); break; } type_tag = rhs.const_type (); switch (type_tag) { case scalar_constant: scalar = rhs.double_value (); break; case matrix_constant: matrix = new Matrix (rhs.matrix_value ()); break; case string_constant: string = strsave (rhs.string_value ()); break; case complex_matrix_constant: complex_matrix = new ComplexMatrix (rhs.complex_matrix_value ()); break; case complex_scalar_constant: complex_scalar = new Complex (rhs.complex_value ()); break; case range_constant: range = new Range (rhs.range_value ()); break; case magic_colon: default: panic_impossible (); break; } } // Functions for deleting rows or columns of a matrix. These are used // to handle statements like // // M (i, j) = [] void TC_REP::delete_row (int idx) { if (type_tag == matrix_constant) { int nr = matrix->rows (); int nc = matrix->columns (); Matrix *new_matrix = new Matrix (nr-1, nc); int ii = 0; for (int i = 0; i < nr; i++) { if (i != idx) { for (int j = 0; j < nc; j++) new_matrix->elem (ii, j) = matrix->elem (i, j); ii++; } } delete matrix; matrix = new_matrix; } else if (type_tag == complex_matrix_constant) { int nr = complex_matrix->rows (); int nc = complex_matrix->columns (); ComplexMatrix *new_matrix = new ComplexMatrix (nr-1, nc); int ii = 0; for (int i = 0; i < nr; i++) { if (i != idx) { for (int j = 0; j < nc; j++) new_matrix->elem (ii, j) = complex_matrix->elem (i, j); ii++; } } delete complex_matrix; complex_matrix = new_matrix; } else panic_impossible (); } void TC_REP::delete_rows (idx_vector& iv) { iv.sort_uniq (); int num_to_delete = iv.length (); int nr = rows (); int nc = columns (); // If deleting all rows of a column vector, make result 0x0. if (nc == 1 && num_to_delete == nr) nc = 0; if (type_tag == matrix_constant) { Matrix *new_matrix = new Matrix (nr-num_to_delete, nc); if (nr > num_to_delete) { int ii = 0; int idx = 0; for (int i = 0; i < nr; i++) { if (i == iv.elem (idx)) idx++; else { for (int j = 0; j < nc; j++) new_matrix->elem (ii, j) = matrix->elem (i, j); ii++; } } } delete matrix; matrix = new_matrix; } else if (type_tag == complex_matrix_constant) { ComplexMatrix *new_matrix = new ComplexMatrix (nr-num_to_delete, nc); if (nr > num_to_delete) { int ii = 0; int idx = 0; for (int i = 0; i < nr; i++) { if (i == iv.elem (idx)) idx++; else { for (int j = 0; j < nc; j++) new_matrix->elem (ii, j) = complex_matrix->elem (i, j); ii++; } } } delete complex_matrix; complex_matrix = new_matrix; } else panic_impossible (); } void TC_REP::delete_rows (Range& ri) { ri.sort (); int num_to_delete = ri.nelem (); int nr = rows (); int nc = columns (); // If deleting all rows of a column vector, make result 0x0. if (nc == 1 && num_to_delete == nr) nc = 0; double ib = ri.base (); double iinc = ri.inc (); int max_idx = tree_to_mat_idx (ri.max ()); if (type_tag == matrix_constant) { Matrix *new_matrix = new Matrix (nr-num_to_delete, nc); if (nr > num_to_delete) { int ii = 0; int idx = 0; for (int i = 0; i < nr; i++) { double itmp = ib + idx * iinc; int row = tree_to_mat_idx (itmp); if (i == row && row <= max_idx) idx++; else { for (int j = 0; j < nc; j++) new_matrix->elem (ii, j) = matrix->elem (i, j); ii++; } } } delete matrix; matrix = new_matrix; } else if (type_tag == complex_matrix_constant) { ComplexMatrix *new_matrix = new ComplexMatrix (nr-num_to_delete, nc); if (nr > num_to_delete) { int ii = 0; int idx = 0; for (int i = 0; i < nr; i++) { double itmp = ib + idx * iinc; int row = tree_to_mat_idx (itmp); if (i == row && row <= max_idx) idx++; else { for (int j = 0; j < nc; j++) new_matrix->elem (ii, j) = complex_matrix->elem (i, j); ii++; } } } delete complex_matrix; complex_matrix = new_matrix; } else panic_impossible (); } void TC_REP::delete_column (int idx) { if (type_tag == matrix_constant) { int nr = matrix->rows (); int nc = matrix->columns (); Matrix *new_matrix = new Matrix (nr, nc-1); int jj = 0; for (int j = 0; j < nc; j++) { if (j != idx) { for (int i = 0; i < nr; i++) new_matrix->elem (i, jj) = matrix->elem (i, j); jj++; } } delete matrix; matrix = new_matrix; } else if (type_tag == complex_matrix_constant) { int nr = complex_matrix->rows (); int nc = complex_matrix->columns (); ComplexMatrix *new_matrix = new ComplexMatrix (nr, nc-1); int jj = 0; for (int j = 0; j < nc; j++) { if (j != idx) { for (int i = 0; i < nr; i++) new_matrix->elem (i, jj) = complex_matrix->elem (i, j); jj++; } } delete complex_matrix; complex_matrix = new_matrix; } else panic_impossible (); } void TC_REP::delete_columns (idx_vector& jv) { jv.sort_uniq (); int num_to_delete = jv.length (); int nr = rows (); int nc = columns (); // If deleting all columns of a row vector, make result 0x0. if (nr == 1 && num_to_delete == nc) nr = 0; if (type_tag == matrix_constant) { Matrix *new_matrix = new Matrix (nr, nc-num_to_delete); if (nc > num_to_delete) { int jj = 0; int idx = 0; for (int j = 0; j < nc; j++) { if (j == jv.elem (idx)) idx++; else { for (int i = 0; i < nr; i++) new_matrix->elem (i, jj) = matrix->elem (i, j); jj++; } } } delete matrix; matrix = new_matrix; } else if (type_tag == complex_matrix_constant) { ComplexMatrix *new_matrix = new ComplexMatrix (nr, nc-num_to_delete); if (nc > num_to_delete) { int jj = 0; int idx = 0; for (int j = 0; j < nc; j++) { if (j == jv.elem (idx)) idx++; else { for (int i = 0; i < nr; i++) new_matrix->elem (i, jj) = complex_matrix->elem (i, j); jj++; } } } delete complex_matrix; complex_matrix = new_matrix; } else panic_impossible (); } void TC_REP::delete_columns (Range& rj) { rj.sort (); int num_to_delete = rj.nelem (); int nr = rows (); int nc = columns (); // If deleting all columns of a row vector, make result 0x0. if (nr == 1 && num_to_delete == nc) nr = 0; double jb = rj.base (); double jinc = rj.inc (); int max_idx = tree_to_mat_idx (rj.max ()); if (type_tag == matrix_constant) { Matrix *new_matrix = new Matrix (nr, nc-num_to_delete); if (nc > num_to_delete) { int jj = 0; int idx = 0; for (int j = 0; j < nc; j++) { double jtmp = jb + idx * jinc; int col = tree_to_mat_idx (jtmp); if (j == col && col <= max_idx) idx++; else { for (int i = 0; i < nr; i++) new_matrix->elem (i, jj) = matrix->elem (i, j); jj++; } } } delete matrix; matrix = new_matrix; } else if (type_tag == complex_matrix_constant) { ComplexMatrix *new_matrix = new ComplexMatrix (nr, nc-num_to_delete); if (nc > num_to_delete) { int jj = 0; int idx = 0; for (int j = 0; j < nc; j++) { double jtmp = jb + idx * jinc; int col = tree_to_mat_idx (jtmp); if (j == col && col <= max_idx) idx++; else { for (int i = 0; i < nr; i++) new_matrix->elem (i, jj) = complex_matrix->elem (i, j); jj++; } } } delete complex_matrix; complex_matrix = new_matrix; } else panic_impossible (); } void TC_REP::maybe_mutate (void) { if (error_state) return; switch (type_tag) { case complex_scalar_constant: if (::imag (*complex_scalar) == 0.0) { double d = ::real (*complex_scalar); delete complex_scalar; scalar = d; type_tag = scalar_constant; } break; case complex_matrix_constant: if (! any_element_is_complex (*complex_matrix)) { Matrix *m = new Matrix (::real (*complex_matrix)); delete complex_matrix; matrix = m; type_tag = matrix_constant; } break; case scalar_constant: case matrix_constant: case string_constant: case range_constant: case magic_colon: break; default: panic_impossible (); break; } // Avoid calling rows() and columns() for things like magic_colon. int nr = 1; int nc = 1; if (type_tag == matrix_constant || type_tag == complex_matrix_constant || type_tag == range_constant) { nr = rows (); nc = columns (); } switch (type_tag) { case matrix_constant: if (nr == 1 && nc == 1) { double d = matrix->elem (0, 0); delete matrix; scalar = d; type_tag = scalar_constant; } break; case complex_matrix_constant: if (nr == 1 && nc == 1) { Complex c = complex_matrix->elem (0, 0); delete complex_matrix; complex_scalar = new Complex (c); type_tag = complex_scalar_constant; } break; case range_constant: if (nr == 1 && nc == 1) { double d = range->base (); delete range; scalar = d; type_tag = scalar_constant; } break; default: break; } } void TC_REP::print (void) { if (error_state) return; if (print) { ostrstream output_buf; switch (type_tag) { case scalar_constant: octave_print_internal (output_buf, scalar); break; case matrix_constant: octave_print_internal (output_buf, *matrix); break; case complex_scalar_constant: octave_print_internal (output_buf, *complex_scalar); break; case complex_matrix_constant: octave_print_internal (output_buf, *complex_matrix); break; case string_constant: output_buf << string << "\n"; break; case range_constant: octave_print_internal (output_buf, *range); break; case magic_colon: default: panic_impossible (); break; } output_buf << ends; maybe_page_output (output_buf); } } static char * undo_string_escapes (char c) { static char retval[2]; retval[1] = '\0'; if (! c) return 0; switch (c) { case '\a': return "\\a"; case '\b': // backspace return "\\b"; case '\f': // formfeed return "\\f"; case '\n': // newline return "\\n"; case '\r': // carriage return return "\\r"; case '\t': // horizontal tab return "\\t"; case '\v': // vertical tab return "\\v"; case '\\': // backslash return "\\\\"; case '"': // double quote return "\\\""; default: retval[0] = c; return retval; } } void TC_REP::print_code (ostream& os) { switch (type_tag) { case scalar_constant: if (orig_text) os << orig_text; else octave_print_internal (os, scalar, 1); break; case matrix_constant: octave_print_internal (os, *matrix, 1); break; case complex_scalar_constant: { double re = complex_scalar->real (); double im = complex_scalar->imag (); // If we have the original text and a pure imaginary, just print the // original text, because this must be a constant that was parsed as // part of a function. if (orig_text && re == 0.0 && im > 0.0) os << orig_text; else octave_print_internal (os, *complex_scalar, 1); } break; case complex_matrix_constant: octave_print_internal (os, *complex_matrix, 1); break; case string_constant: { os << "\""; char *s, *t = string; while (s = undo_string_escapes (*t++)) os << s; os << "\""; } break; case range_constant: octave_print_internal (os, *range, 1); break; case magic_colon: os << ":"; break; default: panic_impossible (); break; } } char * TC_REP::type_as_string (void) const { switch (type_tag) { case scalar_constant: return "real scalar"; case matrix_constant: return "real matrix"; case complex_scalar_constant: return "complex scalar"; case complex_matrix_constant: return "complex matrix"; case string_constant: return "string"; case range_constant: return "range"; default: return "<unknown type>"; } } tree_constant do_binary_op (tree_constant& a, tree_constant& b, tree_expression::type t) { tree_constant ans; int first_empty = (a.rows () == 0 || a.columns () == 0); int second_empty = (b.rows () == 0 || b.columns () == 0); if (first_empty || second_empty) { int flag = user_pref.propagate_empty_matrices; if (flag < 0) warning ("binary operation on empty matrix"); else if (flag == 0) { ::error ("invalid binary operation on empty matrix"); return ans; } } tree_constant tmp_a = a.make_numeric (); tree_constant tmp_b = b.make_numeric (); TC_REP::constant_type a_type = tmp_a.const_type (); TC_REP::constant_type b_type = tmp_b.const_type (); double d1, d2; Matrix m1, m2; Complex c1, c2; ComplexMatrix cm1, cm2; switch (a_type) { case TC_REP::scalar_constant: d1 = tmp_a.double_value (); switch (b_type) { case TC_REP::scalar_constant: d2 = tmp_b.double_value (); ans = do_binary_op (d1, d2, t); break; case TC_REP::matrix_constant: m2 = tmp_b.matrix_value (); ans = do_binary_op (d1, m2, t); break; case TC_REP::complex_scalar_constant: c2 = tmp_b.complex_value (); ans = do_binary_op (d1, c2, t); break; case TC_REP::complex_matrix_constant: cm2 = tmp_b.complex_matrix_value (); ans = do_binary_op (d1, cm2, t); break; case TC_REP::magic_colon: default: panic_impossible (); break; } break; case TC_REP::matrix_constant: m1 = tmp_a.matrix_value (); switch (b_type) { case TC_REP::scalar_constant: d2 = tmp_b.double_value (); ans = do_binary_op (m1, d2, t); break; case TC_REP::matrix_constant: m2 = tmp_b.matrix_value (); ans = do_binary_op (m1, m2, t); break; case TC_REP::complex_scalar_constant: c2 = tmp_b.complex_value (); ans = do_binary_op (m1, c2, t); break; case TC_REP::complex_matrix_constant: cm2 = tmp_b.complex_matrix_value (); ans = do_binary_op (m1, cm2, t); break; case TC_REP::magic_colon: default: panic_impossible (); break; } break; case TC_REP::complex_scalar_constant: c1 = tmp_a.complex_value (); switch (b_type) { case TC_REP::scalar_constant: d2 = tmp_b.double_value (); ans = do_binary_op (c1, d2, t); break; case TC_REP::matrix_constant: m2 = tmp_b.matrix_value (); ans = do_binary_op (c1, m2, t); break; case TC_REP::complex_scalar_constant: c2 = tmp_b.complex_value (); ans = do_binary_op (c1, c2, t); break; case TC_REP::complex_matrix_constant: cm2 = tmp_b.complex_matrix_value (); ans = do_binary_op (c1, cm2, t); break; case TC_REP::magic_colon: default: panic_impossible (); break; } break; case TC_REP::complex_matrix_constant: cm1 = tmp_a.complex_matrix_value (); switch (b_type) { case TC_REP::scalar_constant: d2 = tmp_b.double_value (); ans = do_binary_op (cm1, d2, t); break; case TC_REP::matrix_constant: m2 = tmp_b.matrix_value (); ans = do_binary_op (cm1, m2, t); break; case TC_REP::complex_scalar_constant: c2 = tmp_b.complex_value (); ans = do_binary_op (cm1, c2, t); break; case TC_REP::complex_matrix_constant: cm2 = tmp_b.complex_matrix_value (); ans = do_binary_op (cm1, cm2, t); break; case TC_REP::magic_colon: default: panic_impossible (); break; } break; case TC_REP::magic_colon: default: panic_impossible (); break; } return ans; } tree_constant do_unary_op (tree_constant& a, tree_expression::type t) { tree_constant ans; if (a.rows () == 0 || a.columns () == 0) { int flag = user_pref.propagate_empty_matrices; if (flag < 0) warning ("unary operation on empty matrix"); else if (flag == 0) { ::error ("invalid unary operation on empty matrix"); return ans; } } tree_constant tmp_a = a.make_numeric (); switch (tmp_a.const_type ()) { case TC_REP::scalar_constant: ans = do_unary_op (tmp_a.double_value (), t); break; case TC_REP::matrix_constant: { Matrix m = tmp_a.matrix_value (); ans = do_unary_op (m, t); } break; case TC_REP::complex_scalar_constant: ans = do_unary_op (tmp_a.complex_value (), t); break; case TC_REP::complex_matrix_constant: { ComplexMatrix m = tmp_a.complex_matrix_value (); ans = do_unary_op (m, t); } break; case TC_REP::magic_colon: default: panic_impossible (); break; } return ans; } tree_constant TC_REP::cumprod (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.cumprod (); } tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar); break; case matrix_constant: { Matrix m = matrix->cumprod (); retval = tree_constant (m); } break; case complex_scalar_constant: retval = tree_constant (*complex_scalar); break; case complex_matrix_constant: { ComplexMatrix m = complex_matrix->cumprod (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } tree_constant TC_REP::cumsum (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.cumsum (); } tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar); break; case matrix_constant: { Matrix m = matrix->cumsum (); retval = tree_constant (m); } break; case complex_scalar_constant: retval = tree_constant (*complex_scalar); break; case complex_matrix_constant: { ComplexMatrix m = complex_matrix->cumsum (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } tree_constant TC_REP::prod (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.prod (); } tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar); break; case matrix_constant: { Matrix m = matrix->prod (); retval = tree_constant (m); } break; case complex_scalar_constant: retval = tree_constant (*complex_scalar); break; case complex_matrix_constant: { ComplexMatrix m = complex_matrix->prod (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } tree_constant TC_REP::sum (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.sum (); } tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar); break; case matrix_constant: { Matrix m = matrix->sum (); retval = tree_constant (m); } break; case complex_scalar_constant: retval = tree_constant (*complex_scalar); break; case complex_matrix_constant: { ComplexMatrix m = complex_matrix->sum (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } tree_constant TC_REP::sumsq (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.sumsq (); } tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar * scalar); break; case matrix_constant: { Matrix m = matrix->sumsq (); retval = tree_constant (m); } break; case complex_scalar_constant: { Complex c (*complex_scalar); retval = tree_constant (c * c); } break; case complex_matrix_constant: { ComplexMatrix m = complex_matrix->sumsq (); retval = tree_constant (m); } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } static tree_constant make_diag (const Matrix& v, int k) { int nr = v.rows (); int nc = v.columns (); assert (nc == 1 || nr == 1); tree_constant retval; int roff = 0; int coff = 0; if (k > 0) { roff = 0; coff = k; } else if (k < 0) { roff = -k; coff = 0; } if (nr == 1) { int n = nc + ABS (k); Matrix m (n, n, 0.0); for (int i = 0; i < nc; i++) m.elem (i+roff, i+coff) = v.elem (0, i); retval = tree_constant (m); } else { int n = nr + ABS (k); Matrix m (n, n, 0.0); for (int i = 0; i < nr; i++) m.elem (i+roff, i+coff) = v.elem (i, 0); retval = tree_constant (m); } return retval; } static tree_constant make_diag (const ComplexMatrix& v, int k) { int nr = v.rows (); int nc = v.columns (); assert (nc == 1 || nr == 1); tree_constant retval; int roff = 0; int coff = 0; if (k > 0) { roff = 0; coff = k; } else if (k < 0) { roff = -k; coff = 0; } if (nr == 1) { int n = nc + ABS (k); ComplexMatrix m (n, n, 0.0); for (int i = 0; i < nc; i++) m.elem (i+roff, i+coff) = v.elem (0, i); retval = tree_constant (m); } else { int n = nr + ABS (k); ComplexMatrix m (n, n, 0.0); for (int i = 0; i < nr; i++) m.elem (i+roff, i+coff) = v.elem (i, 0); retval = tree_constant (m); } return retval; } tree_constant TC_REP::diag (void) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.diag (); } tree_constant retval; switch (type_tag) { case scalar_constant: retval = tree_constant (scalar); break; case matrix_constant: { int nr = rows (); int nc = columns (); if (nr == 0 || nc == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else if (nr == 1 || nc == 1) retval = make_diag (matrix_value (), 0); else { ColumnVector v = matrix->diag (); if (v.capacity () > 0) retval = tree_constant (v); } } break; case complex_scalar_constant: retval = tree_constant (*complex_scalar); break; case complex_matrix_constant: { int nr = rows (); int nc = columns (); if (nr == 0 || nc == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else if (nr == 1 || nc == 1) retval = make_diag (complex_matrix_value (), 0); else { ComplexColumnVector v = complex_matrix->diag (); if (v.capacity () > 0) retval = tree_constant (v); } } break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } tree_constant TC_REP::diag (const tree_constant& a) const { if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.diag (a); } tree_constant tmp_a = a.make_numeric (); TC_REP::constant_type a_type = tmp_a.const_type (); tree_constant retval; switch (type_tag) { case scalar_constant: if (a_type == scalar_constant) { int k = NINT (tmp_a.double_value ()); int n = ABS (k) + 1; if (k == 0) retval = tree_constant (scalar); else if (k > 0) { Matrix m (n, n, 0.0); m.elem (0, k) = scalar; retval = tree_constant (m); } else if (k < 0) { Matrix m (n, n, 0.0); m.elem (-k, 0) = scalar; retval = tree_constant (m); } } break; case matrix_constant: if (a_type == scalar_constant) { int k = NINT (tmp_a.double_value ()); int nr = rows (); int nc = columns (); if (nr == 0 || nc == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else if (nr == 1 || nc == 1) retval = make_diag (matrix_value (), k); else { ColumnVector d = matrix->diag (k); retval = tree_constant (d); } } else ::error ("diag: invalid second argument"); break; case complex_scalar_constant: if (a_type == scalar_constant) { int k = NINT (tmp_a.double_value ()); int n = ABS (k) + 1; if (k == 0) retval = tree_constant (*complex_scalar); else if (k > 0) { ComplexMatrix m (n, n, 0.0); m.elem (0, k) = *complex_scalar; retval = tree_constant (m); } else if (k < 0) { ComplexMatrix m (n, n, 0.0); m.elem (-k, 0) = *complex_scalar; retval = tree_constant (m); } } break; case complex_matrix_constant: if (a_type == scalar_constant) { int k = NINT (tmp_a.double_value ()); int nr = rows (); int nc = columns (); if (nr == 0 || nc == 0) { Matrix mtmp; retval = tree_constant (mtmp); } else if (nr == 1 || nc == 1) retval = make_diag (complex_matrix_value (), k); else { ComplexColumnVector d = complex_matrix->diag (k); retval = tree_constant (d); } } else ::error ("diag: invalid second argument"); break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } // XXX FIXME XXX -- this can probably be rewritten efficiently as a // nonmember function... tree_constant TC_REP::mapper (Mapper_fcn& m_fcn, int print) const { tree_constant retval; if (type_tag == string_constant || type_tag == range_constant) { tree_constant tmp = make_numeric (); return tmp.mapper (m_fcn, print); } switch (type_tag) { case scalar_constant: if (m_fcn.can_return_complex_for_real_arg && (scalar < m_fcn.lower_limit || scalar > m_fcn.upper_limit)) { if (m_fcn.c_c_mapper) { Complex c = m_fcn.c_c_mapper (Complex (scalar)); retval = tree_constant (c); } else ::error ("%s: unable to handle real arguments", m_fcn.name); } else { if (m_fcn.d_d_mapper) { double d = m_fcn.d_d_mapper (scalar); retval = tree_constant (d); } else ::error ("%s: unable to handle real arguments", m_fcn.name); } break; case matrix_constant: if (m_fcn.can_return_complex_for_real_arg && (any_element_less_than (*matrix, m_fcn.lower_limit) || any_element_greater_than (*matrix, m_fcn.upper_limit))) { if (m_fcn.c_c_mapper) { ComplexMatrix cm = map (m_fcn.c_c_mapper, ComplexMatrix (*matrix)); retval = tree_constant (cm); } else ::error ("%s: unable to handle real arguments", m_fcn.name); } else { if (m_fcn.d_d_mapper) { Matrix m = map (m_fcn.d_d_mapper, *matrix); retval = tree_constant (m); } else ::error ("%s: unable to handle real arguments", m_fcn.name); } break; case complex_scalar_constant: if (m_fcn.d_c_mapper) { double d; d = m_fcn.d_c_mapper (*complex_scalar); retval = tree_constant (d); } else if (m_fcn.c_c_mapper) { Complex c; c = m_fcn.c_c_mapper (*complex_scalar); retval = tree_constant (c); } else ::error ("%s: unable to handle complex arguments", m_fcn.name); break; case complex_matrix_constant: if (m_fcn.d_c_mapper) { Matrix m; m = map (m_fcn.d_c_mapper, *complex_matrix); retval = tree_constant (m); } else if (m_fcn.c_c_mapper) { ComplexMatrix cm; cm = map (m_fcn.c_c_mapper, *complex_matrix); retval = tree_constant (cm); } else ::error ("%s: unable to handle complex arguments", m_fcn.name); break; case string_constant: case range_constant: case magic_colon: default: panic_impossible (); break; } return retval; } /* ;;; Local Variables: *** ;;; mode: C++ *** ;;; page-delimiter: "^/\\*" *** ;;; End: *** */