# HG changeset patch # User jwe # Date 770863270 0 # Node ID 38cb8809591339357d086ffed40838307d4b7ae9 # Parent 3d4b4f0fa5ba3b6396dfb0b4500fbb02e7234cca [project @ 1994-06-06 00:41:10 by jwe] Initial revision diff --git a/liboctave/CColVector.cc b/liboctave/CColVector.cc new file mode 100644 --- /dev/null +++ b/liboctave/CColVector.cc @@ -0,0 +1,659 @@ +// ColumnVector manipulations. -*- 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 + +#include + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "f77-uscore.h" +#include "lo-error.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (zgemm) (const char*, const char*, const int*, + const int*, const int*, const Complex*, + const Complex*, const int*, const Complex*, + const int*, const Complex*, Complex*, const int*, + long, long); +} + +/* + * Complex Column Vector class + */ + +#define KLUDGE_VECTORS +#define TYPE Complex +#define KL_VEC_TYPE ComplexColumnVector +#include "mx-kludge.cc" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +ComplexColumnVector::ComplexColumnVector (const ColumnVector& a) + : Array (a.length ()) +{ + for (int i = 0; i < length (); i++) + elem (i) = a.elem (i); +} + +#if 0 +ComplexColumnVector& +ComplexColumnVector::resize (int n) +{ + if (n < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimension"); + return *this; + } + + Complex *new_data = (Complex *) NULL; + if (n > 0) + { + new_data = new Complex [n]; + int min_len = len < n ? len : n; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + len = n; + data = new_data; + + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::resize (int n, double val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::resize (int n, const Complex& val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} +#endif + +int +ComplexColumnVector::operator == (const ComplexColumnVector& a) const +{ + int len = length (); + if (len != a.length ()) + return 0; + return equal (data (), a.data (), len); +} + +int +ComplexColumnVector::operator != (const ComplexColumnVector& a) const +{ + return !(*this == a); +} + +// destructive insert/delete/reorder operations + +ComplexColumnVector& +ComplexColumnVector::insert (const ColumnVector& a, int r) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); + + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::insert (const ComplexColumnVector& a, int r) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); + + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::fill (double val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::fill (const Complex& val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::fill (double val, int r1, int r2) +{ + int len = length (); + if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + for (int i = r1; i <= r2; i++) + elem (i) = val; + + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::fill (const Complex& val, int r1, int r2) +{ + int len = length (); + if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + for (int i = r1; i <= r2; i++) + elem (i) = val; + + return *this; +} + +ComplexColumnVector +ComplexColumnVector::stack (const ColumnVector& a) const +{ + int len = length (); + int nr_insert = len; + ComplexColumnVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nr_insert); + return retval; +} + +ComplexColumnVector +ComplexColumnVector::stack (const ComplexColumnVector& a) const +{ + int len = length (); + int nr_insert = len; + ComplexColumnVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nr_insert); + return retval; +} + +ComplexRowVector +ComplexColumnVector::hermitian (void) const +{ + int len = length (); + return ComplexRowVector (conj_dup (data (), len), len); +} + +ComplexRowVector +ComplexColumnVector::transpose (void) const +{ + int len = length (); + return ComplexRowVector (dup (data (), len), len); +} + +ColumnVector +real (const ComplexColumnVector& a) +{ + int a_len = a.length (); + ColumnVector retval; + if (a_len > 0) + retval = ColumnVector (real_dup (a.data (), a_len), a_len); + return retval; +} + +ColumnVector +imag (const ComplexColumnVector& a) +{ + int a_len = a.length (); + ColumnVector retval; + if (a_len > 0) + retval = ColumnVector (imag_dup (a.data (), a_len), a_len); + return retval; +} + +ComplexColumnVector +conj (const ComplexColumnVector& a) +{ + int a_len = a.length (); + ComplexColumnVector retval; + if (a_len > 0) + retval = ComplexColumnVector (conj_dup (a.data (), a_len), a_len); + return retval; +} + +// resize is the destructive equivalent for this one + +ComplexColumnVector +ComplexColumnVector::extract (int r1, int r2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + int new_r = r2 - r1 + 1; + + ComplexColumnVector result (new_r); + + for (int i = 0; i < new_r; i++) + result.elem (i) = elem (r1+i); + + return result; +} + +// column vector by column vector -> column vector operations + +ComplexColumnVector& +ComplexColumnVector::operator += (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::operator -= (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::operator += (const ComplexColumnVector& a) +{ + int len = length (); + + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::operator -= (const ComplexColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +// column vector by scalar -> column vector operations + +ComplexColumnVector +operator + (const ComplexColumnVector& v, double s) +{ + int len = v.length (); + return ComplexColumnVector (add (v.data (), len, s), len); +} + +ComplexColumnVector +operator - (const ComplexColumnVector& v, double s) +{ + int len = v.length (); + return ComplexColumnVector (subtract (v.data (), len, s), len); +} + +ComplexColumnVector +operator * (const ComplexColumnVector& v, double s) +{ + int len = v.length (); + return ComplexColumnVector (multiply (v.data (), len, s), len); +} + +ComplexColumnVector +operator / (const ComplexColumnVector& v, double s) +{ + int len = v.length (); + return ComplexColumnVector (divide (v.data (), len, s), len); +} + +// scalar by column vector -> column vector operations + +ComplexColumnVector +operator + (double s, const ComplexColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (add (a.data (), a_len, s), a_len); +} + +ComplexColumnVector +operator - (double s, const ComplexColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (subtract (s, a.data (), a_len), a_len); +} + +ComplexColumnVector +operator * (double s, const ComplexColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (multiply (a.data (), a_len, s), a_len); +} + +ComplexColumnVector +operator / (double s, const ComplexColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (divide (s, a.data (), a_len), a_len); +} + +// column vector by row vector -> matrix operations + +ComplexMatrix +operator * (const ComplexColumnVector& v, const ComplexRowVector& a) +{ + int len = v.length (); + int a_len = a.length (); + if (len != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return ComplexMatrix (); + } + + if (len == 0) + return ComplexMatrix (len, len, 0.0); + + char transa = 'N'; + char transb = 'N'; + Complex alpha (1.0); + Complex beta (0.0); + int anr = 1; + + Complex *c = new Complex [len * a_len]; + + F77_FCN (zgemm) (&transa, &transb, &len, &a_len, &anr, &alpha, + v.data (), &len, a.data (), &anr, &beta, c, &len, + 1L, 1L); + + return ComplexMatrix (c, len, a_len); +} + +// column vector by column vector -> column vector operations + +ComplexColumnVector +operator + (const ComplexColumnVector& v, const ColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector addition attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (add (v.data (), a.data (), len), len); +} + +ComplexColumnVector +operator - (const ComplexColumnVector& v, const ColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (subtract (v.data (), a.data (), len), len); +} + +ComplexColumnVector +product (const ComplexColumnVector& v, const ColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector product attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (multiply (v.data (), a.data (), len), len); +} + +ComplexColumnVector +quotient (const ComplexColumnVector& v, const ColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector quotient attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (divide (v.data (), a.data (), len), len); +} + +// other operations + +ComplexColumnVector +map (c_c_Mapper f, const ComplexColumnVector& a) +{ + ComplexColumnVector b (a); + b.map (f); + return b; +} + +ColumnVector +map (d_c_Mapper f, const ComplexColumnVector& a) +{ + int a_len = a.length (); + ColumnVector b (a_len); + for (int i = 0; i < a_len; i++) + b.elem (i) = f (a.elem (i)); + return b; +} + +void +ComplexColumnVector::map (c_c_Mapper f) +{ + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); +} + +Complex +ComplexColumnVector::min (void) const +{ + int len = length (); + if (len == 0) + return 0.0; + + Complex res = elem (0); + double absres = abs (res); + + for (int i = 1; i < len; i++) + if (abs (elem (i)) < absres) + { + res = elem (i); + absres = abs (res); + } + + return res; +} + +Complex +ComplexColumnVector::max (void) const +{ + int len = length (); + if (len == 0) + return 0.0; + + Complex res = elem (0); + double absres = abs (res); + + for (int i = 1; i < len; i++) + if (abs (elem (i)) > absres) + { + res = elem (i); + absres = abs (res); + } + + return res; +} + +// i/o + +ostream& +operator << (ostream& os, const ComplexColumnVector& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.length (); i++) + os << /* setw (field_width) << */ a.elem (i) << "\n"; + return os; +} + +istream& +operator >> (istream& is, ComplexColumnVector& a) +{ + int len = a.length(); + + if (len < 1) + is.clear (ios::badbit); + else + { + double tmp; + for (int i = 0; i < len; i++) + { + is >> tmp; + if (is) + a.elem (i) = tmp; + else + break; + } + } +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CColVector.h b/liboctave/CColVector.h new file mode 100644 --- /dev/null +++ b/liboctave/CColVector.h @@ -0,0 +1,171 @@ +// -*- 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. + +*/ + +#if !defined (octave_ComplexColumnVector_h) +#define octave_ComplexColumnVector_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "mx-defs.h" + +extern "C++" { + +class ComplexColumnVector : public Array +{ +friend class ComplexRowVector; + +public: + + ComplexColumnVector (void) : Array () { } + ComplexColumnVector (int n) : Array (n) { } + ComplexColumnVector (int n, const Complex& val) + : Array (n, val) { } + ComplexColumnVector (const ColumnVector& a); + ComplexColumnVector (const Array& a) : Array (a) { } + ComplexColumnVector (const ComplexColumnVector& a) : Array (a) { } +// ComplexColumnVector (double a) : Array (1, a) { } +// ComplexColumnVector (const Complex& a) : Array (1, a) { } + + ComplexColumnVector& operator = (const ComplexColumnVector& a) + { + Array::operator = (a); + return *this; + } + +// operator Array& () const { return *this; } + + int operator == (const ComplexColumnVector& a) const; + int operator != (const ComplexColumnVector& a) const; + +// destructive insert/delete/reorder operations + + ComplexColumnVector& insert (const ColumnVector& a, int r); + ComplexColumnVector& insert (const ComplexColumnVector& a, int r); + + ComplexColumnVector& fill (double val); + ComplexColumnVector& fill (const Complex& val); + ComplexColumnVector& fill (double val, int r1, int r2); + ComplexColumnVector& fill (const Complex& val, int r1, int r2); + + ComplexColumnVector stack (const ColumnVector& a) const; + ComplexColumnVector stack (const ComplexColumnVector& a) const; + + ComplexRowVector hermitian (void) const; // complex conjugate transpose. + ComplexRowVector transpose (void) const; + + friend ColumnVector real (const ComplexColumnVector& a); + friend ColumnVector imag (const ComplexColumnVector& a); + friend ComplexColumnVector conj (const ComplexColumnVector& a); + +// resize is the destructive equivalent for this one + + ComplexColumnVector extract (int r1, int r2) const; + +// column vector by column vector -> column vector operations + + ComplexColumnVector& operator += (const ColumnVector& a); + ComplexColumnVector& operator -= (const ColumnVector& a); + + ComplexColumnVector& operator += (const ComplexColumnVector& a); + ComplexColumnVector& operator -= (const ComplexColumnVector& a); + +// column vector by scalar -> column vector operations + + friend ComplexColumnVector operator + (const ComplexColumnVector& a, + double s); + friend ComplexColumnVector operator - (const ComplexColumnVector& a, + double s); + friend ComplexColumnVector operator * (const ComplexColumnVector& a, + double s); + friend ComplexColumnVector operator / (const ComplexColumnVector& a, + double s); + +// scalar by column vector -> column vector operations + + friend ComplexColumnVector operator + (double s, + const ComplexColumnVector& a); + friend ComplexColumnVector operator - (double s, + const ComplexColumnVector& a); + friend ComplexColumnVector operator * (double s, + const ComplexColumnVector& a); + friend ComplexColumnVector operator / (double s, + const ComplexColumnVector& a); + +// column vector by row vector -> matrix operations + + friend ComplexMatrix operator * (const ComplexColumnVector& a, + const ComplexRowVector& b); + +// column vector by column vector -> column vector operations + + friend ComplexColumnVector operator + (const ComplexColumnVector& a, + const ColumnVector& b); + friend ComplexColumnVector operator - (const ComplexColumnVector& a, + const ColumnVector& b); + + friend ComplexColumnVector product (const ComplexColumnVector& a, + const ColumnVector& b); + friend ComplexColumnVector quotient (const ComplexColumnVector& a, + const ColumnVector& b); + +// other operations + + friend ComplexColumnVector map (c_c_Mapper f, const ComplexColumnVector& a); + friend ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a); + void map (c_c_Mapper f); + + Complex min (void) const; + Complex max (void) const; + +// i/o + + friend ostream& operator << (ostream& os, const ComplexColumnVector& a); + friend ostream& operator >> (ostream& is, ComplexColumnVector& a); + +#define KLUDGE_VECTORS +#define TYPE Complex +#define KL_VEC_TYPE ComplexColumnVector +#include "mx-kludge.h" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +private: + + ComplexColumnVector (Complex *d, int l) : Array (d, l) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CDiagMatrix.cc b/liboctave/CDiagMatrix.cc new file mode 100644 --- /dev/null +++ b/liboctave/CDiagMatrix.cc @@ -0,0 +1,1168 @@ +// DiagMatrix manipulations. -*- 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 + +#include + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "lo-error.h" + +/* + * Complex Diagonal Matrix class + */ + +#define KLUDGE_DIAG_MATRICES +#define TYPE Complex +#define KL_DMAT_TYPE ComplexDiagMatrix +#include "mx-kludge.cc" +#undef KLUDGE_DIAG_MATRICES +#undef TYPE +#undef KL_DMAT_TYPE + +ComplexDiagMatrix::ComplexDiagMatrix (const RowVector& a) + : DiagArray (a.length ()) +{ + for (int i = 0; i < length (); i++) + elem (i, i) = a.elem (i); +} + +ComplexDiagMatrix::ComplexDiagMatrix (const ColumnVector& a) + : DiagArray (a.length ()) +{ + for (int i = 0; i < length (); i++) + elem (i, i) = a.elem (i); +} + +ComplexDiagMatrix::ComplexDiagMatrix (const DiagMatrix& a) + : DiagArray (a.rows (), a.cols ()) +{ + for (int i = 0; i < length (); i++) + elem (i, i) = a.elem (i, i); +} + +#if 0 +ComplexDiagMatrix& +ComplexDiagMatrix::resize (int r, int c) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r < c ? r : c; + Complex *new_data = (Complex *) NULL; + if (new_len > 0) + { + new_data = new Complex [new_len]; + + int min_len = new_len < len ? new_len : len; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::resize (int r, int c, double val) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r < c ? r : c; + Complex *new_data = (Complex *) NULL; + if (new_len > 0) + { + new_data = new Complex [new_len]; + + int min_len = new_len < len ? new_len : len; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + + for (i = min_len; i < new_len; i++) + new_data[i] = val; + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::resize (int r, int c, const Complex& val) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r < c ? r : c; + Complex *new_data = (Complex *) NULL; + if (new_len > 0) + { + new_data = new Complex [new_len]; + + int min_len = new_len < len ? new_len : len; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + + for (i = min_len; i < new_len; i++) + new_data[i] = val; + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} +#endif + +int +ComplexDiagMatrix::operator == (const ComplexDiagMatrix& a) const +{ + if (rows () != a.rows () || cols () != a.cols ()) + return 0; + + return equal (data (), a.data (), length ()); +} + +int +ComplexDiagMatrix::operator != (const ComplexDiagMatrix& a) const +{ + return !(*this == a); +} + +ComplexDiagMatrix +ComplexDiagMatrix::hermitian (void) const +{ + return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ()); +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (double val) +{ + for (int i = 0; i < length (); i++) + elem (i, i) = val; + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const Complex& val) +{ + for (int i = 0; i < length (); i++) + elem (i, i) = val; + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (double val, int beg, int end) +{ + if (beg < 0 || end >= length () || end < beg) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = beg; i < end; i++) + elem (i, i) = val; + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const Complex& val, int beg, int end) +{ + if (beg < 0 || end >= length () || end < beg) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = beg; i < end; i++) + elem (i, i) = val; + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const ColumnVector& a) +{ + int len = length (); + if (a.length () != len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const ComplexColumnVector& a) +{ + int len = length (); + if (a.length () != len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const RowVector& a) +{ + int len = length (); + if (a.length () != len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const ComplexRowVector& a) +{ + int len = length (); + if (a.length () != len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const ColumnVector& a, int beg) +{ + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const ComplexColumnVector& a, int beg) +{ + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const RowVector& a, int beg) +{ + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::fill (const ComplexRowVector& a, int beg) +{ + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + + return *this; +} + +ComplexDiagMatrix +ComplexDiagMatrix::transpose (void) const +{ + return ComplexDiagMatrix (dup (data (), length ()), cols (), rows ()); +} + +DiagMatrix +real (const ComplexDiagMatrix& a) +{ + DiagMatrix retval; + int a_len = a.length (); + if (a_len > 0) + retval = DiagMatrix (real_dup (a.data (), a_len), a.rows (), + a.cols ()); + return retval; +} + +DiagMatrix +imag (const ComplexDiagMatrix& a) +{ + DiagMatrix retval; + int a_len = a.length (); + if (a_len > 0) + retval = DiagMatrix (imag_dup (a.data (), a_len), a.rows (), + a.cols ()); + return retval; +} + +ComplexDiagMatrix +conj (const ComplexDiagMatrix& a) +{ + ComplexDiagMatrix retval; + int a_len = a.length (); + if (a_len > 0) + retval = ComplexDiagMatrix (conj_dup (a.data (), a_len), + a.rows (), a.cols ()); + return retval; +} + +// resize is the destructive analog for this one + +ComplexMatrix +ComplexDiagMatrix::extract (int r1, int c1, int r2, int c2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + int new_r = r2 - r1 + 1; + int new_c = c2 - c1 + 1; + + ComplexMatrix result (new_r, new_c); + + for (int j = 0; j < new_c; j++) + for (int i = 0; i < new_r; i++) + result.elem (i, j) = elem (r1+i, c1+j); + + return result; +} + +// extract row or column i. + +ComplexRowVector +ComplexDiagMatrix::row (int i) const +{ + int nr = rows (); + int nc = cols (); + if (i < 0 || i >= nr) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } + + ComplexRowVector retval (nc, 0.0); + if (nr <= nc || (nr > nc && i < nc)) + retval.elem (i) = elem (i, i); + + return retval; +} + +ComplexRowVector +ComplexDiagMatrix::row (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return ComplexRowVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return row (0); + else if (c == 'l' || c == 'L') + return row (rows () - 1); + else + { + (*current_liboctave_error_handler) ("invalid row selection"); + return ComplexRowVector (); + } +} + +ComplexColumnVector +ComplexDiagMatrix::column (int i) const +{ + int nr = rows (); + int nc = cols (); + if (i < 0 || i >= nc) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } + + ComplexColumnVector retval (nr, 0.0); + if (nr >= nc || (nr < nc && i < nr)) + retval.elem (i) = elem (i, i); + + return retval; +} + +ComplexColumnVector +ComplexDiagMatrix::column (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return column (0); + else if (c == 'l' || c == 'L') + return column (cols () - 1); + else + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } +} + +ComplexDiagMatrix +ComplexDiagMatrix::inverse (void) const +{ + int info; + return inverse (info); +} + +ComplexDiagMatrix +ComplexDiagMatrix::inverse (int& info) const +{ + int nr = rows (); + int nc = cols (); + if (nr != nc) + { + (*current_liboctave_error_handler) ("inverse requires square matrix"); + return DiagMatrix (); + } + + ComplexDiagMatrix retval (nr, nc); + + info = 0; + for (int i = 0; i < length (); i++) + { + if (elem (i, i) == 0.0) + { + info = -1; + return *this; + } + else + retval.elem (i, i) = 1.0 / elem (i, i); + } + + return *this; +} + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + +ComplexDiagMatrix& +ComplexDiagMatrix::operator += (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::operator -= (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::operator += (const ComplexDiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + return *this; +} + +ComplexDiagMatrix& +ComplexDiagMatrix::operator -= (const ComplexDiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + return *this; +} + +// diagonal matrix by scalar -> matrix operations + +ComplexMatrix +operator + (const ComplexDiagMatrix& a, double s) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return a + tmp; +} + +ComplexMatrix +operator - (const ComplexDiagMatrix& a, double s) +{ + ComplexMatrix tmp (a.rows (), a.cols (), -s); + return a + tmp; +} + +ComplexMatrix +operator + (const ComplexDiagMatrix& a, const Complex& s) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return a + tmp; +} + +ComplexMatrix +operator - (const ComplexDiagMatrix& a, const Complex& s) +{ + ComplexMatrix tmp (a.rows (), a.cols (), -s); + return a + tmp; +} + +// diagonal matrix by scalar -> diagonal matrix operations + +ComplexDiagMatrix +operator * (const ComplexDiagMatrix& a, double s) +{ + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexDiagMatrix +operator / (const ComplexDiagMatrix& a, double s) +{ + return ComplexDiagMatrix (divide (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +// scalar by diagonal matrix -> matrix operations + +ComplexMatrix +operator + (double s, const ComplexDiagMatrix& a) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return tmp + a; +} + +ComplexMatrix +operator - (double s, const ComplexDiagMatrix& a) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return tmp - a; +} + +ComplexMatrix +operator + (const Complex& s, const ComplexDiagMatrix& a) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return tmp + a; +} + +ComplexMatrix +operator - (const Complex& s, const ComplexDiagMatrix& a) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return tmp - a; +} + +// scalar by diagonal matrix -> diagonal matrix operations + +ComplexDiagMatrix +operator * (double s, const ComplexDiagMatrix& a) +{ + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +// diagonal matrix by column vector -> column vector operations + +ComplexColumnVector +operator * (const ComplexDiagMatrix& m, const ColumnVector& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_len = a.length (); + if (nc != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant matrix muliplication attempted"); + return ComplexColumnVector (); + } + + if (nc == 0 || nr == 0) + return ComplexColumnVector (0); + + ComplexColumnVector result (nr); + + for (int i = 0; i < a_len; i++) + result.elem (i) = a.elem (i) * m.elem (i, i); + + for (i = a_len; i < nr; i++) + result.elem (i) = 0.0; + + return result; +} + +ComplexColumnVector +operator * (const ComplexDiagMatrix& m, const ComplexColumnVector& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_len = a.length (); + if (nc != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant matrix muliplication attempted"); + return ComplexColumnVector (); + } + + if (nc == 0 || nr == 0) + return ComplexColumnVector (0); + + ComplexColumnVector result (nr); + + for (int i = 0; i < a_len; i++) + result.elem (i) = a.elem (i) * m.elem (i, i); + + for (i = a_len; i < nr; i++) + result.elem (i) = 0.0; + + return result; +} + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + +ComplexDiagMatrix +operator * (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b) +{ + int nr_a = a.rows (); + int nc_a = a.cols (); + int nr_b = b.rows (); + int nc_b = b.cols (); + if (nc_a != nr_b) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexDiagMatrix (); + } + + if (nr_a == 0 || nc_a == 0 || nc_b == 0) + return ComplexDiagMatrix (nr_a, nc_a, 0.0); + + ComplexDiagMatrix c (nr_a, nc_b); + + int len = nr_a < nc_b ? nr_a : nc_b; + + for (int i = 0; i < len; i++) + { + Complex a_element = a.elem (i, i); + Complex b_element = b.elem (i, i); + + if (a_element == 0.0 || b_element == 0.0) + c.elem (i, i) = 0.0; + else if (a_element == 1.0) + c.elem (i, i) = b_element; + else if (b_element == 1.0) + c.elem (i, i) = a_element; + else + c.elem (i, i) = a_element * b_element; + } + + return c; +} + +ComplexDiagMatrix +operator + (const ComplexDiagMatrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexDiagMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexDiagMatrix (nr, nc); + + return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexDiagMatrix +operator - (const ComplexDiagMatrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexDiagMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexDiagMatrix (nr, nc); + + return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()), + nr, nc); +} + +ComplexDiagMatrix +operator * (const ComplexDiagMatrix& a, const DiagMatrix& b) +{ + int nr_a = a.rows (); + int nc_a = a.cols (); + int nr_b = b.rows (); + int nc_b = b.cols (); + if (nc_a != nr_b) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexDiagMatrix (); + } + + if (nr_a == 0 || nc_a == 0 || nc_b == 0) + return ComplexDiagMatrix (nr_a, nc_a, 0.0); + + ComplexDiagMatrix c (nr_a, nc_b); + + int len = nr_a < nc_b ? nr_a : nc_b; + + for (int i = 0; i < len; i++) + { + Complex a_element = a.elem (i, i); + double b_element = b.elem (i, i); + + if (a_element == 0.0 || b_element == 0.0) + c.elem (i, i) = 0.0; + else if (a_element == 1.0) + c.elem (i, i) = b_element; + else if (b_element == 1.0) + c.elem (i, i) = a_element; + else + c.elem (i, i) = a_element * b_element; + } + + return c; +} + +ComplexDiagMatrix +product (const ComplexDiagMatrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix product attempted"); + return ComplexDiagMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexDiagMatrix (nr, nc); + + return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()), + nr, nc); +} + +// diagonal matrix by matrix -> matrix operations + +ComplexMatrix +operator + (const ComplexDiagMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +ComplexMatrix +operator - (const ComplexDiagMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (-a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +ComplexMatrix +operator * (const ComplexDiagMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, a_nc, 0.0); + + ComplexMatrix c (nr, a_nc); + + for (int i = 0; i < m.length (); i++) + { + if (m.elem (i, i) == 1.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = a.elem (i, j); + } + else if (m.elem (i, i) == 0.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = 0.0; + } + else + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = m.elem (i, i) * a.elem (i, j); + } + } + + if (nr > nc) + { + for (int j = 0; j < a_nc; j++) + for (int i = a_nr; i < nr; i++) + c.elem (i, j) = 0.0; + } + + return c; +} + +ComplexMatrix +operator + (const ComplexDiagMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +ComplexMatrix +operator - (const ComplexDiagMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (-a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +ComplexMatrix +operator * (const ComplexDiagMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, a_nc, 0.0); + + ComplexMatrix c (nr, a_nc); + + for (int i = 0; i < m.length (); i++) + { + if (m.elem (i, i) == 1.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = a.elem (i, j); + } + else if (m.elem (i, i) == 0.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = 0.0; + } + else + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = m.elem (i, i) * a.elem (i, j); + } + } + + if (nr > nc) + { + for (int j = 0; j < a_nc; j++) + for (int i = a_nr; i < nr; i++) + c.elem (i, j) = 0.0; + } + + return c; +} + +// other operations + +ComplexColumnVector +ComplexDiagMatrix::diag (void) const +{ + return diag (0); +} + +// Could be optimized... + +ComplexColumnVector +ComplexDiagMatrix::diag (int k) const +{ + int nnr = rows (); + int nnc = cols (); + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + ComplexColumnVector d; + + if (nnr > 0 && nnc > 0) + { + int ndiag = (nnr < nnc) ? nnr : nnc; + + d.resize (ndiag); + + if (k > 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i+k); + } + else if ( k < 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i-k, i); + } + else + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i); + } + } + else + cerr << "diag: requested diagonal out of range\n"; + + return d; +} + +// i/o + +ostream& +operator << (ostream& os, const ComplexDiagMatrix& a) +{ + Complex ZERO (0.0); +// int field_width = os.precision () + 7; + for (int i = 0; i < a.rows (); i++) + { + for (int j = 0; j < a.cols (); j++) + { + if (i == j) + os << " " /* setw (field_width) */ << a.elem (i, i); + else + os << " " /* setw (field_width) */ << ZERO; + } + os << "\n"; + } + return os; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CDiagMatrix.h b/liboctave/CDiagMatrix.h new file mode 100644 --- /dev/null +++ b/liboctave/CDiagMatrix.h @@ -0,0 +1,218 @@ +// -*- 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. + +*/ + +#if !defined (octave_ComplexDiagMatrix_h) +#define octave_ComplexDiagMatrix_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "dRowVector.h" +#include "CRowVector.h" +#include "dColVector.h" +#include "CColVector.h" + +#include "mx-defs.h" + +extern "C++" { + +class ComplexDiagMatrix : public DiagArray +{ +public: + + ComplexDiagMatrix (void) : DiagArray () { } + ComplexDiagMatrix (int n) : DiagArray (n) { } + ComplexDiagMatrix (int n, const Complex& val) + : DiagArray (n, val) { } + ComplexDiagMatrix (int r, int c) : DiagArray (r, c) { } + ComplexDiagMatrix (int r, int c, const Complex& val) + : DiagArray (r, c, val) { } + ComplexDiagMatrix (const RowVector& a); + ComplexDiagMatrix (const ComplexRowVector& a) : DiagArray (a) { } + ComplexDiagMatrix (const ColumnVector& a); + ComplexDiagMatrix (const ComplexColumnVector& a) + : DiagArray (a) { } + ComplexDiagMatrix (const DiagMatrix& a); + ComplexDiagMatrix (const DiagArray& a) + : DiagArray (a) { } + ComplexDiagMatrix (const ComplexDiagMatrix& a) : DiagArray (a) { } +// ComplexDiagMatrix (const Complex& a) : DiagArray (1, a) { } + + ComplexDiagMatrix& operator = (const ComplexDiagMatrix& a) + { + DiagArray::operator = (a); + return *this; + } + +// operator DiagArray& () const { return *this; } + + int operator == (const ComplexDiagMatrix& a) const; + int operator != (const ComplexDiagMatrix& a) const; + + ComplexDiagMatrix& fill (double val); + ComplexDiagMatrix& fill (const Complex& val); + ComplexDiagMatrix& fill (double val, int beg, int end); + ComplexDiagMatrix& fill (const Complex& val, int beg, int end); + ComplexDiagMatrix& fill (const ColumnVector& a); + ComplexDiagMatrix& fill (const ComplexColumnVector& a); + ComplexDiagMatrix& fill (const RowVector& a); + ComplexDiagMatrix& fill (const ComplexRowVector& a); + ComplexDiagMatrix& fill (const ColumnVector& a, int beg); + ComplexDiagMatrix& fill (const ComplexColumnVector& a, int beg); + ComplexDiagMatrix& fill (const RowVector& a, int beg); + ComplexDiagMatrix& fill (const ComplexRowVector& a, int beg); + + ComplexDiagMatrix hermitian (void) const; // complex conjugate transpose + ComplexDiagMatrix transpose (void) const; + + friend DiagMatrix real (const ComplexDiagMatrix& a); + friend DiagMatrix imag (const ComplexDiagMatrix& a); + friend ComplexDiagMatrix conj (const ComplexDiagMatrix& a); + +// resize is the destructive analog for this one + + ComplexMatrix extract (int r1, int c1, int r2, int c2) const; + +// extract row or column i. + + ComplexRowVector row (int i) const; + ComplexRowVector row (char *s) const; + + ComplexColumnVector column (int i) const; + ComplexColumnVector column (char *s) const; + + ComplexDiagMatrix inverse (int& info) const; + ComplexDiagMatrix inverse (void) const; + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + + ComplexDiagMatrix& operator += (const DiagMatrix& a); + ComplexDiagMatrix& operator -= (const DiagMatrix& a); + + ComplexDiagMatrix& operator += (const ComplexDiagMatrix& a); + ComplexDiagMatrix& operator -= (const ComplexDiagMatrix& a); + +// diagonal matrix by scalar -> matrix operations + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, double s); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, double s); + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, + const Complex& s); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, + const Complex& s); + +// diagonal matrix by scalar -> diagonal matrix operations + + friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, double s); + friend ComplexDiagMatrix operator / (const ComplexDiagMatrix& a, double s); + +// scalar by diagonal matrix -> matrix operations + + friend ComplexMatrix operator + (double s, const ComplexDiagMatrix& a); + friend ComplexMatrix operator - (double s, const ComplexDiagMatrix& a); + + friend ComplexMatrix operator + (const Complex& s, + const ComplexDiagMatrix& a); + friend ComplexMatrix operator - (const Complex& s, + const ComplexDiagMatrix& a); + +// scalar by diagonal matrix -> diagonal matrix operations + + friend ComplexDiagMatrix operator * (double s, const ComplexDiagMatrix& a); + +// diagonal matrix by column vector -> column vector operations + + friend ComplexColumnVector operator * (const ComplexDiagMatrix& a, + const ColumnVector& b); + + friend ComplexColumnVector operator * (const ComplexDiagMatrix& a, + const ComplexColumnVector& b); + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + + friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, + const ComplexDiagMatrix& b); + + friend ComplexDiagMatrix operator + (const ComplexDiagMatrix& a, + const DiagMatrix& b); + friend ComplexDiagMatrix operator - (const ComplexDiagMatrix& a, + const DiagMatrix& b); + friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, + const DiagMatrix& b); + + friend ComplexDiagMatrix product (const ComplexDiagMatrix& a, + const DiagMatrix& b); + +// diagonal matrix by matrix -> matrix operations + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, + const Matrix& b); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, + const Matrix& b); + friend ComplexMatrix operator * (const ComplexDiagMatrix& a, + const Matrix& b); + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator * (const ComplexDiagMatrix& a, + const ComplexMatrix& b); + +// other operations + + ComplexColumnVector diag (void) const; + ComplexColumnVector diag (int k) const; + +// i/o + + friend ostream& operator << (ostream& os, const ComplexDiagMatrix& a); + +#define KLUDGE_DIAG_MATRICES +#define TYPE Complex +#define KL_DMAT_TYPE ComplexDiagMatrix +#include "mx-kludge.h" +#undef KLUDGE_DIAG_MATRICES +#undef TYPE +#undef KL_DMAT_TYPE + +private: + + ComplexDiagMatrix (Complex *d, int nr, int nc) + : DiagArray (d, nr, nc) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc new file mode 100644 --- /dev/null +++ b/liboctave/CMatrix.cc @@ -0,0 +1,2613 @@ +// Matrix manipulations. -*- 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 +#include + +#include + +#include "mx-base.h" +#include "CmplxDET.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (zgemm) (const char*, const char*, const int*, + const int*, const int*, const Complex*, + const Complex*, const int*, const Complex*, + const int*, const Complex*, Complex*, const int*, + long, long); + + int F77_FCN (zgemv) (const char*, const int*, const int*, + const Complex*, const Complex*, const int*, + const Complex*, const int*, const Complex*, + Complex*, const int*, long); + + int F77_FCN (zgeco) (Complex*, const int*, const int*, int*, + double*, Complex*); + + int F77_FCN (zgedi) (Complex*, const int*, const int*, int*, + Complex*, Complex*, const int*); + + int F77_FCN (zgesl) (Complex*, const int*, const int*, int*, + Complex*, const int*); + + int F77_FCN (zgelss) (const int*, const int*, const int*, Complex*, + const int*, Complex*, const int*, double*, + const double*, int*, Complex*, const int*, + double*, int*); + +// Note that the original complex fft routines were not written for +// double complex arguments. They have been modified by adding an +// implicit double precision (a-h,o-z) statement at the beginning of +// each subroutine. + + int F77_FCN (cffti) (const int*, Complex*); + + int F77_FCN (cfftf) (const int*, Complex*, Complex*); + + int F77_FCN (cfftb) (const int*, Complex*, Complex*); +} + +#define KLUDGE_MATRICES +#define TYPE Complex +#define KL_MAT_TYPE ComplexMatrix +#include "mx-kludge.cc" +#undef KLUDGE_MATRICES +#undef TYPE +#undef KL_MAT_TYPE + +/* + * Complex Matrix class + */ + +ComplexMatrix::ComplexMatrix (const Matrix& a) + : Array2 (a.rows (), a.cols ()) +{ + for (int j = 0; j < cols (); j++) + for (int i = 0; i < rows (); i++) + elem (i, j) = a.elem (i, j); +} + +ComplexMatrix::ComplexMatrix (const DiagMatrix& a) + : Array2 (a.rows (), a.cols (), 0.0) +{ + for (int i = 0; i < a.length (); i++) + elem (i, i) = a.elem (i, i); +} + +ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a) + : Array2 (a.rows (), a.cols (), 0.0) +{ + for (int i = 0; i < a.length (); i++) + elem (i, i) = a.elem (i, i); +} + +#if 0 +ComplexMatrix& +ComplexMatrix::resize (int r, int c) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r * c; + Complex* new_data = (Complex *) NULL; + if (new_len > 0) + { + new_data = new Complex [new_len]; + + int min_r = nr < r ? nr : r; + int min_c = nc < c ? nc : c; + + for (int j = 0; j < min_c; j++) + for (int i = 0; i < min_r; i++) + new_data[r*j+i] = elem (i, j); + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} + +ComplexMatrix& +ComplexMatrix::resize (int r, int c, double val) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r * c; + Complex *new_data = (Complex *) NULL; + if (new_len > 0) + { + new_data = new Complex [new_len]; + +// There may be faster or cleaner ways to do this. + + if (r > nr || c > nc) + copy (new_data, new_len, val); + + int min_r = nr < r ? nr : r; + int min_c = nc < c ? nc : c; + + for (int j = 0; j < min_c; j++) + for (int i = 0; i < min_r; i++) + new_data[r*j+i] = elem (i, j); + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} + +ComplexMatrix& +ComplexMatrix::resize (int r, int c, const Complex& val) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r * c; + Complex *new_data = (Complex *) NULL; + if (new_len > 0) + { + new_data = new Complex [new_len]; + +// There may be faster or cleaner ways to do this. + + if (r > nr || c > nc) + copy (new_data, new_len, val); + + int min_r = nr < r ? nr : r; + int min_c = nc < c ? nc : c; + + for (int j = 0; j < min_c; j++) + for (int i = 0; i < min_r; i++) + new_data[r*j+i] = elem (i, j); + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} +#endif + +int +ComplexMatrix::operator == (const ComplexMatrix& a) const +{ + if (rows () != a.rows () || cols () != a.cols ()) + return 0; + + return equal (data (), a.data (), length ()); +} + +int +ComplexMatrix::operator != (const ComplexMatrix& a) const +{ + return !(*this == a); +} + +// destructive insert/delete/reorder operations + +ComplexMatrix& +ComplexMatrix::insert (const Matrix& a, int r, int c) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (r < 0 || r + a_nr - 1 > rows () || c < 0 || c + a_nc - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + elem (r+i, c+j) = a.elem (i, j); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const RowVector& a, int r, int c) +{ + int a_len = a.length (); + if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r, c+i) = a.elem (i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const ColumnVector& a, int r, int c) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i, c) = a.elem (i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const DiagMatrix& a, int r, int c) +{ + if (r < 0 || r + a.rows () - 1 > rows () + || c < 0 || c + a.cols () - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a.length (); i++) + elem (r+i, c+i) = a.elem (i, i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const ComplexMatrix& a, int r, int c) +{ + int a_nr = a.rows (); + int a_nc = a.cols (); + if (r < 0 || r + a_nr - 1 > rows () || c < 0 || c + a_nc - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + elem (r+i, c+j) = a.elem (i, j); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const ComplexRowVector& a, int r, int c) +{ + int a_len = a.length (); + if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r, c+i) = a.elem (i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const ComplexColumnVector& a, int r, int c) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i, c) = a.elem (i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::insert (const ComplexDiagMatrix& a, int r, int c) +{ + if (r < 0 || r + a.rows () - 1 > rows () + || c < 0 || c + a.cols () - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a.length (); i++) + elem (r+i, c+i) = a.elem (i, i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::fill (double val) +{ + int nr = rows (); + int nc = cols (); + if (nr > 0 && nc > 0) + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + elem (i, j) = val; + + return *this; +} + +ComplexMatrix& +ComplexMatrix::fill (const Complex& val) +{ + int nr = rows (); + int nc = cols (); + if (nr > 0 && nc > 0) + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + elem (i, j) = val; + + return *this; +} + +ComplexMatrix& +ComplexMatrix::fill (double val, int r1, int c1, int r2, int c2) +{ + int nr = rows (); + int nc = cols (); + if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 + || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + for (int j = c1; j <= c2; j++) + for (int i = r1; i <= r2; i++) + elem (i, j) = val; + + return *this; +} + +ComplexMatrix& +ComplexMatrix::fill (const Complex& val, int r1, int c1, int r2, int c2) +{ + int nr = rows (); + int nc = cols (); + if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 + || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + for (int j = c1; j <= c2; j++) + for (int i = r1; i <= r2; i++) + elem (i, j) = val; + + return *this; +} + +ComplexMatrix +ComplexMatrix::append (const Matrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + a.cols ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const RowVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != 1) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + a.length ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const ColumnVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.length ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + 1); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const DiagMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + a.cols ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const ComplexMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + a.cols ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const ComplexRowVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != 1) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + a.length ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const ComplexColumnVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.length ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + 1); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::append (const ComplexDiagMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + ComplexMatrix retval (nr, nc + a.cols ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const Matrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + a.rows (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const RowVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.length ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + 1, nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const ColumnVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != 1) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + a.length (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const DiagMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + a.rows (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const ComplexMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + a.rows (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const ComplexRowVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.length ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + 1, nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const ComplexColumnVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != 1) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + a.length (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::stack (const ComplexDiagMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return *this; + } + + int nr_insert = nr; + ComplexMatrix retval (nr + a.rows (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +ComplexMatrix +ComplexMatrix::hermitian (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix result; + if (length () > 0) + { + result.resize (nc, nr); + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + result.elem (j, i) = conj (elem (i, j)); + } + return result; +} + +ComplexMatrix +ComplexMatrix::transpose (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix result (nc, nr); + if (length () > 0) + { + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + result.elem (j, i) = elem (i, j); + } + return result; +} + +Matrix +real (const ComplexMatrix& a) +{ + int a_len = a.length (); + Matrix retval; + if (a_len > 0) + retval = Matrix (real_dup (a.data (), a_len), a.rows (), a.cols ()); + return retval; +} + +Matrix +imag (const ComplexMatrix& a) +{ + int a_len = a.length (); + Matrix retval; + if (a_len > 0) + retval = Matrix (imag_dup (a.data (), a_len), a.rows (), a.cols ()); + return retval; +} + +ComplexMatrix +conj (const ComplexMatrix& a) +{ + int a_len = a.length (); + ComplexMatrix retval; + if (a_len > 0) + retval = ComplexMatrix (conj_dup (a.data (), a_len), a.rows (), + a.cols ()); + return retval; +} + +// resize is the destructive equivalent for this one + +ComplexMatrix +ComplexMatrix::extract (int r1, int c1, int r2, int c2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + int new_r = r2 - r1 + 1; + int new_c = c2 - c1 + 1; + + ComplexMatrix result (new_r, new_c); + + for (int j = 0; j < new_c; j++) + for (int i = 0; i < new_r; i++) + result.elem (i, j) = elem (r1+i, c1+j); + + return result; +} + +// extract row or column i. + +ComplexRowVector +ComplexMatrix::row (int i) const +{ + int nc = cols (); + if (i < 0 || i >= rows ()) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return ComplexRowVector (); + } + + ComplexRowVector retval (nc); + for (int j = 0; j < cols (); j++) + retval.elem (j) = elem (i, j); + + return retval; +} + +ComplexRowVector +ComplexMatrix::row (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return ComplexRowVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return row (0); + else if (c == 'l' || c == 'L') + return row (rows () - 1); + else + { + (*current_liboctave_error_handler) ("invalid row selection"); + return ComplexRowVector (); + } +} + +ComplexColumnVector +ComplexMatrix::column (int i) const +{ + int nr = rows (); + if (i < 0 || i >= cols ()) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ComplexColumnVector (); + } + + ComplexColumnVector retval (nr); + for (int j = 0; j < nr; j++) + retval.elem (j) = elem (j, i); + + return retval; +} + +ComplexColumnVector +ComplexMatrix::column (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ComplexColumnVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return column (0); + else if (c == 'l' || c == 'L') + return column (cols () - 1); + else + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ComplexColumnVector (); + } +} + +ComplexMatrix +ComplexMatrix::inverse (void) const +{ + int info; + double rcond; return inverse (info, rcond); +} + +ComplexMatrix +ComplexMatrix::inverse (int& info) const +{ + double rcond; + return inverse (info, rcond); +} + +ComplexMatrix +ComplexMatrix::inverse (int& info, double& rcond) const +{ + int nr = rows (); + int nc = cols (); + int len = length (); + if (nr != nc) + { + (*current_liboctave_error_handler) ("inverse requires square matrix"); + return ComplexMatrix (); + } + + info = 0; + + int *ipvt = new int [nr]; + Complex *z = new Complex [nr]; + Complex *tmp_data = dup (data (), len); + + F77_FCN (zgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -1; + copy (tmp_data, data (), len); // Restore contents. + } + else + { + int job = 1; + Complex dummy; + + F77_FCN (zgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job); + } + + delete [] ipvt; + delete [] z; + + return ComplexMatrix (tmp_data, nr, nc); +} + +ComplexMatrix +ComplexMatrix::fourier (void) const +{ + int nr = rows (); + int nc = cols (); + int npts, nsamples; + if (nr == 1 || nc == 1) + { + npts = nr > nc ? nr : nc; + nsamples = 1; + } + else + { + npts = nr; + nsamples = nc; + } + + int nn = 4*npts+15; + Complex *wsave = new Complex [nn]; + Complex *tmp_data = dup (data (), length ()); + + F77_FCN (cffti) (&npts, wsave); + + for (int j = 0; j < nsamples; j++) + F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave); + + delete [] wsave; + + return ComplexMatrix (tmp_data, nr, nc); +} + +ComplexMatrix +ComplexMatrix::ifourier (void) const +{ + int nr = rows (); + int nc = cols (); + int npts, nsamples; + if (nr == 1 || nc == 1) + { + npts = nr > nc ? nr : nc; + nsamples = 1; + } + else + { + npts = nr; + nsamples = nc; + } + + int nn = 4*npts+15; + Complex *wsave = new Complex [nn]; + Complex *tmp_data = dup (data (), length ()); + + F77_FCN (cffti) (&npts, wsave); + + for (int j = 0; j < nsamples; j++) + F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave); + + for (j = 0; j < npts*nsamples; j++) + tmp_data[j] = tmp_data[j] / (double) npts; + + delete [] wsave; + + return ComplexMatrix (tmp_data, nr, nc); +} + +ComplexDET +ComplexMatrix::determinant (void) const +{ + int info; + double rcond; + return determinant (info, rcond); +} + +ComplexDET +ComplexMatrix::determinant (int& info) const +{ + double rcond; + return determinant (info, rcond); +} + +ComplexDET +ComplexMatrix::determinant (int& info, double& rcond) const +{ + ComplexDET retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 0 || nc == 0) + { + Complex d[2]; + d[0] = 1.0; + d[1] = 0.0; + retval = ComplexDET (d); + } + else + { + info = 0; + int *ipvt = new int [nr]; + + Complex *z = new Complex [nr]; + Complex *tmp_data = dup (data (), length ()); + + F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -1; + retval = ComplexDET (); + } + else + { + int job = 10; + Complex d[2]; + F77_FCN (zgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job); + retval = ComplexDET (d); + } + + delete [] tmp_data; + delete [] ipvt; + delete [] z; + } + + return retval; +} + +ComplexMatrix +ComplexMatrix::solve (const Matrix& b) const +{ + int info; + double rcond; + return solve (b, info, rcond); +} + +ComplexMatrix +ComplexMatrix::solve (const Matrix& b, int& info) const +{ + double rcond; + return solve (b, info, rcond); +} + +ComplexMatrix +ComplexMatrix::solve (const Matrix& b, int& info, double& rcond) const +{ + ComplexMatrix tmp (b); + return solve (tmp, info, rcond); +} + +ComplexMatrix +ComplexMatrix::solve (const ComplexMatrix& b) const +{ + int info; + double rcond; + return solve (b, info, rcond); +} + +ComplexMatrix +ComplexMatrix::solve (const ComplexMatrix& b, int& info) const +{ + double rcond; + return solve (b, info, rcond); +} +ComplexMatrix +ComplexMatrix::solve (const ComplexMatrix& b, int& info, double& rcond) const +{ + ComplexMatrix retval; + + int nr = rows (); + int nc = cols (); + int b_nr = b.rows (); + int b_nc = b.cols (); + if (nr == 0 || nc == 0 || nr != nc || nr != b_nr) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of linear equations"); + return ComplexMatrix (); + } + + info = 0; + int *ipvt = new int [nr]; + + Complex *z = new Complex [nr]; + Complex *tmp_data = dup (data (), length ()); + + F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -2; + } + else + { + int job = 0; + + Complex *result = dup (b.data (), b.length ()); + + for (int j = 0; j < b_nc; j++) + F77_FCN (zgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job); + + retval = ComplexMatrix (result, b_nr, b_nc); + } + + delete [] tmp_data; + delete [] ipvt; + delete [] z; + + return retval; +} + +ComplexColumnVector +ComplexMatrix::solve (const ComplexColumnVector& b) const +{ + int info; + double rcond; + return solve (b, info, rcond); +} + +ComplexColumnVector +ComplexMatrix::solve (const ComplexColumnVector& b, int& info) const +{ + double rcond; + return solve (b, info, rcond); +} + +ComplexColumnVector +ComplexMatrix::solve (const ComplexColumnVector& b, int& info, + double& rcond) const +{ + ComplexColumnVector retval; + + int nr = rows (); + int nc = cols (); + int b_len = b.length (); + if (nr == 0 || nc == 0 || nr != nc || nr != b_len) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of linear equations"); + return ComplexColumnVector (); + } + + info = 0; + int *ipvt = new int [nr]; + + Complex *z = new Complex [nr]; + Complex *tmp_data = dup (data (), length ()); + + F77_FCN (zgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -2; + } + else + { + int job = 0; + + Complex *result = dup (b.data (), b_len); + + F77_FCN (zgesl) (tmp_data, &nr, &nr, ipvt, result, &job); + + retval = ComplexColumnVector (result, b_len); + } + + delete [] tmp_data; + delete [] ipvt; + delete [] z; + + return retval; +} + +ComplexMatrix +ComplexMatrix::lssolve (const ComplexMatrix& b) const +{ + int info; + int rank; + return lssolve (b, info, rank); +} + +ComplexMatrix +ComplexMatrix::lssolve (const ComplexMatrix& b, int& info) const +{ + int rank; + return lssolve (b, info, rank); +} + +ComplexMatrix +ComplexMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const +{ + int nrhs = b.cols (); + + int m = rows (); + int n = cols (); + + if (m == 0 || n == 0 || m != b.rows ()) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch solution of linear equations"); + return Matrix (); + } + + Complex *tmp_data = dup (data (), length ()); + + int nrr = m > n ? m : n; + ComplexMatrix result (nrr, nrhs); + + int i, j; + for (j = 0; j < nrhs; j++) + for (i = 0; i < m; i++) + result.elem (i, j) = b.elem (i, j); + + Complex *presult = result.fortran_vec (); + + int len_s = m < n ? m : n; + double *s = new double [len_s]; + double rcond = -1.0; + int lwork; + if (m < n) + lwork = 2*m + (nrhs > n ? nrhs : n); + else + lwork = 2*n + (nrhs > m ? nrhs : m); + + Complex *work = new Complex [lwork]; + + int lrwork = (5 * (m < n ? m : n)) - 4; + lrwork = lrwork > 1 ? lrwork : 1; + double *rwork = new double [lrwork]; + + F77_FCN (zgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, + &rcond, &rank, work, &lwork, rwork, &info); + + ComplexMatrix retval (n, nrhs); + for (j = 0; j < nrhs; j++) + for (i = 0; i < n; i++) + retval.elem (i, j) = result.elem (i, j); + + delete [] tmp_data; + delete [] s; + delete [] work; + delete [] rwork; + + return retval; +} + +ComplexColumnVector +ComplexMatrix::lssolve (const ComplexColumnVector& b) const +{ + int info; + int rank; + return lssolve (b, info, rank); +} + +ComplexColumnVector +ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info) const +{ + int rank; + return lssolve (b, info, rank); +} + +ComplexColumnVector +ComplexMatrix::lssolve (const ComplexColumnVector& b, int& info, + int& rank) const +{ + int nrhs = 1; + + int m = rows (); + int n = cols (); + + if (m == 0 || n == 0 || m != b.length ()) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch solution of least squares problem"); + return ComplexColumnVector (); + } + + Complex *tmp_data = dup (data (), length ()); + + int nrr = m > n ? m : n; + ComplexColumnVector result (nrr); + + int i; + for (i = 0; i < m; i++) + result.elem (i) = b.elem (i); + + Complex *presult = result.fortran_vec (); + + int len_s = m < n ? m : n; + double *s = new double [len_s]; + double rcond = -1.0; + int lwork; + if (m < n) + lwork = 2*m + (nrhs > n ? nrhs : n); + else + lwork = 2*n + (nrhs > m ? nrhs : m); + + Complex *work = new Complex [lwork]; + + int lrwork = (5 * (m < n ? m : n)) - 4; + lrwork = lrwork > 1 ? lrwork : 1; + double *rwork = new double [lrwork]; + + F77_FCN (zgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, + &rcond, &rank, work, &lwork, rwork, &info); + + ComplexColumnVector retval (n); + for (i = 0; i < n; i++) + retval.elem (i) = result.elem (i); + + delete [] tmp_data; + delete [] s; + delete [] work; + delete [] rwork; + + return retval; +} + +// matrix by diagonal matrix -> matrix operations + +ComplexMatrix& +ComplexMatrix::operator += (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return ComplexMatrix (); + } + + for (int i = 0; i < a.length (); i++) + elem (i, i) += a.elem (i, i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::operator -= (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return ComplexMatrix (); + } + + for (int i = 0; i < a.length (); i++) + elem (i, i) -= a.elem (i, i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::operator += (const ComplexDiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return ComplexMatrix (); + } + + for (int i = 0; i < a.length (); i++) + elem (i, i) += a.elem (i, i); + + return *this; +} + +ComplexMatrix& +ComplexMatrix::operator -= (const ComplexDiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return ComplexMatrix (); + } + + for (int i = 0; i < a.length (); i++) + elem (i, i) -= a.elem (i, i); + + return *this; +} + +// matrix by matrix -> matrix operations + +ComplexMatrix& +ComplexMatrix::operator += (const Matrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + return *this; +} + +ComplexMatrix& +ComplexMatrix::operator -= (const Matrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + return *this; +} + +ComplexMatrix& +ComplexMatrix::operator += (const ComplexMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + return *this; +} + +ComplexMatrix& +ComplexMatrix::operator -= (const ComplexMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + return *this; +} + +// unary operations + +Matrix +ComplexMatrix::operator ! (void) const +{ + return Matrix (not (data (), length ()), rows (), cols ()); +} + +// matrix by scalar -> matrix operations + +ComplexMatrix +operator + (const ComplexMatrix& a, double s) +{ + return ComplexMatrix (add (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator - (const ComplexMatrix& a, double s) +{ + return ComplexMatrix (subtract (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator * (const ComplexMatrix& a, double s) +{ + return ComplexMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator / (const ComplexMatrix& a, double s) +{ + return ComplexMatrix (divide (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +// scalar by matrix -> matrix operations + +ComplexMatrix +operator + (double s, const ComplexMatrix& a) +{ + return ComplexMatrix (add (a.data (), a.length (), s), a.rows (), + a.cols ()); +} + +ComplexMatrix +operator - (double s, const ComplexMatrix& a) +{ + return ComplexMatrix (subtract (s, a.data (), a.length ()), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator * (double s, const ComplexMatrix& a) +{ + return ComplexMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator / (double s, const ComplexMatrix& a) +{ + return ComplexMatrix (divide (s, a.data (), a.length ()), + a.rows (), a.cols ()); +} + +// matrix by column vector -> column vector operations + +ComplexColumnVector +operator * (const ComplexMatrix& m, const ColumnVector& a) +{ + ComplexColumnVector tmp (a); + return m * tmp; +} + +ComplexColumnVector +operator * (const ComplexMatrix& m, const ComplexColumnVector& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nc != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexColumnVector (); + } + + if (nc == 0 || nr == 0) + return ComplexColumnVector (0); + + char trans = 'N'; + int ld = nr; + Complex alpha (1.0); + Complex beta (0.0); + int i_one = 1; + + Complex *y = new Complex [nr]; + + F77_FCN (zgemv) (&trans, &nr, &nc, &alpha, m.data (), &ld, a.data (), + &i_one, &beta, y, &i_one, 1L); + + return ComplexColumnVector (y, nr); +} + +// matrix by diagonal matrix -> matrix operations + +ComplexMatrix +operator + (const ComplexMatrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (m); + for (int i = 0; i < a.length (); i++) + result.elem (i, i) += a.elem (i, i); + + return result; +} + +ComplexMatrix +operator - (const ComplexMatrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (m); + for (int i = 0; i < a.length (); i++) + result.elem (i, i) -= a.elem (i, i); + + return result; +} + +ComplexMatrix +operator * (const ComplexMatrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nc = a.cols (); + if (nc != a.rows ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, nc, 0.0); + + Complex *c = new Complex [nr*a_nc]; + Complex *ctmp = (Complex *) NULL; + + for (int j = 0; j < a.length (); j++) + { + int idx = j * nr; + ctmp = c + idx; + if (a.elem (j, j) == 1.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = m.elem (i, j); + } + else if (a.elem (j, j) == 0.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = 0.0; + } + else + { + for (int i = 0; i < nr; i++) + ctmp[i] = a.elem (j, j) * m.elem (i, j); + } + } + + if (a.rows () < a_nc) + { + for (int i = nr * nc; i < nr * a_nc; i++) + ctmp[i] = 0.0; + } + + return ComplexMatrix (c, nr, a_nc); +} + +ComplexMatrix +operator + (const ComplexMatrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (m); + for (int i = 0; i < a.length (); i++) + result.elem (i, i) += a.elem (i, i); + + return result; +} + +ComplexMatrix +operator - (const ComplexMatrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (m); + for (int i = 0; i < a.length (); i++) + result.elem (i, i) -= a.elem (i, i); + + return result; +} + +ComplexMatrix +operator * (const ComplexMatrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nc = a.cols (); + if (nc != a.rows ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, nc, 0.0); + + Complex *c = new Complex [nr*a_nc]; + Complex *ctmp = (Complex *) NULL; + + for (int j = 0; j < a.length (); j++) + { + int idx = j * nr; + ctmp = c + idx; + if (a.elem (j, j) == 1.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = m.elem (i, j); + } + else if (a.elem (j, j) == 0.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = 0.0; + } + else + { + for (int i = 0; i < nr; i++) + ctmp[i] = a.elem (j, j) * m.elem (i, j); + } + } + + if (a.rows () < a_nc) + { + for (int i = nr * nc; i < nr * a_nc; i++) + ctmp[i] = 0.0; + } + + return ComplexMatrix (c, nr, a_nc); +} + +// matrix by matrix -> matrix operations + +ComplexMatrix +operator + (const ComplexMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexMatrix +operator - (const ComplexMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexMatrix +operator * (const ComplexMatrix& m, const Matrix& a) +{ + ComplexMatrix tmp (a); + return m * tmp; +} + +ComplexMatrix +operator * (const ComplexMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nc = a.cols (); + if (nc != a.rows ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, nc, 0.0); + + char trans = 'N'; + char transa = 'N'; + + int ld = nr; + int lda = a.rows (); + + Complex alpha (1.0); + Complex beta (0.0); + + Complex *c = new Complex [nr*a_nc]; + + F77_FCN (zgemm) (&trans, &transa, &nr, &a_nc, &nc, &alpha, m.data (), + &ld, a.data (), &lda, &beta, c, &nr, 1L, 1L); + + return ComplexMatrix (c, nr, a_nc); +} + +ComplexMatrix +product (const ComplexMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix product attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexMatrix +quotient (const ComplexMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix quotient attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc); +} + +// other operations + +ComplexMatrix +map (c_c_Mapper f, const ComplexMatrix& a) +{ + ComplexMatrix b (a); + b.map (f); + return b; +} + +Matrix +map (d_c_Mapper f, const ComplexMatrix& a) +{ + int a_nc = a.cols (); + int a_nr = a.rows (); + Matrix b (a_nr, a_nc); + for (int j = 0; j < a_nc; j++) + for (int i = 0; i < a_nr; i++) + b.elem (i, j) = f (a.elem (i, j)); + return b; +} + +void +ComplexMatrix::map (c_c_Mapper f) +{ + for (int j = 0; j < cols (); j++) + for (int i = 0; i < rows (); i++) + elem (i, j) = f (elem (i, j)); +} + +Matrix +ComplexMatrix::all (void) const +{ + int nr = rows (); + int nc = cols (); + Matrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int j = 0; j < nc; j++) + { + if (elem (0, j) == 0.0) + { + retval.elem (0, 0) = 0.0; + break; + } + } + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, 0) == 0.0) + { + retval.elem (0, 0) = 0.0; + break; + } + } + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 1.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, j) == 0.0) + { + retval.elem (0, j) = 0.0; + break; + } + } + } + } + } + return retval; +} + +Matrix +ComplexMatrix::any (void) const +{ + int nr = rows (); + int nc = cols (); + Matrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int j = 0; j < nc; j++) + { + if (elem (0, j) != 0.0) + { + retval.elem (0, 0) = 1.0; + break; + } + } + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, 0) != 0.0) + { + retval.elem (0, 0) = 1.0; + break; + } + } + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 0.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, j) != 0.0) + { + retval.elem (0, j) = 1.0; + break; + } + } + } + } + } + return retval; +} + +ComplexMatrix +ComplexMatrix::cumprod (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, nc); + Complex prod = elem (0, 0); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = prod; + if (j < nc - 1) + prod *= elem (0, j+1); + } + } + else if (nc == 1) + { + retval.resize (nr, 1); + Complex prod = elem (0, 0); + for (int i = 0; i < nr; i++) + { + retval.elem (i, 0) = prod; + if (i < nr - 1) + prod *= elem (i+1, 0); + } + } + else + { + retval.resize (nr, nc); + for (int j = 0; j < nc; j++) + { + Complex prod = elem (0, j); + for (int i = 0; i < nr; i++) + { + retval.elem (i, j) = prod; + if (i < nr - 1) + prod *= elem (i+1, j); + } + } + } + } + return retval; +} + +ComplexMatrix +ComplexMatrix::cumsum (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, nc); + Complex sum = elem (0, 0); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = sum; + if (j < nc - 1) + sum += elem (0, j+1); + } + } + else if (nc == 1) + { + retval.resize (nr, 1); + Complex sum = elem (0, 0); + for (int i = 0; i < nr; i++) + { + retval.elem (i, 0) = sum; + if (i < nr - 1) + sum += elem (i+1, 0); + } + } + else + { + retval.resize (nr, nc); + for (int j = 0; j < nc; j++) + { + Complex sum = elem (0, j); + for (int i = 0; i < nr; i++) + { + retval.elem (i, j) = sum; + if (i < nr - 1) + sum += elem (i+1, j); + } + } + } + } + return retval; +} + +ComplexMatrix +ComplexMatrix::prod (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int j = 0; j < nc; j++) + retval.elem (0, 0) *= elem (0, j); + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int i = 0; i < nr; i++) + retval.elem (0, 0) *= elem (i, 0); + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 1.0; + for (int i = 0; i < nr; i++) + retval.elem (0, j) *= elem (i, j); + } + } + } + return retval; +} + +ComplexMatrix +ComplexMatrix::sum (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int j = 0; j < nc; j++) + retval.elem (0, 0) += elem (0, j); + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int i = 0; i < nr; i++) + retval.elem (0, 0) += elem (i, 0); + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 0.0; + for (int i = 0; i < nr; i++) + retval.elem (0, j) += elem (i, j); + } + } + } + return retval; +} + +ComplexMatrix +ComplexMatrix::sumsq (void) const +{ + int nr = rows (); + int nc = cols (); + ComplexMatrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int j = 0; j < nc; j++) + { + Complex d = elem (0, j); + retval.elem (0, 0) += d * d; + } + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int i = 0; i < nr; i++) + { + Complex d = elem (i, 0); + retval.elem (0, 0) += d * d; + } + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 0.0; + for (int i = 0; i < nr; i++) + { + Complex d = elem (i, j); + retval.elem (0, j) += d * d; + } + } + } + } + return retval; +} + +ComplexColumnVector +ComplexMatrix::diag (void) const +{ + return diag (0); +} + +ComplexColumnVector +ComplexMatrix::diag (int k) const +{ + int nnr = rows (); + int nnc = cols (); + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + ComplexColumnVector d; + + if (nnr > 0 && nnc > 0) + { + int ndiag = (nnr < nnc) ? nnr : nnc; + + d.resize (ndiag); + + if (k > 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i+k); + } + else if ( k < 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i-k, i); + } + else + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i); + } + } + else + cerr << "diag: requested diagonal out of range\n"; + + return d; +} + +ComplexColumnVector +ComplexMatrix::row_min (void) const +{ + ComplexColumnVector result; + + int nr = rows (); + int nc = cols (); + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + Complex res = elem (i, 0); + double absres = abs (res); + for (int j = 1; j < nc; j++) + if (abs (elem (i, j)) < absres) + { + res = elem (i, j); + absres = abs (res); + } + result.elem (i) = res; + } + } + + return result; +} + +ComplexColumnVector +ComplexMatrix::row_min_loc (void) const +{ + ComplexColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + Complex res = 0; + double absres = abs (elem (i, 0)); + for (int j = 0; j < nc; j++) + if (abs (elem (i, j)) < absres) + { + res = j; + absres = abs (elem (i, j)); + } + result.elem (i) = res + 1; + } + } + + return result; +} + +ComplexColumnVector +ComplexMatrix::row_max (void) const +{ + ComplexColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + Complex res = elem (i, 0); + double absres = abs (res); + for (int j = 1; j < nc; j++) + if (abs (elem (i, j)) > absres) + { + res = elem (i, j); + absres = abs (res); + } + result.elem (i) = res; + } + } + + return result; +} + +ComplexColumnVector +ComplexMatrix::row_max_loc (void) const +{ + ComplexColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + Complex res = 0; + double absres = abs (elem (i, 0)); + for (int j = 0; j < nc; j++) + if (abs (elem (i, j)) > absres) + { + res = j; + absres = abs (elem (i, j)); + } + result.elem (i) = res + 1; + } + } + + return result; +} + +ComplexRowVector +ComplexMatrix::column_min (void) const +{ + ComplexRowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + Complex res = elem (0, j); + double absres = abs (res); + for (int i = 1; i < nr; i++) + if (abs (elem (i, j)) < absres) + { + res = elem (i, j); + absres = abs (res); + } + result.elem (j) = res; + } + } + + return result; +} + +ComplexRowVector +ComplexMatrix::column_min_loc (void) const +{ + ComplexRowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + Complex res = 0; + double absres = abs (elem (0, j)); + for (int i = 0; i < nr; i++) + if (abs (elem (i, j)) < absres) + { + res = i; + absres = abs (elem (i, j)); + } + result.elem (j) = res + 1; + } + } + + return result; +} + +ComplexRowVector +ComplexMatrix::column_max (void) const +{ + ComplexRowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + Complex res = elem (0, j); + double absres = abs (res); + for (int i = 1; i < nr; i++) + if (abs (elem (i, j)) > absres) + { + res = elem (i, j); + absres = abs (res); + } + result.elem (j) = res; + } + } + + return result; +} + +ComplexRowVector +ComplexMatrix::column_max_loc (void) const +{ + ComplexRowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + Complex res = 0; + double absres = abs (elem (0, j)); + for (int i = 0; i < nr; i++) + if (abs (elem (i, j)) > absres) + { + res = i; + absres = abs (elem (i, j)); + } + result.elem (j) = res + 1; + } + } + + return result; +} + +// i/o + +ostream& +operator << (ostream& os, const ComplexMatrix& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.rows (); i++) + { + for (int j = 0; j < a.cols (); j++) + os << " " /* setw (field_width) */ << a.elem (i, j); + os << "\n"; + } + return os; +} + +istream& +operator >> (istream& is, ComplexMatrix& a) +{ + int nr = a.rows (); + int nc = a.cols (); + + if (nr < 1 || nc < 1) + is.clear (ios::badbit); + else + { + Complex tmp; + for (int i = 0; i < nr; i++) + for (int j = 0; j < nc; j++) + { + is >> tmp; + if (is) + a.elem (i, j) = tmp; + else + break; + } + } + + return is; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h new file mode 100644 --- /dev/null +++ b/liboctave/CMatrix.h @@ -0,0 +1,287 @@ +// -*- 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. + +*/ + +#if !defined (octave_ComplexMatrix_h) +#define octave_ComplexMatrix_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include + +#include "Array.h" + +#include "mx-defs.h" + +extern "C++" { + +class ComplexMatrix : public Array2 +{ +friend class ComplexLU; +friend class ComplexSVD; + +public: + + ComplexMatrix (void) : Array2 () { } + ComplexMatrix (int r, int c) : Array2 (r, c) { } + ComplexMatrix (int r, int c, const Complex& val) + : Array2 (r, c, val) { } + ComplexMatrix (const Matrix& a); + ComplexMatrix (const Array2& a) : Array2 (a) { } + ComplexMatrix (const ComplexMatrix& a) : Array2 (a) { } + ComplexMatrix (const DiagMatrix& a); + ComplexMatrix (const DiagArray& a) : Array2 (a) { } + ComplexMatrix (const ComplexDiagMatrix& a); +// ComplexMatrix (double a) : Array2 (1, 1, a) { } +// ComplexMatrix (const Complex& a) : Array2 (1, 1, a) { } + + ComplexMatrix& operator = (const ComplexMatrix& a) + { + Array2::operator = (a); + return *this; + } + +// operator Array2& () const { return *this; } + + int operator == (const ComplexMatrix& a) const; + int operator != (const ComplexMatrix& a) const; + +// destructive insert/delete/reorder operations + + ComplexMatrix& insert (const Matrix& a, int r, int c); + ComplexMatrix& insert (const RowVector& a, int r, int c); + ComplexMatrix& insert (const ColumnVector& a, int r, int c); + ComplexMatrix& insert (const DiagMatrix& a, int r, int c); + + ComplexMatrix& insert (const ComplexMatrix& a, int r, int c); + ComplexMatrix& insert (const ComplexRowVector& a, int r, int c); + ComplexMatrix& insert (const ComplexColumnVector& a, int r, int c); + ComplexMatrix& insert (const ComplexDiagMatrix& a, int r, int c); + + ComplexMatrix& fill (double val); + ComplexMatrix& fill (const Complex& val); + ComplexMatrix& fill (double val, int r1, int c1, int r2, int c2); + ComplexMatrix& fill (const Complex& val, int r1, int c1, int r2, int c2); + + ComplexMatrix append (const Matrix& a) const; + ComplexMatrix append (const RowVector& a) const; + ComplexMatrix append (const ColumnVector& a) const; + ComplexMatrix append (const DiagMatrix& a) const; + + ComplexMatrix append (const ComplexMatrix& a) const; + ComplexMatrix append (const ComplexRowVector& a) const; + ComplexMatrix append (const ComplexColumnVector& a) const; + ComplexMatrix append (const ComplexDiagMatrix& a) const; + + ComplexMatrix stack (const Matrix& a) const; + ComplexMatrix stack (const RowVector& a) const; + ComplexMatrix stack (const ColumnVector& a) const; + ComplexMatrix stack (const DiagMatrix& a) const; + + ComplexMatrix stack (const ComplexMatrix& a) const; + ComplexMatrix stack (const ComplexRowVector& a) const; + ComplexMatrix stack (const ComplexColumnVector& a) const; + ComplexMatrix stack (const ComplexDiagMatrix& a) const; + + ComplexMatrix hermitian (void) const; // complex conjugate transpose + ComplexMatrix transpose (void) const; + + friend Matrix real (const ComplexMatrix& a); + friend Matrix imag (const ComplexMatrix& a); + friend ComplexMatrix conj (const ComplexMatrix& a); + +// resize is the destructive equivalent for this one + + ComplexMatrix extract (int r1, int c1, int r2, int c2) const; + +// extract row or column i. + + ComplexRowVector row (int i) const; + ComplexRowVector row (char *s) const; + + ComplexColumnVector column (int i) const; + ComplexColumnVector column (char *s) const; + + ComplexMatrix inverse (void) const; + ComplexMatrix inverse (int& info) const; + ComplexMatrix inverse (int& info, double& rcond) const; + + ComplexMatrix fourier (void) const; + ComplexMatrix ifourier (void) const; + + ComplexDET determinant (void) const; + ComplexDET determinant (int& info) const; + ComplexDET determinant (int& info, double& rcond) const; + + ComplexMatrix solve (const Matrix& b) const; + ComplexMatrix solve (const Matrix& b, int& info) const; + ComplexMatrix solve (const Matrix& b, int& info, double& rcond) const; + + ComplexMatrix solve (const ComplexMatrix& b) const; + ComplexMatrix solve (const ComplexMatrix& b, int& info) const; + ComplexMatrix solve (const ComplexMatrix& b, int& info, double& rcond) const; + + ComplexColumnVector solve (const ComplexColumnVector& b) const; + ComplexColumnVector solve (const ComplexColumnVector& b, int& info) const; + ComplexColumnVector solve (const ComplexColumnVector& b, int& info, + double& rcond) const; + + ComplexMatrix lssolve (const ComplexMatrix& b) const; + ComplexMatrix lssolve (const ComplexMatrix& b, int& info) const; + ComplexMatrix lssolve (const ComplexMatrix& b, int& info, + int& rank) const; + + ComplexColumnVector lssolve (const ComplexColumnVector& b) const; + ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info) const; + ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info, + int& rank) const; + +// matrix by diagonal matrix -> matrix operations + + ComplexMatrix& operator += (const DiagMatrix& a); + ComplexMatrix& operator -= (const DiagMatrix& a); + + ComplexMatrix& operator += (const ComplexDiagMatrix& a); + ComplexMatrix& operator -= (const ComplexDiagMatrix& a); + +// matrix by matrix -> matrix operations + + ComplexMatrix& operator += (const Matrix& a); + ComplexMatrix& operator -= (const Matrix& a); + + ComplexMatrix& operator += (const ComplexMatrix& a); + ComplexMatrix& operator -= (const ComplexMatrix& a); + +// unary operations + + Matrix operator ! (void) const; + +// matrix by scalar -> matrix operations + + friend ComplexMatrix operator + (const ComplexMatrix& a, double s); + friend ComplexMatrix operator - (const ComplexMatrix& a, double s); + friend ComplexMatrix operator * (const ComplexMatrix& a, double s); + friend ComplexMatrix operator / (const ComplexMatrix& a, double s); + +// scalar by matrix -> matrix operations + + friend ComplexMatrix operator + (double s, const ComplexMatrix& a); + friend ComplexMatrix operator - (double s, const ComplexMatrix& a); + friend ComplexMatrix operator * (double s, const ComplexMatrix& a); + friend ComplexMatrix operator / (double s, const ComplexMatrix& a); + +// matrix by column vector -> column vector operations + + friend ComplexColumnVector operator * (const ComplexMatrix& a, + const ColumnVector& b); + + friend ComplexColumnVector operator * (const ComplexMatrix& a, + const ComplexColumnVector& b); + +// matrix by diagonal matrix -> matrix operations + + friend ComplexMatrix operator + (const ComplexMatrix& a, + const DiagMatrix& b); + friend ComplexMatrix operator - (const ComplexMatrix& a, + const DiagMatrix& b); + friend ComplexMatrix operator * (const ComplexMatrix& a, + const DiagMatrix& b); + + friend ComplexMatrix operator + (const ComplexMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexMatrix operator - (const ComplexMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexMatrix operator * (const ComplexMatrix& a, + const ComplexDiagMatrix& b); + +// matrix by matrix -> matrix operations + + friend ComplexMatrix operator + (const ComplexMatrix& a, const Matrix& b); + friend ComplexMatrix operator - (const ComplexMatrix& a, const Matrix& b); + + friend ComplexMatrix operator * (const ComplexMatrix& a, const Matrix& b); + friend ComplexMatrix operator * (const ComplexMatrix& a, + const ComplexMatrix& b); + + friend ComplexMatrix product (const ComplexMatrix& a, const Matrix& b); + friend ComplexMatrix quotient (const ComplexMatrix& a, const Matrix& b); + +// other operations + + friend ComplexMatrix map (c_c_Mapper f, const ComplexMatrix& a); + friend Matrix map (d_c_Mapper f, const ComplexMatrix& a); + void map (c_c_Mapper f); + + Matrix all (void) const; + Matrix any (void) const; + + ComplexMatrix cumprod (void) const; + ComplexMatrix cumsum (void) const; + ComplexMatrix prod (void) const; + ComplexMatrix sum (void) const; + ComplexMatrix sumsq (void) const; + + ComplexColumnVector diag (void) const; + ComplexColumnVector diag (int k) const; + + ComplexColumnVector row_min (void) const; + ComplexColumnVector row_min_loc (void) const; + + ComplexColumnVector row_max (void) const; + ComplexColumnVector row_max_loc (void) const; + + ComplexRowVector column_min (void) const; + ComplexRowVector column_min_loc (void) const; + + ComplexRowVector column_max (void) const; + ComplexRowVector column_max_loc (void) const; + +// i/o + + friend ostream& operator << (ostream& os, const ComplexMatrix& a); + friend istream& operator >> (istream& is, ComplexMatrix& a); + +#define KLUDGE_MATRICES +#define TYPE Complex +#define KL_MAT_TYPE ComplexMatrix +#include "mx-kludge.h" +#undef KLUDGE_MATRICES +#undef TYPE +#undef KL_MAT_TYPE + +private: + + ComplexMatrix (Complex *d, int r, int c) : Array2 (d, r, c) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CRowVector.cc b/liboctave/CRowVector.cc new file mode 100644 --- /dev/null +++ b/liboctave/CRowVector.cc @@ -0,0 +1,688 @@ +// RowVector manipulations. -*- 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 + +#include + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (zgemv) (const char*, const int*, const int*, + const Complex*, const Complex*, const int*, + const Complex*, const int*, const Complex*, + Complex*, const int*, long); +} + +/* + * Complex Row Vector class + */ + +#define KLUDGE_VECTORS +#define TYPE Complex +#define KL_VEC_TYPE ComplexRowVector +#include "mx-kludge.cc" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +ComplexRowVector::ComplexRowVector (const RowVector& a) + : Array (a.length ()) +{ + for (int i = 0; i < length (); i++) + elem (i) = a.elem (i); +} + +#if 0 +ComplexRowVector& +ComplexRowVector::resize (int n) +{ + if (n < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimension"); + return *this; + } + + Complex *new_data = (Complex *) NULL; + if (n > 0) + { + new_data = new Complex [n]; + int min_len = len < n ? len : n; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + len = n; + data = new_data; + + return *this; +} + +ComplexRowVector& +ComplexRowVector::resize (int n, double val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} + +ComplexRowVector& +ComplexRowVector::resize (int n, const Complex& val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} +#endif + +int +ComplexRowVector::operator == (const ComplexRowVector& a) const +{ + int len = length (); + if (len != a.length ()) + return 0; + return equal (data (), a.data (), len); +} + +int +ComplexRowVector::operator != (const ComplexRowVector& a) const +{ + return !(*this == a); +} + +// destructive insert/delete/reorder operations + +ComplexRowVector& +ComplexRowVector::insert (const RowVector& a, int c) +{ + int a_len = a.length (); + if (c < 0 || c + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (c+i) = a.elem (i); + + return *this; +} + +ComplexRowVector& +ComplexRowVector::insert (const ComplexRowVector& a, int c) +{ + int a_len = a.length (); + if (c < 0 || c + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (c+i) = a.elem (i); + + return *this; +} + +ComplexRowVector& +ComplexRowVector::fill (double val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +ComplexRowVector& +ComplexRowVector::fill (const Complex& val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +ComplexRowVector& +ComplexRowVector::fill (double val, int c1, int c2) +{ + int len = length (); + if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + for (int i = c1; i <= c2; i++) + elem (i) = val; + + return *this; +} + +ComplexRowVector& +ComplexRowVector::fill (const Complex& val, int c1, int c2) +{ + int len = length (); + if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + for (int i = c1; i <= c2; i++) + elem (i) = val; + + return *this; +} + +ComplexRowVector +ComplexRowVector::append (const RowVector& a) const +{ + int len = length (); + int nc_insert = len; + ComplexRowVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nc_insert); + return retval; +} + +ComplexRowVector +ComplexRowVector::append (const ComplexRowVector& a) const +{ + int len = length (); + int nc_insert = len; + ComplexRowVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nc_insert); + return retval; +} + +ComplexColumnVector +ComplexRowVector::hermitian (void) const +{ + int len = length (); + return ComplexColumnVector (conj_dup (data (), len), len); +} + +ComplexColumnVector +ComplexRowVector::transpose (void) const +{ + int len = length (); + return ComplexColumnVector (dup (data (), len), len); +} + +RowVector +real (const ComplexRowVector& a) +{ + int a_len = a.length (); + RowVector retval; + if (a_len > 0) + retval = RowVector (real_dup (a.data (), a_len), a_len); + return retval; +} + +RowVector +imag (const ComplexRowVector& a) +{ + int a_len = a.length (); + RowVector retval; + if (a_len > 0) + retval = RowVector (imag_dup (a.data (), a_len), a_len); + return retval; +} + +ComplexRowVector +conj (const ComplexRowVector& a) +{ + int a_len = a.length (); + ComplexRowVector retval; + if (a_len > 0) + retval = ComplexRowVector (conj_dup (a.data (), a_len), a_len); + return retval; +} + +// resize is the destructive equivalent for this one + +ComplexRowVector +ComplexRowVector::extract (int c1, int c2) const +{ + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + int new_c = c2 - c1 + 1; + + ComplexRowVector result (new_c); + + for (int i = 0; i < new_c; i++) + result.elem (i) = elem (c1+i); + + return result; +} + +// row vector by row vector -> row vector operations + +ComplexRowVector& +ComplexRowVector::operator += (const RowVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ComplexRowVector& +ComplexRowVector::operator -= (const RowVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +ComplexRowVector& +ComplexRowVector::operator += (const ComplexRowVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ComplexRowVector& +ComplexRowVector::operator -= (const ComplexRowVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +// row vector by scalar -> row vector operations + +ComplexRowVector +operator + (const ComplexRowVector& v, double s) +{ + int len = v.length (); + return ComplexRowVector (add (v.data (), len, s), len); +} + +ComplexRowVector +operator - (const ComplexRowVector& v, double s) +{ + int len = v.length (); + return ComplexRowVector (subtract (v.data (), len, s), len); +} + +ComplexRowVector +operator * (const ComplexRowVector& v, double s) +{ + int len = v.length (); + return ComplexRowVector (multiply (v.data (), len, s), len); +} + +ComplexRowVector +operator / (const ComplexRowVector& v, double s) +{ + int len = v.length (); + return ComplexRowVector (divide (v.data (), len, s), len); +} + +// scalar by row vector -> row vector operations + +ComplexRowVector +operator + (double s, const ComplexRowVector& a) +{ + int a_len = a.length (); + return ComplexRowVector (add (a.data (), a_len, s), a_len); +} + +ComplexRowVector +operator - (double s, const ComplexRowVector& a) +{ + int a_len = a.length (); + return ComplexRowVector (subtract (s, a.data (), a_len), a_len); +} + +ComplexRowVector +operator * (double s, const ComplexRowVector& a) +{ + int a_len = a.length (); + return ComplexRowVector (multiply (a.data (), a_len, s), a_len); +} + +ComplexRowVector +operator / (double s, const ComplexRowVector& a) +{ + int a_len = a.length (); + return ComplexRowVector (divide (s, a.data (), a_len), a_len); +} + +// row vector by column vector -> scalar + +Complex +operator * (const ComplexRowVector& v, const ColumnVector& a) +{ + ComplexColumnVector tmp (a); + return v * tmp; +} + +Complex +operator * (const ComplexRowVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return 0.0; + } + + Complex retval (0.0, 0.0); + + for (int i = 0; i < len; i++) + retval += v.elem (i) * a.elem (i); + + return retval; +} + +// row vector by matrix -> row vector + +ComplexRowVector +operator * (const ComplexRowVector& v, const ComplexMatrix& a) +{ + int len = v.length (); + if (a.rows () != len) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return ComplexRowVector (); + } + + if (len == 0 || a.cols () == 0) + return ComplexRowVector (0); + +// Transpose A to form A'*x == (x'*A)' + + int a_nr = a.rows (); + int a_nc = a.cols (); + + char trans = 'T'; + int ld = a_nr; + Complex alpha (1.0); + Complex beta (0.0); + int i_one = 1; + + Complex *y = new Complex [len]; + + F77_FCN (zgemv) (&trans, &a_nc, &a_nr, &alpha, a.data (), &ld, + v.data (), &i_one, &beta, y, &i_one, 1L); + + return ComplexRowVector (y, len); +} + +// row vector by row vector -> row vector operations + +ComplexRowVector +operator + (const ComplexRowVector& v, const RowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector addition attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (add (v.data (), a.data (), len), len); +} + +ComplexRowVector +operator - (const ComplexRowVector& v, const RowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (subtract (v.data (), a.data (), len), len); +} + +ComplexRowVector +product (const ComplexRowVector& v, const RowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector product attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (multiply (v.data (), a.data (), len), len); +} + +ComplexRowVector +quotient (const ComplexRowVector& v, const RowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector quotient attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (divide (v.data (), a.data (), len), len); +} + +// other operations + +ComplexRowVector +map (c_c_Mapper f, const ComplexRowVector& a) +{ + ComplexRowVector b (a); + b.map (f); + return b; +} + +RowVector +map (d_c_Mapper f, const ComplexRowVector& a) +{ + int a_len = a.length (); + RowVector b (a_len); + for (int i = 0; i < a_len; i++) + b.elem (i) = f (a.elem (i)); + return b; +} + +void +ComplexRowVector::map (c_c_Mapper f) +{ + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); +} + +Complex +ComplexRowVector::min (void) const +{ + int len = length (); + if (len == 0) + return Complex (0.0); + + Complex res = elem (0); + double absres = abs (res); + + for (int i = 1; i < len; i++) + if (abs (elem (i)) < absres) + { + res = elem (i); + absres = abs (res); + } + + return res; +} + +Complex +ComplexRowVector::max (void) const +{ + int len = length (); + if (len == 0) + return Complex (0.0); + + Complex res = elem (0); + double absres = abs (res); + + for (int i = 1; i < len; i++) + if (abs (elem (i)) > absres) + { + res = elem (i); + absres = abs (res); + } + + return res; +} + +// i/o + +ostream& +operator << (ostream& os, const ComplexRowVector& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.length (); i++) + os << " " /* setw (field_width) */ << a.elem (i); + return os; +} + +istream& +operator >> (istream& is, ComplexRowVector& a) +{ + int len = a.length(); + + if (len < 1) + is.clear (ios::badbit); + else + { + Complex tmp; + for (int i = 0; i < len; i++) + { + is >> tmp; + if (is) + a.elem (i) = tmp; + else + break; + } + } +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/CRowVector.h b/liboctave/CRowVector.h new file mode 100644 --- /dev/null +++ b/liboctave/CRowVector.h @@ -0,0 +1,169 @@ +// -*- 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. + +*/ + +#if !defined (octave_ComplexRowVector_h) +#define octave_ComplexRowVector_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "mx-defs.h" + +extern "C++" { + +class ComplexRowVector : public Array +{ +friend class ComplexColumnVector; + +public: + + ComplexRowVector (void) : Array () { } + ComplexRowVector (int n) : Array (n) { } + ComplexRowVector (int n, const Complex& val) : Array (n, val) { } + ComplexRowVector (const RowVector& a); + ComplexRowVector (const Array& a) : Array (a) { } + ComplexRowVector (const ComplexRowVector& a) : Array (a) { } +// ComplexRowVector (double a) : Array (1, a) { } +// ComplexRowVector (const Complex& a) : Array (1, a) { } + + ComplexRowVector& operator = (const ComplexRowVector& a) + { + Array::operator = (a); + return *this; + } + +// operator Array& () const { return *this; } + + int operator == (const ComplexRowVector& a) const; + int operator != (const ComplexRowVector& a) const; + +// destructive insert/delete/reorder operations + + ComplexRowVector& insert (const RowVector& a, int c); + ComplexRowVector& insert (const ComplexRowVector& a, int c); + + ComplexRowVector& fill (double val); + ComplexRowVector& fill (const Complex& val); + ComplexRowVector& fill (double val, int c1, int c2); + ComplexRowVector& fill (const Complex& val, int c1, int c2); + + ComplexRowVector append (const RowVector& a) const; + ComplexRowVector append (const ComplexRowVector& a) const; + + ComplexColumnVector hermitian (void) const; // complex conjugate transpose. + ComplexColumnVector transpose (void) const; + + friend RowVector real (const ComplexRowVector& a); + friend RowVector imag (const ComplexRowVector& a); + friend ComplexRowVector conj (const ComplexRowVector& a); + +// resize is the destructive equivalent for this one + + ComplexRowVector extract (int c1, int c2) const; + +// row vector by row vector -> row vector operations + + ComplexRowVector& operator += (const RowVector& a); + ComplexRowVector& operator -= (const RowVector& a); + + ComplexRowVector& operator += (const ComplexRowVector& a); + ComplexRowVector& operator -= (const ComplexRowVector& a); + +// row vector by scalar -> row vector operations + + friend ComplexRowVector operator + (const ComplexRowVector& a, double s); + friend ComplexRowVector operator - (const ComplexRowVector& a, double s); + friend ComplexRowVector operator * (const ComplexRowVector& a, double s); + friend ComplexRowVector operator / (const ComplexRowVector& a, double s); + +// scalar by row vector -> row vector operations + + friend ComplexRowVector operator + (double s, const ComplexRowVector& a); + friend ComplexRowVector operator - (double s, const ComplexRowVector& a); + friend ComplexRowVector operator * (double s, const ComplexRowVector& a); + friend ComplexRowVector operator / (double s, const ComplexRowVector& a); + +// row vector by column vector -> scalar + + friend Complex operator * (const ComplexRowVector& a, const ColumnVector& b); + + friend Complex operator * (const ComplexRowVector& a, + const ComplexColumnVector& b); + +// row vector by matrix -> row vector + + friend ComplexRowVector operator * (const ComplexRowVector& a, + const ComplexMatrix& b); + +// row vector by row vector -> row vector operations + + friend ComplexRowVector operator + (const ComplexRowVector& a, + const RowVector& b); + friend ComplexRowVector operator - (const ComplexRowVector& a, + const RowVector& b); + + friend ComplexRowVector product (const ComplexRowVector& a, + const RowVector& b); + friend ComplexRowVector quotient (const ComplexRowVector& a, + const RowVector& b); + +// other operations + + friend ComplexRowVector map (c_c_Mapper f, const ComplexRowVector& a); + friend RowVector map (d_c_Mapper f, const ComplexRowVector& a); + void map (c_c_Mapper f); + + Complex min (void) const; + Complex max (void) const; + +// i/o + + friend ostream& operator << (ostream& os, const ComplexRowVector& a); + friend ostream& operator >> (ostream& is, ComplexRowVector& a); + +#define KLUDGE_VECTORS +#define TYPE Complex +#define KL_VEC_TYPE ComplexRowVector +#include "mx-kludge.h" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +private: + + ComplexRowVector (Complex *d, int l) : Array (d, l) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dColVector.cc b/liboctave/dColVector.cc new file mode 100644 --- /dev/null +++ b/liboctave/dColVector.cc @@ -0,0 +1,491 @@ +// ColumnVector manipulations. -*- 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 + +#include + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "f77-uscore.h" +#include "lo-error.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (dgemm) (const char*, const char*, const int*, + const int*, const int*, const double*, + const double*, const int*, const double*, + const int*, const double*, double*, const int*, + long, long); +} + +/* + * Column Vector class. + */ + +#define KLUDGE_VECTORS +#define TYPE double +#define KL_VEC_TYPE ColumnVector +#include "mx-kludge.cc" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +#if 0 +ColumnVector& +ColumnVector::resize (int n) +{ + if (n < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimension"); + return *this; + } + + double *new_data = (double *) NULL; + if (n > 0) + { + new_data = new double [n]; + int min_len = len < n ? len : n; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + len = n; + data = new_data; + + return *this; +} + +ColumnVector& +ColumnVector::resize (int n, double val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} +#endif + +int +ColumnVector::operator == (const ColumnVector& a) const +{ + int len = length (); + if (len != a.length ()) + return 0; + return equal (data (), a.data (), len); +} + +int +ColumnVector::operator != (const ColumnVector& a) const +{ + return !(*this == a); +} + +ColumnVector& +ColumnVector::insert (const ColumnVector& a, int r) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); + + return *this; +} + +ColumnVector& +ColumnVector::fill (double val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +ColumnVector& +ColumnVector::fill (double val, int r1, int r2) +{ + int len = length (); + if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + for (int i = r1; i <= r2; i++) + elem (i) = val; + + return *this; +} + +ColumnVector +ColumnVector::stack (const ColumnVector& a) const +{ + int len = length (); + int nr_insert = len; + ColumnVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nr_insert); + return retval; +} + +RowVector +ColumnVector::transpose (void) const +{ + int len = length (); + return RowVector (dup (data (), len), len); +} + +// resize is the destructive equivalent for this one + +ColumnVector +ColumnVector::extract (int r1, int r2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + int new_r = r2 - r1 + 1; + + ColumnVector result (new_r); + + for (int i = 0; i < new_r; i++) + result.elem (i) = elem (r1+i); + + return result; +} + +// column vector by column vector -> column vector operations + +ColumnVector& +ColumnVector::operator += (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return ColumnVector (); + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ColumnVector& +ColumnVector::operator -= (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return ColumnVector (); + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +// scalar by column vector -> column vector operations + +ComplexColumnVector +operator + (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (add (a.data (), len, s), len); +} + +ComplexColumnVector +operator - (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (subtract (a.data (), len, s), len); +} + +ComplexColumnVector +operator * (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (multiply (a.data (), len, s), len); +} + +ComplexColumnVector +operator / (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (divide (a.data (), len, s), len); +} + +// scalar by column vector -> column vector operations + +ComplexColumnVector +operator + (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (add (a.data (), a_len, s), a_len); +} + +ComplexColumnVector +operator - (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (subtract (s, a.data (), a_len), a_len); +} + +ComplexColumnVector +operator * (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (multiply (a.data (), a_len, s), a_len); +} + +ComplexColumnVector +operator / (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (divide (s, a.data (), a_len), a_len); +} + +// column vector by row vector -> matrix operations + +Matrix +operator * (const ColumnVector& v, const RowVector& a) +{ + int len = v.length (); + int a_len = a.length (); + if (len != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return Matrix (); + } + + if (len == 0) + return Matrix (len, len, 0.0); + + char transa = 'N'; + char transb = 'N'; + double alpha = 1.0; + double beta = 0.0; + int anr = 1; + + double *c = new double [len * a_len]; + + F77_FCN (dgemm) (&transa, &transb, &len, &a_len, &anr, &alpha, + v.data (), &len, a.data (), &anr, &beta, c, &len, + 1L, 1L); + + return Matrix (c, len, a_len); +} + +ComplexMatrix +operator * (const ColumnVector& v, const ComplexRowVector& a) +{ + ComplexColumnVector tmp (v); + return tmp * a; +} + +ComplexColumnVector +operator + (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (add (v.data (), a.data (), len), len); +} + +ComplexColumnVector +operator - (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (subtract (v.data (), a.data (), len), len); +} + +ComplexColumnVector +product (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector product attempted"); + return ColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (multiply (v.data (), a.data (), len), len); +} + +ComplexColumnVector +quotient (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector quotient attempted"); + return ColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (divide (v.data (), a.data (), len), len); +} + +// other operations + +ColumnVector +map (d_d_Mapper f, const ColumnVector& a) +{ + ColumnVector b (a); + b.map (f); + return b; +} + +void +ColumnVector::map (d_d_Mapper f) +{ + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); +} + +double +ColumnVector::min (void) const +{ + int len = length (); + if (len == 0) + return 0.0; + + double res = elem (0); + + for (int i = 1; i < len; i++) + if (elem (i) < res) + res = elem (i); + + return res; +} + +double +ColumnVector::max (void) const +{ + int len = length (); + if (len == 0) + return 0.0; + + double res = elem (0); + + for (int i = 1; i < len; i++) + if (elem (i) > res) + res = elem (i); + + return res; +} + +ostream& +operator << (ostream& os, const ColumnVector& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.length (); i++) + os << /* setw (field_width) << */ a.elem (i) << "\n"; + return os; +} + +istream& +operator >> (istream& is, ColumnVector& a) +{ + int len = a.length(); + + if (len < 1) + is.clear (ios::badbit); + else + { + double tmp; + for (int i = 0; i < len; i++) + { + is >> tmp; + if (is) + a.elem (i) = tmp; + else + break; + } + } +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dColVector.h b/liboctave/dColVector.h new file mode 100644 --- /dev/null +++ b/liboctave/dColVector.h @@ -0,0 +1,159 @@ +// -*- 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. + +*/ + +#if !defined (octave_ColumnVector_h) +#define octave_ColumnVector_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "mx-defs.h" + +extern "C++" { + +class ColumnVector : public Array +{ +friend class RowVector; + +public: + + ColumnVector (void) : Array () { } + ColumnVector (int n) : Array (n) { } + ColumnVector (int n, double val) : Array (n, val) { } + ColumnVector (const Array& a) : Array (a) { } + ColumnVector (const ColumnVector& a) : Array (a) { } +// ColumnVector (double a) : Array (1, a) { } + + ColumnVector& operator = (const ColumnVector& a) + { + Array::operator = (a); + return *this; + } + +// operator Array& () const { return *this; } + + int operator == (const ColumnVector& a) const; + int operator != (const ColumnVector& a) const; + +// destructive insert/delete/reorder operations + + ColumnVector& insert (const ColumnVector& a, int r); + + ColumnVector& fill (double val); + ColumnVector& fill (double val, int r1, int r2); + + ColumnVector stack (const ColumnVector& a) const; + + RowVector transpose (void) const; + +// resize is the destructive equivalent for this one + + ColumnVector extract (int r1, int r2) const; + +// column vector by column vector -> column vector operations + + ColumnVector& operator += (const ColumnVector& a); + ColumnVector& operator -= (const ColumnVector& a); + +// column vector by scalar -> column vector operations + + friend ComplexColumnVector operator + (const ColumnVector& a, + const Complex& s); + friend ComplexColumnVector operator - (const ColumnVector& a, + const Complex& s); + friend ComplexColumnVector operator * (const ColumnVector& a, + const Complex& s); + friend ComplexColumnVector operator / (const ColumnVector& a, + const Complex& s); + +// scalar by column vector -> column vector operations + + friend ComplexColumnVector operator + (const Complex& s, + const ColumnVector& a); + friend ComplexColumnVector operator - (const Complex& s, + const ColumnVector& a); + friend ComplexColumnVector operator * (const Complex& s, + const ColumnVector& a); + friend ComplexColumnVector operator / (const Complex& s, + const ColumnVector& a); + +// column vector by row vector -> matrix operations + + friend Matrix operator * (const ColumnVector& a, const RowVector& a); + + friend ComplexMatrix operator * (const ColumnVector& a, + const ComplexRowVector& b); + +// column vector by column vector -> column vector operations + + friend ComplexColumnVector operator + (const ComplexColumnVector& a, + const ComplexColumnVector& b); + + friend ComplexColumnVector operator - (const ComplexColumnVector& a, + const ComplexColumnVector& b); + + friend ComplexColumnVector product (const ComplexColumnVector& a, + const ComplexColumnVector& b); + + friend ComplexColumnVector quotient (const ComplexColumnVector& a, + const ComplexColumnVector& b); + +// other operations + + friend ColumnVector map (d_d_Mapper f, const ColumnVector& a); + void map (d_d_Mapper f); + + double min (void) const; + double max (void) const; + +// i/o + + friend ostream& operator << (ostream& os, const ColumnVector& a); + friend ostream& operator >> (ostream& is, ColumnVector& a); + +#define KLUDGE_VECTORS +#define TYPE double +#define KL_VEC_TYPE ColumnVector +#include "mx-kludge.h" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +private: + + ColumnVector (double *d, int l) : Array (d, l) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dDiagMatrix.cc b/liboctave/dDiagMatrix.cc new file mode 100644 --- /dev/null +++ b/liboctave/dDiagMatrix.cc @@ -0,0 +1,943 @@ +// DiagMatrix manipulations. -*- 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 + +#include + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "lo-error.h" + +/* + * Diagonal Matrix class. + */ + +#define KLUDGE_DIAG_MATRICES +#define TYPE double +#define KL_DMAT_TYPE DiagMatrix +#include "mx-kludge.cc" +#undef KLUDGE_DIAG_MATRICES +#undef TYPE +#undef KL_DMAT_TYPE + +#if 0 +DiagMatrix& +DiagMatrix::resize (int r, int c) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r < c ? r : c; + double *new_data = (double *) NULL; + if (new_len > 0) + { + new_data = new double [new_len]; + + int min_len = new_len < len ? new_len : len; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} + +DiagMatrix& +DiagMatrix::resize (int r, int c, double val) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r < c ? r : c; + double *new_data = (double *) NULL; + if (new_len > 0) + { + new_data = new double [new_len]; + + int min_len = new_len < len ? new_len : len; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + + for (i = min_len; i < new_len; i++) + new_data[i] = val; + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} +#endif + +int +DiagMatrix::operator == (const DiagMatrix& a) const +{ + if (rows () != a.rows () || cols () != a.cols ()) + return 0; + + return equal (data (), a.data (), length ()); +} + +int +DiagMatrix::operator != (const DiagMatrix& a) const +{ + return !(*this == a); +} + +DiagMatrix& +DiagMatrix::fill (double val) +{ + for (int i = 0; i < length (); i++) + elem (i, i) = val; + return *this; +} + +DiagMatrix& +DiagMatrix::fill (double val, int beg, int end) +{ + if (beg < 0 || end >= length () || end < beg) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = beg; i < end; i++) + elem (i, i) = val; + + return *this; +} + +DiagMatrix& +DiagMatrix::fill (const ColumnVector& a) +{ + int len = length (); + if (a.length () != len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + + return *this; +} + +DiagMatrix& +DiagMatrix::fill (const RowVector& a) +{ + int len = length (); + if (a.length () != len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + + return *this; +} + +DiagMatrix& +DiagMatrix::fill (const ColumnVector& a, int beg) +{ + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + + return *this; +} + +DiagMatrix& +DiagMatrix::fill (const RowVector& a, int beg) +{ + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + + return *this; +} + +DiagMatrix +DiagMatrix::transpose (void) const +{ + return DiagMatrix (dup (data (), length ()), cols (), rows ()); +} + +Matrix +DiagMatrix::extract (int r1, int c1, int r2, int c2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + int new_r = r2 - r1 + 1; + int new_c = c2 - c1 + 1; + + Matrix result (new_r, new_c); + + for (int j = 0; j < new_c; j++) + for (int i = 0; i < new_r; i++) + result.elem (i, j) = elem (r1+i, c1+j); + + return result; +} + +// extract row or column i. + +RowVector +DiagMatrix::row (int i) const +{ + int nr = rows (); + int nc = cols (); + if (i < 0 || i >= nr) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } + + RowVector retval (nc, 0.0); + if (nr <= nc || (nr > nc && i < nc)) + retval.elem (i) = elem (i, i); + + return retval; +} + +RowVector +DiagMatrix::row (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return row (0); + else if (c == 'l' || c == 'L') + return row (rows () - 1); + else + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } +} + +ColumnVector +DiagMatrix::column (int i) const +{ + int nr = rows (); + int nc = cols (); + if (i < 0 || i >= nc) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } + + ColumnVector retval (nr, 0.0); + if (nr >= nc || (nr < nc && i < nr)) + retval.elem (i) = elem (i, i); + + return retval; +} + +ColumnVector +DiagMatrix::column (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return column (0); + else if (c == 'l' || c == 'L') + return column (cols () - 1); + else + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } +} + +DiagMatrix +DiagMatrix::inverse (void) const +{ + int info; + return inverse (info); +} + +DiagMatrix +DiagMatrix::inverse (int &info) const +{ + int nr = rows (); + int nc = cols (); + int len = length (); + if (nr != nc) + { + (*current_liboctave_error_handler) ("inverse requires square matrix"); + return DiagMatrix (); + } + + info = 0; + double *tmp_data = dup (data (), len); + for (int i = 0; i < len; i++) + { + if (elem (i, i) == 0.0) + { + info = -1; + copy (tmp_data, data (), len); // Restore contents. + break; + } + else + { + tmp_data[i] = 1.0 / elem (i, i); + } + } + + return DiagMatrix (tmp_data, nr, nc); +} + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + +DiagMatrix& +DiagMatrix::operator += (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nc == 0 || nr == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + return *this; +} + +DiagMatrix& +DiagMatrix::operator -= (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + return *this; +} + +// diagonal matrix by scalar -> matrix operations + +Matrix +operator + (const DiagMatrix& a, double s) +{ + Matrix tmp (a.rows (), a.cols (), s); + return a + tmp; +} + +Matrix +operator - (const DiagMatrix& a, double s) +{ + Matrix tmp (a.rows (), a.cols (), -s); + return a + tmp; +} + +ComplexMatrix +operator + (const DiagMatrix& a, const Complex& s) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return a + tmp; +} + +ComplexMatrix +operator - (const DiagMatrix& a, const Complex& s) +{ + ComplexMatrix tmp (a.rows (), a.cols (), -s); + return a + tmp; +} + +// diagonal matrix by scalar -> diagonal matrix operations + +ComplexDiagMatrix +operator * (const DiagMatrix& a, const Complex& s) +{ + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexDiagMatrix +operator / (const DiagMatrix& a, const Complex& s) +{ + return ComplexDiagMatrix (divide (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +// scalar by diagonal matrix -> matrix operations + +Matrix +operator + (double s, const DiagMatrix& a) +{ + Matrix tmp (a.rows (), a.cols (), s); + return tmp + a; +} + +Matrix +operator - (double s, const DiagMatrix& a) +{ + Matrix tmp (a.rows (), a.cols (), s); + return tmp - a; +} + +ComplexMatrix +operator + (const Complex& s, const DiagMatrix& a) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return tmp + a; +} + +ComplexMatrix +operator - (const Complex& s, const DiagMatrix& a) +{ + ComplexMatrix tmp (a.rows (), a.cols (), s); + return tmp - a; +} + +// scalar by diagonal matrix -> diagonal matrix operations + +ComplexDiagMatrix +operator * (const Complex& s, const DiagMatrix& a) +{ + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +// diagonal matrix by column vector -> column vector operations + +ColumnVector +operator * (const DiagMatrix& m, const ColumnVector& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_len = a.length (); + if (nc != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ColumnVector (); + } + + if (nc == 0 || nr == 0) + return ColumnVector (0); + + ColumnVector result (nr); + + for (int i = 0; i < a_len; i++) + result.elem (i) = a.elem (i) * m.elem (i, i); + + for (i = a_len; i < nr; i++) + result.elem (i) = 0.0; + + return result; +} + +ComplexColumnVector +operator * (const DiagMatrix& m, const ComplexColumnVector& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_len = a.length (); + if (nc != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ColumnVector (); + } + + if (nc == 0 || nr == 0) + return ComplexColumnVector (0); + + ComplexColumnVector result (nr); + + for (int i = 0; i < a_len; i++) + result.elem (i) = a.elem (i) * m.elem (i, i); + + for (i = a_len; i < nr; i++) + result.elem (i) = 0.0; + + return result; +} + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + +DiagMatrix +operator * (const DiagMatrix& a, const DiagMatrix& b) +{ + int nr_a = a.rows (); + int nc_a = a.cols (); + int nr_b = b.rows (); + int nc_b = b.cols (); + if (nc_a != nr_b) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return DiagMatrix (); + } + + if (nr_a == 0 || nc_a == 0 || nc_b == 0) + return DiagMatrix (nr_a, nc_a, 0.0); + + DiagMatrix c (nr_a, nc_b); + + int len = nr_a < nc_b ? nr_a : nc_b; + + for (int i = 0; i < len; i++) + { + double a_element = a.elem (i, i); + double b_element = b.elem (i, i); + + if (a_element == 0.0 || b_element == 0.0) + c.elem (i, i) = 0.0; + else if (a_element == 1.0) + c.elem (i, i) = b_element; + else if (b_element == 1.0) + c.elem (i, i) = a_element; + else + c.elem (i, i) = a_element * b_element; + } + + return c; +} + +ComplexDiagMatrix +operator + (const DiagMatrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexDiagMatrix (); + } + + if (nc == 0 || nr == 0) + return ComplexDiagMatrix (nr, nc); + + return ComplexDiagMatrix (add (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexDiagMatrix +operator - (const DiagMatrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexDiagMatrix (); + } + + if (nc == 0 || nr == 0) + return ComplexDiagMatrix (nr, nc); + + return ComplexDiagMatrix (subtract (m.data (), a.data (), m.length ()), + nr, nc); +} + +ComplexDiagMatrix +operator * (const DiagMatrix& a, const ComplexDiagMatrix& b) +{ + int nr_a = a.rows (); + int nc_a = a.cols (); + int nr_b = b.rows (); + int nc_b = b.cols (); + if (nc_a != nr_b) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexDiagMatrix (); + } + + if (nr_a == 0 || nc_a == 0 || nc_b == 0) + return ComplexDiagMatrix (nr_a, nc_a, 0.0); + + ComplexDiagMatrix c (nr_a, nc_b); + + int len = nr_a < nc_b ? nr_a : nc_b; + + for (int i = 0; i < len; i++) + { + double a_element = a.elem (i, i); + Complex b_element = b.elem (i, i); + + if (a_element == 0.0 || b_element == 0.0) + c.elem (i, i) = 0.0; + else if (a_element == 1.0) + c.elem (i, i) = b_element; + else if (b_element == 1.0) + c.elem (i, i) = a_element; + else + c.elem (i, i) = a_element * b_element; + } + + return c; +} + +ComplexDiagMatrix +product (const DiagMatrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix product attempted"); + return ComplexDiagMatrix (); + } + + if (nc == 0 || nr == 0) + return ComplexDiagMatrix (nr, nc); + + return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()), + nr, nc); +} + +// diagonal matrix by matrix -> matrix operations + +Matrix +operator + (const DiagMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0) + return Matrix (nr, nc); + + Matrix result (a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +Matrix +operator - (const DiagMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0) + return Matrix (nr, nc); + + Matrix result (-a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +Matrix +operator * (const DiagMatrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return Matrix (nr, a_nc, 0.0); + + Matrix c (nr, a_nc); + + for (int i = 0; i < m.length (); i++) + { + if (m.elem (i, i) == 1.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = a.elem (i, j); + } + else if (m.elem (i, i) == 0.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = 0.0; + } + else + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = m.elem (i, i) * a.elem (i, j); + } + } + + if (nr > nc) + { + for (int j = 0; j < a_nc; j++) + for (int i = a_nr; i < nr; i++) + c.elem (i, j) = 0.0; + } + + return c; +} + +ComplexMatrix +operator + (const DiagMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +ComplexMatrix +operator - (const DiagMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (-a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); + + return result; +} + +ComplexMatrix +operator * (const DiagMatrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, nc, 0.0); + + ComplexMatrix c (nr, a_nc); + + for (int i = 0; i < m.length (); i++) + { + if (m.elem (i, i) == 1.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = a.elem (i, j); + } + else if (m.elem (i, i) == 0.0) + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = 0.0; + } + else + { + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = m.elem (i, i) * a.elem (i, j); + } + } + + if (nr > nc) + { + for (int j = 0; j < a_nc; j++) + for (int i = a_nr; i < nr; i++) + c.elem (i, j) = 0.0; + } + + return c; +} + +// other operations + +ColumnVector +DiagMatrix::diag (void) const +{ + return diag (0); +} + +// Could be optimized... + +ColumnVector +DiagMatrix::diag (int k) const +{ + int nnr = rows (); + int nnc = cols (); + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + ColumnVector d; + + if (nnr > 0 && nnc > 0) + { + int ndiag = (nnr < nnc) ? nnr : nnc; + + d.resize (ndiag); + + if (k > 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i+k); + } + else if ( k < 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i-k, i); + } + else + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i); + } + } + else + cerr << "diag: requested diagonal out of range\n"; + + return d; +} + +ostream& +operator << (ostream& os, const DiagMatrix& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.rows (); i++) + { + for (int j = 0; j < a.cols (); j++) + { + if (i == j) + os << " " /* setw (field_width) */ << a.elem (i, i); + else + os << " " /* setw (field_width) */ << 0.0; + } + os << "\n"; + } + return os; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dDiagMatrix.h b/liboctave/dDiagMatrix.h new file mode 100644 --- /dev/null +++ b/liboctave/dDiagMatrix.h @@ -0,0 +1,189 @@ +// -*- 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. + +*/ + +#if !defined (octave_DiagMatrix_h) +#define octave_DiagMatrix_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "dRowVector.h" +#include "dColVector.h" + +#include "mx-defs.h" + +extern "C++" { + +class DiagMatrix : public DiagArray +{ +friend class SVD; +friend class ComplexSVD; + +public: + + DiagMatrix (void) : DiagArray () { } + DiagMatrix (int n) : DiagArray (n) { } + DiagMatrix (int n, double val) : DiagArray (n, val) { } + DiagMatrix (int r, int c) : DiagArray (r, c) { } + DiagMatrix (int r, int c, double val) : DiagArray (r, c, val) { } + DiagMatrix (const RowVector& a) : DiagArray (a) { } + DiagMatrix (const ColumnVector& a) : DiagArray (a) { } + DiagMatrix (const DiagArray& a) : DiagArray (a) { } + DiagMatrix (const DiagMatrix& a) : DiagArray (a) { } +// DiagMatrix (double a) : DiagArray (1, a) { } + + DiagMatrix& operator = (const DiagMatrix& a) + { + DiagArray::operator = (a); + return *this; + } + +// operator DiagArray& () const { return *this; } + + int operator == (const DiagMatrix& a) const; + int operator != (const DiagMatrix& a) const; + + DiagMatrix& fill (double val); + DiagMatrix& fill (double val, int beg, int end); + DiagMatrix& fill (const ColumnVector& a); + DiagMatrix& fill (const RowVector& a); + DiagMatrix& fill (const ColumnVector& a, int beg); + DiagMatrix& fill (const RowVector& a, int beg); + + DiagMatrix transpose (void) const; + +// resize is the destructive analog for this one + + Matrix extract (int r1, int c1, int r2, int c2) const; + +// extract row or column i. + + RowVector row (int i) const; + RowVector row (char *s) const; + + ColumnVector column (int i) const; + ColumnVector column (char *s) const; + + DiagMatrix inverse (void) const; + DiagMatrix inverse (int& info) const; + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + + DiagMatrix& operator += (const DiagMatrix& a); + DiagMatrix& operator -= (const DiagMatrix& a); + +// diagonal matrix by scalar -> matrix operations + + friend Matrix operator + (const DiagMatrix& a, double s); + friend Matrix operator - (const DiagMatrix& a, double s); + + friend ComplexMatrix operator + (const DiagMatrix& a, const Complex& s); + friend ComplexMatrix operator - (const DiagMatrix& a, const Complex& s); + +// diagonal matrix by scalar -> diagonal matrix operations + + friend ComplexDiagMatrix operator * (const DiagMatrix& a, const Complex& s); + friend ComplexDiagMatrix operator / (const DiagMatrix& a, const Complex& s); + +// scalar by diagonal matrix -> matrix operations + + friend Matrix operator + (double s, const DiagMatrix& a); + friend Matrix operator - (double s, const DiagMatrix& a); + + friend ComplexMatrix operator + (const Complex& s, const DiagMatrix& a); + friend ComplexMatrix operator - (const Complex& s, const DiagMatrix& a); + +// scalar by diagonal matrix -> diagonal matrix operations + + friend ComplexDiagMatrix operator * (const Complex& s, const DiagMatrix& a); + +// diagonal matrix by column vector -> column vector operations + + friend ColumnVector operator * (const DiagMatrix& a, const ColumnVector& b); + + friend ComplexColumnVector operator * (const DiagMatrix& a, + const ComplexColumnVector& b); + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + + friend DiagMatrix operator * (const DiagMatrix& a, + const DiagMatrix& b); + + friend ComplexDiagMatrix operator + (const DiagMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexDiagMatrix operator - (const DiagMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexDiagMatrix operator * (const DiagMatrix& a, + const ComplexDiagMatrix& b); + + friend ComplexDiagMatrix product (const DiagMatrix& a, + const ComplexDiagMatrix& b); + +// diagonal matrix by matrix -> matrix operations + + friend Matrix operator + (const DiagMatrix& a, const Matrix& b); + friend Matrix operator - (const DiagMatrix& a, const Matrix& b); + friend Matrix operator * (const DiagMatrix& a, const Matrix& b); + + friend ComplexMatrix operator + (const DiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator - (const DiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator * (const DiagMatrix& a, + const ComplexMatrix& b); + +// other operations + + ColumnVector diag (void) const; + ColumnVector diag (int k) const; + +// i/o + + friend ostream& operator << (ostream& os, const DiagMatrix& a); + +#define KLUDGE_DIAG_MATRICES +#define TYPE double +#define KL_DMAT_TYPE DiagMatrix +#include "mx-kludge.h" +#undef KLUDGE_DIAG_MATRICES +#undef TYPE +#undef KL_DMAT_TYPE + +private: + + DiagMatrix (double *d, int nr, int nc) : DiagArray (d, nr, nc) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc new file mode 100644 --- /dev/null +++ b/liboctave/dMatrix.cc @@ -0,0 +1,2432 @@ +// Matrix manipulations. -*- 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 +#include +#include + +#include + +#include "mx-base.h" +#include "dbleDET.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (dgemm) (const char*, const char*, const int*, + const int*, const int*, const double*, + const double*, const int*, const double*, + const int*, const double*, double*, const int*, + long, long); + + int F77_FCN (dgemv) (const char*, const int*, const int*, + const double*, const double*, const int*, + const double*, const int*, const double*, + double*, const int*, long); + + int F77_FCN (dgeco) (double*, const int*, const int*, int*, double*, + double*); + + int F77_FCN (dgesl) (const double*, const int*, const int*, + const int*, double*, const int*); + + int F77_FCN (dgedi) (double*, const int*, const int*, const int*, + double*, double*, const int*); + + int F77_FCN (dgelss) (const int*, const int*, const int*, double*, + const int*, double*, const int*, double*, + const double*, int*, double*, const int*, + int*); + +// Note that the original complex fft routines were not written for +// double complex arguments. They have been modified by adding an +// implicit double precision (a-h,o-z) statement at the beginning of +// each subroutine. + + int F77_FCN (cffti) (const int*, Complex*); + + int F77_FCN (cfftf) (const int*, Complex*, Complex*); + + int F77_FCN (cfftb) (const int*, Complex*, Complex*); +} + +#define KLUDGE_MATRICES +#define TYPE double +#define KL_MAT_TYPE Matrix +#include "mx-kludge.cc" +#undef KLUDGE_MATRICES +#undef TYPE +#undef KL_MAT_TYPE + +/* + * Matrix class. + */ + +Matrix::Matrix (const DiagMatrix& a) + : Array2 (a.rows (), a.cols (), 0.0) +{ + for (int i = 0; i < a.length (); i++) + elem (i, i) = a.elem (i, i); +} + +#if 0 +Matrix& +Matrix::resize (int r, int c) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r * c; + double* new_data = (double *) NULL; + if (new_len > 0) + { + new_data = new double [new_len]; + + int min_r = nr < r ? nr : r; + int min_c = nc < c ? nc : c; + + for (int j = 0; j < min_c; j++) + for (int i = 0; i < min_r; i++) + new_data[r*j+i] = elem (i, j); + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} + +Matrix& +Matrix::resize (int r, int c, double val) +{ + if (r < 0 || c < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimensions"); + return *this; + } + + int new_len = r * c; + double *new_data = (double *) NULL; + if (new_len > 0) + { + new_data = new double [new_len]; + +// There may be faster or cleaner ways to do this. + + if (r > nr || c > nc) + copy (new_data, new_len, val); + + int min_r = nr < r ? nr : r; + int min_c = nc < c ? nc : c; + + for (int j = 0; j < min_c; j++) + for (int i = 0; i < min_r; i++) + new_data[r*j+i] = elem (i, j); + } + + delete [] data; + nr = r; + nc = c; + len = new_len; + data = new_data; + + return *this; +} +#endif + +int +Matrix::operator == (const Matrix& a) const +{ + if (rows () != a.rows () || cols () != a.cols ()) + return 0; + + return equal (data (), a.data (), length ()); +} + +int +Matrix::operator != (const Matrix& a) const +{ + return !(*this == a); +} + +Matrix& +Matrix::insert (const Matrix& a, int r, int c) +{ + int a_rows = a.rows (); + int a_cols = a.cols (); + if (r < 0 || r + a_rows - 1 > rows () + || c < 0 || c + a_cols - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int j = 0; j < a_cols; j++) + for (int i = 0; i < a_rows; i++) + elem (r+i, c+j) = a.elem (i, j); + + return *this; +} + +Matrix& +Matrix::insert (const RowVector& a, int r, int c) +{ + int a_len = a.length (); + if (r < 0 || r >= rows () || c < 0 || c + a_len - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r, c+i) = a.elem (i); + + return *this; +} + +Matrix& +Matrix::insert (const ColumnVector& a, int r, int c) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > rows () || c < 0 || c >= cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i, c) = a.elem (i); + + return *this; +} + +Matrix& +Matrix::insert (const DiagMatrix& a, int r, int c) +{ + if (r < 0 || r + a.rows () - 1 > rows () + || c < 0 || c + a.cols () - 1 > cols ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a.length (); i++) + elem (r+i, c+i) = a.elem (i, i); + + return *this; +} + +Matrix& +Matrix::fill (double val) +{ + int nr = rows (); + int nc = cols (); + if (nr > 0 && nc > 0) + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + elem (i, j) = val; + + return *this; +} + +Matrix& +Matrix::fill (double val, int r1, int c1, int r2, int c2) +{ + int nr = rows (); + int nc = cols (); + if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 + || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + for (int j = c1; j <= c2; j++) + for (int i = r1; i <= r2; i++) + elem (i, j) = val; + + return *this; +} + +Matrix +Matrix::append (const Matrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return Matrix (); + } + + int nc_insert = nc; + Matrix retval (nr, nc + a.cols ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +Matrix +Matrix::append (const RowVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != 1) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return Matrix (); + } + + int nc_insert = nc; + Matrix retval (nr, nc + a.length ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +Matrix +Matrix::append (const ColumnVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.length ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return Matrix (); + } + + int nc_insert = nc; + Matrix retval (nr, nc + 1); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +Matrix +Matrix::append (const DiagMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows ()) + { + (*current_liboctave_error_handler) ("row dimension mismatch for append"); + return *this; + } + + int nc_insert = nc; + Matrix retval (nr, nc + a.cols ()); + retval.insert (*this, 0, 0); + retval.insert (a, 0, nc_insert); + return retval; +} + +Matrix +Matrix::stack (const Matrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return Matrix (); + } + + int nr_insert = nr; + Matrix retval (nr + a.rows (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +Matrix +Matrix::stack (const RowVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.length ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return Matrix (); + } + + int nr_insert = nr; + Matrix retval (nr + 1, nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +Matrix +Matrix::stack (const ColumnVector& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != 1) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return Matrix (); + } + + int nr_insert = nr; + Matrix retval (nr + a.length (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +Matrix +Matrix::stack (const DiagMatrix& a) const +{ + int nr = rows (); + int nc = cols (); + if (nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("column dimension mismatch for stack"); + return Matrix (); + } + + int nr_insert = nr; + Matrix retval (nr + a.rows (), nc); + retval.insert (*this, 0, 0); + retval.insert (a, nr_insert, 0); + return retval; +} + +Matrix +Matrix::transpose (void) const +{ + int nr = rows (); + int nc = cols (); + Matrix result (nc, nr); + if (length () > 0) + { + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + result.elem (j, i) = elem (i, j); + } + return result; +} + +Matrix +Matrix::extract (int r1, int c1, int r2, int c2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + int new_r = r2 - r1 + 1; + int new_c = c2 - c1 + 1; + + Matrix result (new_r, new_c); + + for (int j = 0; j < new_c; j++) + for (int i = 0; i < new_r; i++) + result.elem (i, j) = elem (r1+i, c1+j); + + return result; +} + +// extract row or column i. + +RowVector +Matrix::row (int i) const +{ + int nc = cols (); + if (i < 0 || i >= rows ()) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } + + RowVector retval (nc); + for (int j = 0; j < nc; j++) + retval.elem (j) = elem (i, j); + + return retval; +} + +RowVector +Matrix::row (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return row (0); + else if (c == 'l' || c == 'L') + return row (rows () - 1); + else + { + (*current_liboctave_error_handler) ("invalid row selection"); + return RowVector (); + } +} + +ColumnVector +Matrix::column (int i) const +{ + int nr = rows (); + if (i < 0 || i >= cols ()) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } + + ColumnVector retval (nr); + for (int j = 0; j < nr; j++) + retval.elem (j) = elem (j, i); + + return retval; +} + +ColumnVector +Matrix::column (char *s) const +{ + if (s == (char *) NULL) + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } + + char c = *s; + if (c == 'f' || c == 'F') + return column (0); + else if (c == 'l' || c == 'L') + return column (cols () - 1); + else + { + (*current_liboctave_error_handler) ("invalid column selection"); + return ColumnVector (); + } +} + +Matrix +Matrix::inverse (void) const +{ + int info; + double rcond; + return inverse (info, rcond); +} + +Matrix +Matrix::inverse (int& info) const +{ + double rcond; + return inverse (info, rcond); +} + +Matrix +Matrix::inverse (int& info, double& rcond) const +{ + int nr = rows (); + int nc = cols (); + int len = length (); + if (nr != nc || nr == 0 || nc == 0) + { + (*current_liboctave_error_handler) ("inverse requires square matrix"); + return Matrix (); + } + + info = 0; + + int *ipvt = new int [nr]; + double *z = new double [nr]; + double *tmp_data = dup (data (), len); + + F77_FCN (dgeco) (tmp_data, &nr, &nc, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -1; + copy (tmp_data, data (), len); // Restore matrix contents. + } + else + { + int job = 1; + double dummy; + + F77_FCN (dgedi) (tmp_data, &nr, &nc, ipvt, &dummy, z, &job); + } + + delete [] ipvt; + delete [] z; + + return Matrix (tmp_data, nr, nc); +} + +ComplexMatrix +Matrix::fourier (void) const +{ + int nr = rows (); + int nc = cols (); + int npts, nsamples; + if (nr == 1 || nc == 1) + { + npts = nr > nc ? nr : nc; + nsamples = 1; + } + else + { + npts = nr; + nsamples = nc; + } + + int nn = 4*npts+15; + Complex *wsave = new Complex [nn]; + Complex *tmp_data = make_complex (data (), length ()); + + F77_FCN (cffti) (&npts, wsave); + + for (int j = 0; j < nsamples; j++) + F77_FCN (cfftf) (&npts, &tmp_data[npts*j], wsave); + + delete [] wsave; + + return ComplexMatrix (tmp_data, nr, nc); +} + +ComplexMatrix +Matrix::ifourier (void) const +{ + int nr = rows (); + int nc = cols (); + int npts, nsamples; + if (nr == 1 || nc == 1) + { + npts = nr > nc ? nr : nc; + nsamples = 1; + } + else + { + npts = nr; + nsamples = nc; + } + + int nn = 4*npts+15; + Complex *wsave = new Complex [nn]; + Complex *tmp_data = make_complex (data (), length ()); + + F77_FCN (cffti) (&npts, wsave); + + for (int j = 0; j < nsamples; j++) + F77_FCN (cfftb) (&npts, &tmp_data[npts*j], wsave); + + for (j = 0; j < npts*nsamples; j++) + tmp_data[j] = tmp_data[j] / (double) npts; + + delete [] wsave; + + return ComplexMatrix (tmp_data, nr, nc); +} + +DET +Matrix::determinant (void) const +{ + int info; + double rcond; + return determinant (info, rcond); +} + +DET +Matrix::determinant (int& info) const +{ + double rcond; + return determinant (info, rcond); +} + +DET +Matrix::determinant (int& info, double& rcond) const +{ + DET retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 0 || nc == 0) + { + double d[2]; + d[0] = 1.0; + d[1] = 0.0; + retval = DET (d); + } + else + { + info = 0; + int *ipvt = new int [nr]; + + double *z = new double [nr]; + double *tmp_data = dup (data (), length ()); + + F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -1; + retval = DET (); + } + else + { + int job = 10; + double d[2]; + F77_FCN (dgedi) (tmp_data, &nr, &nr, ipvt, d, z, &job); + retval = DET (d); + } + + delete [] tmp_data; + delete [] ipvt; + delete [] z; + } + + return retval; +} + +Matrix +Matrix::solve (const Matrix& b) const +{ + int info; + double rcond; + return solve (b, info, rcond); +} + +Matrix +Matrix::solve (const Matrix& b, int& info) const +{ + double rcond; + return solve (b, info, rcond); +} + +Matrix +Matrix::solve (const Matrix& b, int& info, double& rcond) const +{ + Matrix retval; + + int nr = rows (); + int nc = cols (); + if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch solution of linear equations"); + return Matrix (); + } + + info = 0; + int *ipvt = new int [nr]; + + double *z = new double [nr]; + double *tmp_data = dup (data (), length ()); + + F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -2; + } + else + { + int job = 0; + + double *result = dup (b.data (), b.length ()); + + int b_nc = b.cols (); + for (int j = 0; j < b_nc; j++) + F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, &result[nr*j], &job); + + retval = Matrix (result, b.rows (), b_nc); + } + + delete [] tmp_data; + delete [] ipvt; + delete [] z; + + return retval; +} + +ComplexMatrix +Matrix::solve (const ComplexMatrix& b) const +{ + ComplexMatrix tmp (*this); + return tmp.solve (b); +} + +ComplexMatrix +Matrix::solve (const ComplexMatrix& b, int& info) const +{ + ComplexMatrix tmp (*this); + return tmp.solve (b, info); +} + +ComplexMatrix +Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const +{ + ComplexMatrix tmp (*this); + return tmp.solve (b, info, rcond); +} + +ColumnVector +Matrix::solve (const ColumnVector& b) const +{ + int info; double rcond; + return solve (b, info, rcond); +} + +ColumnVector +Matrix::solve (const ColumnVector& b, int& info) const +{ + double rcond; + return solve (b, info, rcond); +} + +ColumnVector +Matrix::solve (const ColumnVector& b, int& info, double& rcond) const +{ + ColumnVector retval; + + int nr = rows (); + int nc = cols (); + if (nr == 0 || nc == 0 || nr != nc || nr != b.length ()) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch solution of linear equations"); + return ColumnVector (); + } + + info = 0; + int *ipvt = new int [nr]; + + double *z = new double [nr]; + double *tmp_data = dup (data (), length ()); + + F77_FCN (dgeco) (tmp_data, &nr, &nr, ipvt, &rcond, z); + + if (rcond + 1.0 == 1.0) + { + info = -2; + } + else + { + int job = 0; + + int b_len = b.length (); + + double *result = dup (b.data (), b_len); + + F77_FCN (dgesl) (tmp_data, &nr, &nr, ipvt, result, &job); + + retval = ColumnVector (result, b_len); + } + + delete [] tmp_data; + delete [] ipvt; + delete [] z; + + return retval; +} + +ComplexColumnVector +Matrix::solve (const ComplexColumnVector& b) const +{ + ComplexMatrix tmp (*this); + return tmp.solve (b); +} + +ComplexColumnVector +Matrix::solve (const ComplexColumnVector& b, int& info) const +{ + ComplexMatrix tmp (*this); + return tmp.solve (b, info); +} + +ComplexColumnVector +Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const +{ + ComplexMatrix tmp (*this); + return tmp.solve (b, info, rcond); +} + +Matrix +Matrix::lssolve (const Matrix& b) const +{ + int info; + int rank; + return lssolve (b, info, rank); +} + +Matrix +Matrix::lssolve (const Matrix& b, int& info) const +{ + int rank; + return lssolve (b, info, rank); +} + +Matrix +Matrix::lssolve (const Matrix& b, int& info, int& rank) const +{ + int nrhs = b.cols (); + + int m = rows (); + int n = cols (); + + if (m == 0 || n == 0 || m != b.rows ()) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of least squares problem"); + return Matrix (); + } + + double *tmp_data = dup (data (), length ()); + + int nrr = m > n ? m : n; + Matrix result (nrr, nrhs); + + int i, j; + for (j = 0; j < nrhs; j++) + for (i = 0; i < m; i++) + result.elem (i, j) = b.elem (i, j); + + double *presult = result.fortran_vec (); + + int len_s = m < n ? m : n; + double *s = new double [len_s]; + double rcond = -1.0; + int lwork; + if (m < n) + lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n)); + else + lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m)); + + double *work = new double [lwork]; + + F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, + &rcond, &rank, work, &lwork, &info); + + Matrix retval (n, nrhs); + for (j = 0; j < nrhs; j++) + for (i = 0; i < n; i++) + retval.elem (i, j) = result.elem (i, j); + + delete [] tmp_data; + delete [] s; + delete [] work; + + return retval; +} + +ComplexMatrix +Matrix::lssolve (const ComplexMatrix& b) const +{ + ComplexMatrix tmp (*this); + return tmp.lssolve (b); +} + +ComplexMatrix +Matrix::lssolve (const ComplexMatrix& b, int& info) const +{ + ComplexMatrix tmp (*this); + return tmp.lssolve (b); +} + +ComplexMatrix +Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const +{ + ComplexMatrix tmp (*this); + return tmp.lssolve (b); +} + +ColumnVector +Matrix::lssolve (const ColumnVector& b) const +{ + int info; + int rank; return lssolve (b, info, rank); +} + +ColumnVector +Matrix::lssolve (const ColumnVector& b, int& info) const +{ + int rank; + return lssolve (b, info, rank); +} + +ColumnVector +Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const +{ + int nrhs = 1; + + int m = rows (); + int n = cols (); + + if (m == 0 || n == 0 || m != b.length ()) + { + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of least squares problem"); + return ColumnVector (); + } + + double *tmp_data = dup (data (), length ()); + + int nrr = m > n ? m : n; + ColumnVector result (nrr); + + int i; + for (i = 0; i < m; i++) + result.elem (i) = b.elem (i); + + double *presult = result.fortran_vec (); + + int len_s = m < n ? m : n; + double *s = new double [len_s]; + double rcond = -1.0; + int lwork; + if (m < n) + lwork = 3*m + (2*m > nrhs ? (2*m > n ? 2*m : n) : (nrhs > n ? nrhs : n)); + else + lwork = 3*n + (2*n > nrhs ? (2*n > m ? 2*n : m) : (nrhs > m ? nrhs : m)); + + double *work = new double [lwork]; + + F77_FCN (dgelss) (&m, &n, &nrhs, tmp_data, &m, presult, &nrr, s, + &rcond, &rank, work, &lwork, &info); + + ColumnVector retval (n); + for (i = 0; i < n; i++) + retval.elem (i) = result.elem (i); + + delete [] tmp_data; + delete [] s; + delete [] work; + + return retval; +} + +ComplexColumnVector +Matrix::lssolve (const ComplexColumnVector& b) const +{ + ComplexMatrix tmp (*this); + return tmp.lssolve (b); +} + +ComplexColumnVector +Matrix::lssolve (const ComplexColumnVector& b, int& info) const +{ + ComplexMatrix tmp (*this); + return tmp.lssolve (b, info); +} + +ComplexColumnVector +Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const +{ + ComplexMatrix tmp (*this); + return tmp.lssolve (b, info, rank); +} + +Matrix& +Matrix::operator += (const Matrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + + return *this; +} + +Matrix& +Matrix::operator -= (const Matrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + + return *this; +} + +Matrix& +Matrix::operator += (const DiagMatrix& a) +{ + if (rows () != a.rows () || cols () != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + for (int i = 0; i < a.length (); i++) + elem (i, i) += a.elem (i, i); + + return *this; +} + +Matrix& +Matrix::operator -= (const DiagMatrix& a) +{ + if (rows () != a.rows () || cols () != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + for (int i = 0; i < a.length (); i++) + elem (i, i) -= a.elem (i, i); + + return *this; +} + +// unary operations + +Matrix +Matrix::operator ! (void) const +{ + int nr = rows (); + int nc = cols (); + + Matrix b (nr, nc); + + for (int j = 0; j < nc; j++) + for (int i = 0; i < nr; i++) + b.elem (i, j) = ! elem (i, j); + + return b; +} + +// matrix by scalar -> matrix operations. + +ComplexMatrix +operator + (const Matrix& a, const Complex& s) +{ + return ComplexMatrix (add (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator - (const Matrix& a, const Complex& s) +{ + return ComplexMatrix (subtract (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator * (const Matrix& a, const Complex& s) +{ + return ComplexMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator / (const Matrix& a, const Complex& s) +{ + return ComplexMatrix (divide (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +// scalar by matrix -> matrix operations. + +ComplexMatrix +operator + (const Complex& s, const Matrix& a) +{ + return ComplexMatrix (add (s, a.data (), a.length ()), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator - (const Complex& s, const Matrix& a) +{ + return ComplexMatrix (subtract (s, a.data (), a.length ()), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator * (const Complex& s, const Matrix& a) +{ + return ComplexMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); +} + +ComplexMatrix +operator / (const Complex& s, const Matrix& a) +{ + return ComplexMatrix (divide (s, a.data (), a.length ()), + a.rows (), a.cols ()); +} + +// matrix by column vector -> column vector operations + +ColumnVector +operator * (const Matrix& m, const ColumnVector& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nc != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ColumnVector (); + } + + if (nr == 0 || nc == 0) + return ColumnVector (0); + + char trans = 'N'; + int ld = nr; + double alpha = 1.0; + double beta = 0.0; + int i_one = 1; + + double *y = new double [nr]; + + F77_FCN (dgemv) (&trans, &nr, &nc, &alpha, m.data (), &ld, a.data (), + &i_one, &beta, y, &i_one, 1L); + + return ColumnVector (y, nr); +} + +ComplexColumnVector +operator * (const Matrix& m, const ComplexColumnVector& a) +{ + ComplexMatrix tmp (m); + return tmp * a; +} + +// matrix by diagonal matrix -> matrix operations + +Matrix +operator + (const Matrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0) + return Matrix (nr, nc); + + Matrix result (m); + int a_len = a.length (); + for (int i = 0; i < a_len; i++) + result.elem (i, i) += a.elem (i, i); + + return result; +} + +Matrix +operator - (const Matrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0) + return Matrix (nr, nc); + + Matrix result (m); + int a_len = a.length (); + for (int i = 0; i < a_len; i++) + result.elem (i, i) -= a.elem (i, i); + + return result; +} + +Matrix +operator * (const Matrix& m, const DiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return Matrix (nr, a_nc, 0.0); + + double *c = new double [nr*a_nc]; + double *ctmp = (double *) NULL; + + int a_len = a.length (); + for (int j = 0; j < a_len; j++) + { + int idx = j * nr; + ctmp = c + idx; + if (a.elem (j, j) == 1.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = m.elem (i, j); + } + else if (a.elem (j, j) == 0.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = 0.0; + } + else + { + for (int i = 0; i < nr; i++) + ctmp[i] = a.elem (j, j) * m.elem (i, j); + } + } + + if (a_nr < a_nc) + { + for (int i = nr * nc; i < nr * a_nc; i++) + ctmp[i] = 0.0; + } + + return Matrix (c, nr, a_nc); +} + +ComplexMatrix +operator + (const Matrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (m); + for (int i = 0; i < a.length (); i++) + result.elem (i, i) += a.elem (i, i); + + return result; +} + +ComplexMatrix +operator - (const Matrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + ComplexMatrix result (m); + for (int i = 0; i < a.length (); i++) + result.elem (i, i) -= a.elem (i, i); + + return result; +} + +ComplexMatrix +operator * (const Matrix& m, const ComplexDiagMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return ComplexMatrix (nr, a_nc, 0.0); + + Complex *c = new Complex [nr*a_nc]; + Complex *ctmp = (Complex *) NULL; + + for (int j = 0; j < a.length (); j++) + { + int idx = j * nr; + ctmp = c + idx; + if (a.elem (j, j) == 1.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = m.elem (i, j); + } + else if (a.elem (j, j) == 0.0) + { + for (int i = 0; i < nr; i++) + ctmp[i] = 0.0; + } + else + { + for (int i = 0; i < nr; i++) + ctmp[i] = a.elem (j, j) * m.elem (i, j); + } + } + + if (a_nr < a_nc) + { + for (int i = nr * nc; i < nr * a_nc; i++) + ctmp[i] = 0.0; + } + + return ComplexMatrix (c, nr, a_nc); +} + +// matrix by matrix -> matrix operations + +Matrix +operator * (const Matrix& m, const Matrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + int a_nr = a.rows (); + int a_nc = a.cols (); + if (nc != a_nr) + { + (*current_liboctave_error_handler) + ("nonconformant matrix multiplication attempted"); + return Matrix (); + } + + if (nr == 0 || nc == 0 || a_nc == 0) + return Matrix (nr, a_nc, 0.0); + + char trans = 'N'; + char transa = 'N'; + + int ld = nr; + int lda = a_nr; + + double alpha = 1.0; + double beta = 0.0; + + double *c = new double [nr*a_nc]; + + F77_FCN (dgemm) (&trans, &transa, &nr, &a_nc, &nc, &alpha, m.data (), + &ld, a.data (), &lda, &beta, c, &nr, 1L, 1L); + + return Matrix (c, nr, a_nc); +} + +ComplexMatrix +operator * (const Matrix& m, const ComplexMatrix& a) +{ + ComplexMatrix tmp (m); + return tmp * a; +} + +ComplexMatrix +operator + (const Matrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix addition attempted"); + return ComplexMatrix (); + } + + return ComplexMatrix (add (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexMatrix +operator - (const Matrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix subtraction attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (subtract (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexMatrix +product (const Matrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix product attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (multiply (m.data (), a.data (), m.length ()), nr, nc); +} + +ComplexMatrix +quotient (const Matrix& m, const ComplexMatrix& a) +{ + int nr = m.rows (); + int nc = m.cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix quotient attempted"); + return ComplexMatrix (); + } + + if (nr == 0 || nc == 0) + return ComplexMatrix (nr, nc); + + return ComplexMatrix (divide (m.data (), a.data (), m.length ()), nr, nc); +} + +// other operations. + +Matrix +map (d_d_Mapper f, const Matrix& a) +{ + Matrix b (a); + b.map (f); + return b; +} + +void +Matrix::map (d_d_Mapper f) +{ + double *d = fortran_vec (); // Ensures only one reference to my privates! + + for (int i = 0; i < length (); i++) + d[i] = f (d[i]); +} + +// XXX FIXME XXX Do these really belong here? They should maybe be +// cleaned up a bit, no? What about corresponding functions for the +// Vectors? + +Matrix +Matrix::all (void) const +{ + int nr = rows (); + int nc = cols (); + Matrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int j = 0; j < nc; j++) + { + if (elem (0, j) == 0.0) + { + retval.elem (0, 0) = 0.0; + break; + } + } + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, 0) == 0.0) + { + retval.elem (0, 0) = 0.0; + break; + } + } + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 1.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, j) == 0.0) + { + retval.elem (0, j) = 0.0; + break; + } + } + } + } + } + return retval; +} + +Matrix +Matrix::any (void) const +{ + int nr = rows (); + int nc = cols (); + Matrix retval; + if (nr > 0 && nc > 0) + { + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int j = 0; j < nc; j++) + { + if (elem (0, j) != 0.0) + { + retval.elem (0, 0) = 1.0; + break; + } + } + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, 0) != 0.0) + { + retval.elem (0, 0) = 1.0; + break; + } + } + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 0.0; + for (int i = 0; i < nr; i++) + { + if (elem (i, j) != 0.0) + { + retval.elem (0, j) = 1.0; + break; + } + } + } + } + } + return retval; +} + +Matrix +Matrix::cumprod (void) const +{ + Matrix retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 1) + { + retval.resize (1, nc); + if (nc > 0) + { + double prod = elem (0, 0); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = prod; + if (j < nc - 1) + prod *= elem (0, j+1); + } + } + } + else if (nc == 1) + { + retval.resize (nr, 1); + if (nr > 0) + { + double prod = elem (0, 0); + for (int i = 0; i < nr; i++) + { + retval.elem (i, 0) = prod; + if (i < nr - 1) + prod *= elem (i+1, 0); + } + } + } + else + { + retval.resize (nr, nc); + if (nr > 0 && nc > 0) + { + for (int j = 0; j < nc; j++) + { + double prod = elem (0, j); + for (int i = 0; i < nr; i++) + { + retval.elem (i, j) = prod; + if (i < nr - 1) + prod *= elem (i+1, j); + } + } + } + } + return retval; +} + +Matrix +Matrix::cumsum (void) const +{ + Matrix retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 1) + { + retval.resize (1, nc); + if (nc > 0) + { + double sum = elem (0, 0); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = sum; + if (j < nc - 1) + sum += elem (0, j+1); + } + } + } + else if (nc == 1) + { + retval.resize (nr, 1); + if (nr > 0) + { + double sum = elem (0, 0); + for (int i = 0; i < nr; i++) + { + retval.elem (i, 0) = sum; + if (i < nr - 1) + sum += elem (i+1, 0); + } + } + } + else + { + retval.resize (nr, nc); + if (nr > 0 && nc > 0) + { + for (int j = 0; j < nc; j++) + { + double sum = elem (0, j); + for (int i = 0; i < nr; i++) + { + retval.elem (i, j) = sum; + if (i < nr - 1) + sum += elem (i+1, j); + } + } + } + } + return retval; +} + +Matrix +Matrix::prod (void) const +{ + Matrix retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int j = 0; j < nc; j++) + retval.elem (0, 0) *= elem (0, j); + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + for (int i = 0; i < nr; i++) + retval.elem (0, 0) *= elem (i, 0); + } + else + { + if (nc == 0) + { + retval.resize (1, 1); + retval.elem (0, 0) = 1.0; + } + else + retval.resize (1, nc); + + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 1.0; + for (int i = 0; i < nr; i++) + retval.elem (0, j) *= elem (i, j); + } + } + return retval; +} + +Matrix +Matrix::sum (void) const +{ + Matrix retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int j = 0; j < nc; j++) + retval.elem (0, 0) += elem (0, j); + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int i = 0; i < nr; i++) + retval.elem (0, 0) += elem (i, 0); + } + else + { + if (nc == 0) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + } + else + retval.resize (1, nc); + + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 0.0; + for (int i = 0; i < nr; i++) + retval.elem (0, j) += elem (i, j); + } + } + return retval; +} + +Matrix +Matrix::sumsq (void) const +{ + Matrix retval; + + int nr = rows (); + int nc = cols (); + + if (nr == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int j = 0; j < nc; j++) + { + double d = elem (0, j); + retval.elem (0, 0) += d * d; + } + } + else if (nc == 1) + { + retval.resize (1, 1); + retval.elem (0, 0) = 0.0; + for (int i = 0; i < nr; i++) + { + double d = elem (i, 0); + retval.elem (0, 0) += d * d; + } + } + else + { + retval.resize (1, nc); + for (int j = 0; j < nc; j++) + { + retval.elem (0, j) = 0.0; + for (int i = 0; i < nr; i++) + { + double d = elem (i, j); + retval.elem (0, j) += d * d; + } + } + } + return retval; +} + +ColumnVector +Matrix::diag (void) const +{ + return diag (0); +} + +ColumnVector +Matrix::diag (int k) const +{ + int nnr = rows (); + int nnc = cols (); + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + ColumnVector d; + + if (nnr > 0 && nnc > 0) + { + int ndiag = (nnr < nnc) ? nnr : nnc; + + d.resize (ndiag); + + if (k > 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i+k); + } + else if ( k < 0) + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i-k, i); + } + else + { + for (int i = 0; i < ndiag; i++) + d.elem (i) = elem (i, i); + } + } + else + cerr << "diag: requested diagonal out of range\n"; + + return d; +} + +ColumnVector +Matrix::row_min (void) const +{ + ColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + double res = elem (i, 0); + for (int j = 1; j < nc; j++) + if (elem (i, j) < res) + res = elem (i, j); + result.elem (i) = res; + } + } + + return result; +} + +ColumnVector +Matrix::row_min_loc (void) const +{ + ColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + int res = 0; + for (int j = 0; j < nc; j++) + if (elem (i, j) < elem (i, res)) + res = j; + result.elem (i) = (double) (res + 1); + } + } + + return result; +} + +ColumnVector +Matrix::row_max (void) const +{ + ColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + double res = elem (i, 0); + for (int j = 1; j < nc; j++) + if (elem (i, j) > res) + res = elem (i, j); + result.elem (i) = res; + } + } + + return result; +} + +ColumnVector +Matrix::row_max_loc (void) const +{ + ColumnVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nr); + + for (int i = 0; i < nr; i++) + { + int res = 0; + for (int j = 0; j < nc; j++) + if (elem (i, j) > elem (i, res)) + res = j; + result.elem (i) = (double) (res + 1); + } + } + + return result; +} + +RowVector +Matrix::column_min (void) const +{ + RowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + double res = elem (0, j); + for (int i = 1; i < nr; i++) + if (elem (i, j) < res) + res = elem (i, j); + result.elem (j) = res; + } + } + + return result; +} +RowVector +Matrix::column_min_loc (void) const +{ + RowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + int res = 0; + for (int i = 0; i < nr; i++) + if (elem (i, j) < elem (res, j)) + res = i; + result.elem (j) = (double) (res + 1); + } + } + + return result; +} + + +RowVector +Matrix::column_max (void) const +{ + RowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + double res = elem (0, j); + for (int i = 1; i < nr; i++) + if (elem (i, j) > res) + res = elem (i, j); + result.elem (j) = res; + } + } + + return result; +} + +RowVector +Matrix::column_max_loc (void) const +{ + RowVector result; + + int nr = rows (); + int nc = cols (); + + if (nr > 0 && nc > 0) + { + result.resize (nc); + + for (int j = 0; j < nc; j++) + { + int res = 0; + for (int i = 0; i < nr; i++) + if (elem (i, j) > elem (res, j)) + res = i; + result.elem (j) = (double) (res + 1); + } + } + + return result; +} + +ostream& +operator << (ostream& os, const Matrix& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.rows (); i++) + { + for (int j = 0; j < a.cols (); j++) + os << " " /* setw (field_width) */ << a.elem (i, j); + os << "\n"; + } + return os; +} + +istream& +operator >> (istream& is, Matrix& a) +{ + int nr = a.rows (); + int nc = a.cols (); + + if (nr < 1 || nc < 1) + is.clear (ios::badbit); + else + { + double tmp; + for (int i = 0; i < nr; i++) + for (int j = 0; j < nc; j++) + { + is >> tmp; + if (is) + a.elem (i, j) = tmp; + else + break; + } + } + + return is; +} + +/* + * Read an array of data froma file in binary format. + */ +int +Matrix::read (FILE *fptr, int size, Matrix::conversion conv) +{ +// Allocate buffer pointers. + + union + { + void *vd; + char *ch; + u_char *uc; +// s_char *sc; // Some systems may need this? + short *sh; + u_short *us; + int *in; + u_int *ui; + long *ln; + u_long *ul; + float *fl; + double *db; + } + buf; + + buf.db = fortran_vec (); + +// Read data directly into matrix data array. + + int count = fread (buf.ch, size, length (), fptr); + +// Convert data to double. + + int k; + + switch (conv) + { + case CNV_DOUBLE: + break; + + case CNV_CHAR: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.ch[k]; + break; + + case CNV_UCHAR: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.uc[k]; + break; + +// Some systems may need this?? +// case CNV_SCHAR: +// for (k = count - 1; k > -1; k--) +// buf.db[k] = buf.sc[k]; +// break; + + case CNV_SHORT: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.sh[k]; + break; + + case CNV_USHORT: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.us[k]; + break; + + case CNV_INT: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.in[k]; + break; + + case CNV_UINT: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.ui[k]; + break; + + case CNV_LONG: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.ln[k]; + break; + + case CNV_ULONG: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.ul[k]; + break; + + case CNV_FLOAT: + for (k = count - 1; k > -1; k--) + buf.db[k] = buf.fl[k]; + break; + + default: + return 0; + } + + return count; +} + +/* + * Write the data array to a file in binary format. + */ +int +Matrix::write (FILE *fptr, int size, Matrix::conversion conv) +{ +// Allocate buffer pointers. + + union + { + void *vd; + char *ch; + u_char *uc; +// s_char *sc; // Some systems may need this? + short *sh; + u_short *us; + int *in; + u_int *ui; + long *ln; + u_long *ul; + float *fl; + double *db; + } + buf; + + int len = length (); + + if (conv != CNV_DOUBLE) + buf.db = new double [len]; + + double *bufi = fortran_vec (); + +// Convert from double to correct size. + + int k; + + switch (conv) + { + case CNV_DOUBLE: + buf.db = bufi; + break; + + case CNV_CHAR: + for (k = 0; k < len; k++) + buf.ch[k] = (char) bufi[k]; + break; + + case CNV_UCHAR: + for (k = 0; k < len; k++) + buf.uc[k] = (u_char) bufi[k]; + break; + +// Some systems may need this? +// case CNV_SCHAR: +// for (k = 0; k < len; k++) +// buf.uc[k] = (s_char) bufi[k]; +// break; + + case CNV_SHORT: + for (k = 0; k < len; k++) + buf.sh[k] = (short) bufi[k]; + break; + + case CNV_USHORT: + for (k = 0; k < len; k++) + buf.us[k] = (u_short) bufi[k]; + break; + + case CNV_INT: + for (k = 0; k < len; k++) + buf.in[k] = (int) bufi[k]; + break; + + case CNV_UINT: + for (k = 0; k < len; k++) + buf.ui[k] = (u_int) bufi[k]; + break; + + case CNV_LONG: + for (k = 0; k < len; k++) + buf.ln[k] = (long) bufi[k]; + break; + + case CNV_ULONG: + for (k = 0; k < len; k++) + buf.ul[k] = (u_long) bufi[k]; + break; + + case CNV_FLOAT: + for (k = 0; k < len; k++) + buf.fl[k] = (float) bufi[k]; + break; + + default: + return 0; + } + +// Write data from converted matrix data array. + + int count = fwrite (buf.ch, size, length (), fptr); + + if (conv != CNV_DOUBLE) + delete [] buf.db; + + return count; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h new file mode 100644 --- /dev/null +++ b/liboctave/dMatrix.h @@ -0,0 +1,274 @@ +// -*- 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. + +*/ + +#if !defined (octave_Matrix_int_h) +#define octave_Matrix_int_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "mx-defs.h" + +extern "C++" { + +class Matrix : public Array2 +{ +friend class LU; +friend class SVD; + +public: + + Matrix (void) : Array2 () { } + Matrix (int r, int c) : Array2 (r, c) { } + Matrix (int r, int c, double val) : Array2 (r, c, val) { } + Matrix (const Array2& a) : Array2 (a) { } + Matrix (const Matrix& a) : Array2 (a) { } + Matrix (const DiagArray& a) : Array2 (a) { } + Matrix (const DiagMatrix& a); +// Matrix (double a) : Array2 (1, 1, a) { } + + Matrix& operator = (const Matrix& a) + { + Array2::operator = (a); + return *this; + } + + int operator == (const Matrix& a) const; + int operator != (const Matrix& a) const; + +// destructive insert/delete/reorder operations + + Matrix& insert (const Matrix& a, int r, int c); + Matrix& insert (const RowVector& a, int r, int c); + Matrix& insert (const ColumnVector& a, int r, int c); + Matrix& insert (const DiagMatrix& a, int r, int c); + + Matrix& fill (double val); + Matrix& fill (double val, int r1, int c1, int r2, int c2); + + Matrix append (const Matrix& a) const; + Matrix append (const RowVector& a) const; + Matrix append (const ColumnVector& a) const; + Matrix append (const DiagMatrix& a) const; + + Matrix stack (const Matrix& a) const; + Matrix stack (const RowVector& a) const; + Matrix stack (const ColumnVector& a) const; + Matrix stack (const DiagMatrix& a) const; + + Matrix transpose (void) const; + +// resize is the destructive equivalent for this one + + Matrix extract (int r1, int c1, int r2, int c2) const; + +// extract row or column i. + + RowVector row (int i) const; + RowVector row (char *s) const; + + ColumnVector column (int i) const; + ColumnVector column (char *s) const; + + Matrix inverse (void) const; + Matrix inverse (int& info) const; + Matrix inverse (int& info, double& rcond) const; + + ComplexMatrix fourier (void) const; + ComplexMatrix ifourier (void) const; + + DET determinant (void) const; + DET determinant (int& info) const; + DET determinant (int& info, double& rcond) const; + + Matrix solve (const Matrix& b) const; + Matrix solve (const Matrix& b, int& info) const; + Matrix solve (const Matrix& b, int& info, double& rcond) const; + + ComplexMatrix solve (const ComplexMatrix& b) const; + ComplexMatrix solve (const ComplexMatrix& b, int& info) const; + ComplexMatrix solve (const ComplexMatrix& b, int& info, double& rcond) const; + + ColumnVector solve (const ColumnVector& b) const; + ColumnVector solve (const ColumnVector& b, int& info) const; + ColumnVector solve (const ColumnVector& b, int& info, double& rcond) const; + + ComplexColumnVector solve (const ComplexColumnVector& b) const; + ComplexColumnVector solve (const ComplexColumnVector& b, int& info) const; + ComplexColumnVector solve (const ComplexColumnVector& b, int& info, + double& rcond) const; + + Matrix lssolve (const Matrix& b) const; + Matrix lssolve (const Matrix& b, int& info) const; + Matrix lssolve (const Matrix& b, int& info, int& rank) const; + + ComplexMatrix lssolve (const ComplexMatrix& b) const; + ComplexMatrix lssolve (const ComplexMatrix& b, int& info) const; + ComplexMatrix lssolve (const ComplexMatrix& b, int& info, + int& rank) const; + + ColumnVector lssolve (const ColumnVector& b) const; + ColumnVector lssolve (const ColumnVector& b, int& info) const; + ColumnVector lssolve (const ColumnVector& b, int& info, int& rank) const; + + ComplexColumnVector lssolve (const ComplexColumnVector& b) const; + ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info) const; + ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info, + int& rank) const; + + Matrix& operator += (const Matrix& a); + Matrix& operator -= (const Matrix& a); + + Matrix& operator += (const DiagMatrix& a); + Matrix& operator -= (const DiagMatrix& a); + +// unary operations + + Matrix operator ! (void) const; + +// matrix by scalar -> matrix operations + + friend ComplexMatrix operator + (const Matrix& a, const Complex& s); + friend ComplexMatrix operator - (const Matrix& a, const Complex& s); + friend ComplexMatrix operator * (const Matrix& a, const Complex& s); + friend ComplexMatrix operator / (const Matrix& a, const Complex& s); + +// scalar by matrix -> matrix operations + + friend ComplexMatrix operator + (const Complex& s, const Matrix& a); + friend ComplexMatrix operator - (const Complex& s, const Matrix& a); + friend ComplexMatrix operator * (const Complex& s, const Matrix& a); + friend ComplexMatrix operator / (const Complex& s, const Matrix& a); + +// matrix by column vector -> column vector operations + + friend ColumnVector operator * (const Matrix& a, const ColumnVector& b); + friend ComplexColumnVector operator * (const Matrix& a, + const ComplexColumnVector& b); + +// matrix by diagonal matrix -> matrix operations + + friend Matrix operator + (const Matrix& a, const DiagMatrix& b); + friend Matrix operator - (const Matrix& a, const DiagMatrix& b); + friend Matrix operator * (const Matrix& a, const DiagMatrix& b); + + friend ComplexMatrix operator + (const Matrix& a, + const ComplexDiagMatrix& b); + friend ComplexMatrix operator - (const Matrix& a, + const ComplexDiagMatrix& b); + friend ComplexMatrix operator * (const Matrix& a, + const ComplexDiagMatrix& b); + +// matrix by matrix -> matrix operations + + friend Matrix operator * (const Matrix& a, const Matrix& b); + friend ComplexMatrix operator * (const Matrix& a, const ComplexMatrix& b); + + friend ComplexMatrix operator + (const Matrix& a, const ComplexMatrix& b); + friend ComplexMatrix operator - (const Matrix& a, const ComplexMatrix& b); + + friend ComplexMatrix product (const Matrix& a, const ComplexMatrix& b); + friend ComplexMatrix quotient (const Matrix& a, const ComplexMatrix& b); + +// other operations + + friend Matrix map (d_d_Mapper f, const Matrix& a); + void map (d_d_Mapper f); + + Matrix all (void) const; + Matrix any (void) const; + + Matrix cumprod (void) const; + Matrix cumsum (void) const; + Matrix prod (void) const; + Matrix sum (void) const; + Matrix sumsq (void) const; + + ColumnVector diag (void) const; + ColumnVector diag (int k) const; + + ColumnVector row_min (void) const; + ColumnVector row_min_loc (void) const; + + ColumnVector row_max (void) const; + ColumnVector row_max_loc (void) const; + + RowVector column_min (void) const; + RowVector column_min_loc (void) const; + + RowVector column_max (void) const; + RowVector column_max_loc (void) const; + +// i/o + + friend ostream& operator << (ostream& os, const Matrix& a); + friend istream& operator >> (istream& is, Matrix& a); + + enum conversion + { + CNV_UNKNOWN, + CNV_UCHAR, + CNV_CHAR, + CNV_SCHAR, + CNV_SHORT, + CNV_USHORT, + CNV_INT, + CNV_UINT, + CNV_LONG, + CNV_ULONG, + CNV_FLOAT, + CNV_DOUBLE, + }; + + + int read (FILE *fptr, int size, Matrix::conversion); + int write (FILE *fptr, int size, Matrix::conversion); + +// Until templates really work with g++: + +#define KLUDGE_MATRICES +#define TYPE double +#define KL_MAT_TYPE Matrix +#include "mx-kludge.h" +#undef KLUDGE_MATRICES +#undef TYPE +#undef KL_MAT_TYPE + +private: + + Matrix (double *d, int r, int c) : Array2 (d, r, c) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dRowVector.cc b/liboctave/dRowVector.cc new file mode 100644 --- /dev/null +++ b/liboctave/dRowVector.cc @@ -0,0 +1,516 @@ +// RowVector manipulations. -*- 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 + +#include + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "lo-error.h" +#include "f77-uscore.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (dgemv) (const char*, const int*, const int*, + const double*, const double*, const int*, + const double*, const int*, const double*, + double*, const int*, long); + + double F77_FCN (ddot) (const int*, const double*, const int*, + const double*, const int*); +} + +/* + * Row Vector class. + */ + +#define KLUDGE_VECTORS +#define TYPE double +#define KL_VEC_TYPE RowVector +#include "mx-kludge.cc" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +#if 0 +RowVector& +RowVector::resize (int n) +{ + if (n < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimension"); + return *this; + } + + double *new_data = (double *) NULL; + if (n > 0) + { + new_data = new double [n]; + int min_len = len < n ? len : n; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + len = n; + data = new_data; + + return *this; +} + +RowVector& +RowVector::resize (int n, double val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} +#endif + +int +RowVector::operator == (const RowVector& a) const +{ + int len = length (); + if (len != a.length ()) + return 0; + return equal (data (), a.data (), len); +} + +int +RowVector::operator != (const RowVector& a) const +{ + return !(*this == a); +} + +RowVector& +RowVector::insert (const RowVector& a, int c) +{ + int a_len = a.length (); + if (c < 0 || c + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (c+i) = a.elem (i); + + return *this; +} + +RowVector& +RowVector::fill (double val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +RowVector& +RowVector::fill (double val, int c1, int c2) +{ + int len = length (); + if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + for (int i = c1; i <= c2; i++) + elem (i) = val; + + return *this; +} + +RowVector +RowVector::append (const RowVector& a) const +{ + int len = length (); + int nc_insert = len; + RowVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nc_insert); + return retval; +} + +ColumnVector +RowVector::transpose (void) const +{ + int len = length (); + return ColumnVector (dup (data (), len), len); +} + +RowVector +RowVector::extract (int c1, int c2) const +{ + if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } + + int new_c = c2 - c1 + 1; + + RowVector result (new_c); + + for (int i = 0; i < new_c; i++) + result.elem (i) = elem (c1+i); + + return result; +} + +// row vector by row vector -> row vector operations + +RowVector& +RowVector::operator += (const RowVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +RowVector& +RowVector::operator -= (const RowVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +// row vector by scalar -> row vector operations + +ComplexRowVector +operator + (const RowVector& v, const Complex& s) +{ + int len = v.length (); + return ComplexRowVector (add (v.data (), len, s), len); +} + +ComplexRowVector +operator - (const RowVector& v, const Complex& s) +{ + int len = v.length (); + return ComplexRowVector (subtract (v.data (), len, s), len); +} + +ComplexRowVector +operator * (const RowVector& v, const Complex& s) +{ + int len = v.length (); + return ComplexRowVector (multiply (v.data (), len, s), len); +} + +ComplexRowVector +operator / (const RowVector& v, const Complex& s) +{ + int len = v.length (); + return ComplexRowVector (divide (v.data (), len, s), len); +} + +// scalar by row vector -> row vector operations + +ComplexRowVector +operator + (const Complex& s, const RowVector& a) +{ + return ComplexRowVector (); +} + +ComplexRowVector +operator - (const Complex& s, const RowVector& a) +{ + return ComplexRowVector (); +} + +ComplexRowVector +operator * (const Complex& s, const RowVector& a) +{ + return ComplexRowVector (); +} + +ComplexRowVector +operator / (const Complex& s, const RowVector& a) +{ + return ComplexRowVector (); +} + +// row vector by column vector -> scalar + +double +operator * (const RowVector& v, const ColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return 0.0; + } + + int i_one = 1; + return F77_FCN (ddot) (&len, v.data (), &i_one, a.data (), &i_one); +} + +Complex +operator * (const RowVector& v, const ComplexColumnVector& a) +{ + ComplexRowVector tmp (v); + return tmp * a; +} + +// row vector by matrix -> row vector + +RowVector +operator * (const RowVector& v, const Matrix& a) +{ + int len = v.length (); + if (a.rows () != len) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return RowVector (); + } + + if (len == 0 || a.cols () == 0) + return RowVector (0); + +// Transpose A to form A'*x == (x'*A)' + + int a_nr = a.rows (); + int a_nc = a.cols (); + + char trans = 'T'; + int ld = a_nr; + double alpha = 1.0; + double beta = 0.0; + int i_one = 1; + + double *y = new double [len]; + + F77_FCN (dgemv) (&trans, &a_nc, &a_nr, &alpha, a.data (), &ld, + v.data (), &i_one, &beta, y, &i_one, 1L); + + return RowVector (y, len); +} + +ComplexRowVector +operator * (const RowVector& v, const ComplexMatrix& a) +{ + ComplexRowVector tmp (v); + return tmp * a; +} + +// row vector by row vector -> row vector operations + +ComplexRowVector +operator + (const RowVector& v, const ComplexRowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector addition attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (add (v.data (), a.data (), len), len); +} + +ComplexRowVector +operator - (const RowVector& v, const ComplexRowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (subtract (v.data (), a.data (), len), len); +} + +ComplexRowVector +product (const RowVector& v, const ComplexRowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector product attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (multiply (v.data (), a.data (), len), len); +} + +ComplexRowVector +quotient (const RowVector& v, const ComplexRowVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector quotient attempted"); + return ComplexRowVector (); + } + + if (len == 0) + return ComplexRowVector (0); + + return ComplexRowVector (divide (v.data (), a.data (), len), len); +} + +// other operations + +RowVector +map (d_d_Mapper f, const RowVector& a) +{ + RowVector b (a); + b.map (f); + return b; +} + +void +RowVector::map (d_d_Mapper f) +{ + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); +} + +double +RowVector::min (void) const +{ + int len = length (); + if (len == 0) + return 0; + + double res = elem (0); + + for (int i = 1; i < len; i++) + if (elem (i) < res) + res = elem (i); + + return res; +} + +double +RowVector::max (void) const +{ + int len = length (); + if (len == 0) + return 0; + + double res = elem (0); + + for (int i = 1; i < len; i++) + if (elem (i) > res) + res = elem (i); + + return res; +} + +ostream& +operator << (ostream& os, const RowVector& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.length (); i++) + os << " " /* setw (field_width) */ << a.elem (i); + return os; +} + +istream& +operator >> (istream& is, RowVector& a) +{ + int len = a.length(); + + if (len < 1) + is.clear (ios::badbit); + else + { + double tmp; + for (int i = 0; i < len; i++) + { + is >> tmp; + if (is) + a.elem (i) = tmp; + else + break; + } + } +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ diff --git a/liboctave/dRowVector.h b/liboctave/dRowVector.h new file mode 100644 --- /dev/null +++ b/liboctave/dRowVector.h @@ -0,0 +1,155 @@ +// -*- 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. + +*/ + +#if !defined (octave_RowVector_h) +#define octave_RowVector_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +#include "Array.h" + +#include "mx-defs.h" + +extern "C++" { + +class RowVector : public Array +{ +friend class ColumnVector; + +public: + + RowVector (void) : Array () { } + RowVector (int n) : Array (n) { } + RowVector (int n, double val) : Array (n, val) { } + RowVector (const Array& a) : Array (a) { } + RowVector (const RowVector& a) : Array (a) { } +// RowVector (double a) : Array (1, a) { } + + RowVector& operator = (const RowVector& a) + { + Array::operator = (a); + return *this; + } + +// operator Array& () const { return *this; } + + int operator == (const RowVector& a) const; + int operator != (const RowVector& a) const; + +// destructive insert/delete/reorder operations + + RowVector& insert (const RowVector& a, int c); + + RowVector& fill (double val); + RowVector& fill (double val, int c1, int c2); + + RowVector append (const RowVector& a) const; + + ColumnVector transpose (void) const; + +// resize is the destructive equivalent for this one + + RowVector extract (int c1, int c2) const; + +// row vector by row vector -> row vector operations + + RowVector& operator += (const RowVector& a); + RowVector& operator -= (const RowVector& a); + +// row vector by scalar -> row vector operations + + friend ComplexRowVector operator + (const RowVector& a, const Complex& s); + friend ComplexRowVector operator - (const RowVector& a, const Complex& s); + friend ComplexRowVector operator * (const RowVector& a, const Complex& s); + friend ComplexRowVector operator / (const RowVector& a, const Complex& s); + +// scalar by row vector -> row vector operations + + friend ComplexRowVector operator + (const Complex& s, const RowVector& a); + friend ComplexRowVector operator - (const Complex& s, const RowVector& a); + friend ComplexRowVector operator * (const Complex& s, const RowVector& a); + friend ComplexRowVector operator / (const Complex& s, const RowVector& a); + +// row vector by column vector -> scalar + + friend double operator * (const RowVector& a, const ColumnVector& b); + + friend Complex operator * (const RowVector& a, const ComplexColumnVector& b); + +// row vector by matrix -> row vector + + friend RowVector operator * (const RowVector& a, const Matrix& b); + + friend ComplexRowVector operator * (const RowVector& a, + const ComplexMatrix& b); + +// row vector by row vector -> row vector operations + + friend ComplexRowVector operator + (const RowVector& a, + const ComplexRowVector& b); + friend ComplexRowVector operator - (const RowVector& a, + const ComplexRowVector& b); + + friend ComplexRowVector product (const RowVector& a, + const ComplexRowVector& b); + friend ComplexRowVector quotient (const RowVector& a, + const ComplexRowVector& b); + +// other operations + + friend RowVector map (d_d_Mapper f, const RowVector& a); + void map (d_d_Mapper f); + + double min (void) const; + double max (void) const; + +// i/o + + friend ostream& operator << (ostream& os, const RowVector& a); + friend ostream& operator >> (ostream& is, RowVector& a); + +#define KLUDGE_VECTORS +#define TYPE double +#define KL_VEC_TYPE RowVector +#include "mx-kludge.h" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +private: + + RowVector (double *d, int l) : Array (d, l) { } +}; + +} // extern "C++" + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/