# HG changeset patch # User jwe # Date 845158298 0 # Node ID 5b04dfde8d55fb9420d35c7b6080ff6b455ecb7b # Parent 39a6500cde4fb54380763151f06f03a6fcbe260e [project @ 1996-10-12 22:11:37 by jwe] diff --git a/liboctave/ColVector.cc b/liboctave/ColVector.cc deleted file mode 100644 --- a/liboctave/ColVector.cc +++ /dev/null @@ -1,1092 +0,0 @@ -// 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 - -#include - -#include "Matrix.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); - -/* - * f2c translates complex*16 as - * - * typedef struct { doublereal re, im; } doublecomplex; - * - * and Complex.h from libg++ uses - * - * protected: - * double re; - * double im; - * - * as the only data members, so this should work (fingers crossed that - * things don't change). - */ - - 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); -} - -/* - * Column Vector class. - */ - -#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; - } - } -} - -/* - * Complex Column Vector class - */ - -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/DiagMatrix.cc b/liboctave/DiagMatrix.cc deleted file mode 100644 --- a/liboctave/DiagMatrix.cc +++ /dev/null @@ -1,2044 +0,0 @@ -// 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 - -#include - -#include "Matrix.h" -#include "mx-inlines.cc" -#include "lo-error.h" - -/* - * Diagonal Matrix class. - */ - -#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; -} - -/* - * Complex Diagonal Matrix class - */ - -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/Matrix.cc b/liboctave/Matrix.cc deleted file mode 100644 --- a/liboctave/Matrix.cc +++ /dev/null @@ -1,5032 +0,0 @@ -// 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 - -#include -#include - -#include "Matrix.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*); - -/* - * f2c translates complex*16 as - * - * typedef struct { doublereal re, im; } doublecomplex; - * - * and Complex.h from libg++ uses - * - * protected: - * double re; - * double im; - * - * as the only data members, so this should work (fingers crossed that - * things don't change). - */ - - 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*); -} - -// Since this is only temporary, put all of this here, rather than -// putting each type where it logically belongs. This way, it will be -// easier to delete. - -#define KLUDGE_MATRICES -#define TYPE double -#define KL_MAT_TYPE Matrix -#include "mx-kludge.cc" -#undef KLUDGE_MATRICES -#undef TYPE -#undef KL_MAT_TYPE - -#define KLUDGE_VECTORS -#define TYPE double -#define KL_VEC_TYPE ColumnVector -#include "mx-kludge.cc" -#undef KLUDGE_VECTORS -#undef TYPE -#undef KL_VEC_TYPE - -#define KLUDGE_VECTORS -#define TYPE double -#define KL_VEC_TYPE RowVector -#include "mx-kludge.cc" -#undef KLUDGE_VECTORS -#undef TYPE -#undef KL_VEC_TYPE - -#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 - -#define KLUDGE_MATRICES -#define TYPE Complex -#define KL_MAT_TYPE ComplexMatrix -#include "mx-kludge.cc" -#undef KLUDGE_MATRICES -#undef TYPE -#undef KL_MAT_TYPE - -#define KLUDGE_VECTORS -#define TYPE Complex -#define KL_VEC_TYPE ComplexColumnVector -#include "mx-kludge.cc" -#undef KLUDGE_VECTORS -#undef TYPE -#undef KL_VEC_TYPE - -#define KLUDGE_VECTORS -#define TYPE Complex -#define KL_VEC_TYPE ComplexRowVector -#include "mx-kludge.cc" -#undef KLUDGE_VECTORS -#undef TYPE -#undef KL_VEC_TYPE - -#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 - -/* - * 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; - return DET (d); - } - - 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; - } - 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; -} - -/* - * 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; - return ComplexDET (d); - } - - 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; - } - 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; -} - -/* - * 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/RowVector.cc b/liboctave/RowVector.cc deleted file mode 100644 --- a/liboctave/RowVector.cc +++ /dev/null @@ -1,1144 +0,0 @@ -// 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 - -#include "Matrix.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*); - -/* - * f2c translates complex*16 as - * - * typedef struct { doublereal re, im; } doublecomplex; - * - * and Complex.h from libg++ uses - * - * protected: - * double re; - * double im; - * - * as the only data members, so this should work (fingers crossed that - * things don't change). - */ - - 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); -} - -/* - * Row Vector class. - */ - -#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; - } - } -} - -/* - * Complex Row Vector class - */ - -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/mx-kludge.cc b/liboctave/mx-kludge.cc deleted file mode 100644 --- a/liboctave/mx-kludge.cc +++ /dev/null @@ -1,301 +0,0 @@ -// kludge.cc -*- C++ -*- -/* - -Copyright (C) 1996 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Nothing like a little CPP abuse to brighten everyone's day. Would -// have been nice to do this with template functions but as of 2.5.x, -// g++ seems to fail in various ways, either not resolving general -// template functions, or not instantiating non-member template -// functions. -// -// When templates work more reliably in g++, this will be replaced by -// the MArray class. - -#define DO_VS_OP(OP) \ - int l = a.length (); \ - TYPE *result = 0; \ - if (l > 0) \ - { \ - result = new TYPE [l]; \ - const TYPE *x = a.data (); \ - for (int i = 0; i < l; i++) \ - result[i] = x[i] OP s; \ - } - -#define DO_SV_OP(OP) \ - int l = a.length (); \ - TYPE *result = 0; \ - if (l > 0) \ - { \ - result = new TYPE [l]; \ - const TYPE *x = a.data (); \ - for (int i = 0; i < l; i++) \ - result[i] = s OP x[i]; \ - } - -#define DO_VV_OP(OP) \ - TYPE *result = 0; \ - if (l > 0) \ - { \ - result = new TYPE [l]; \ - const TYPE *x = a.data (); \ - const TYPE *y = b.data (); \ - for (int i = 0; i < l; i++) \ - result[i] = x[i] OP y[i]; \ - } - -#define NEG_V \ - int l = a.length (); \ - TYPE *result = 0; \ - if (l > 0) \ - { \ - result = new TYPE [l]; \ - const TYPE *x = a.data (); \ - for (int i = 0; i < l; i++) \ - result[i] = -x[i]; \ - } - -#ifdef KLUDGE_VECTORS - -// Like type operations for vectors. - -// Element by element vector by scalar ops. - -#define KL_VS_OP(OP) \ - KL_VEC_TYPE \ - operator OP (const KL_VEC_TYPE& a, const TYPE& s) \ - { \ - DO_VS_OP (OP); \ - return KL_VEC_TYPE (result, l); \ - } - -KL_VS_OP (+) -KL_VS_OP (-) -KL_VS_OP (*) -KL_VS_OP (/) - -// Element by element scalar by vector ops. - -#define KL_SV_OP(OP) \ - KL_VEC_TYPE \ - operator OP (const TYPE& s, const KL_VEC_TYPE& a) \ - { \ - DO_SV_OP (OP); \ - return KL_VEC_TYPE (result, l); \ - } - -KL_SV_OP (+) -KL_SV_OP (-) -KL_SV_OP (*) -KL_SV_OP (/) - -// Element by element vector by vector ops. - -#define KL_VV_OP(FCN, OP, OP_STR) \ - KL_VEC_TYPE \ - FCN (const KL_VEC_TYPE& a, const KL_VEC_TYPE& b) \ - { \ - int l = a.length (); \ - if (l != b.length ()) \ - { \ - (*current_liboctave_error_handler) \ - ("nonconformant array " OP_STR " attempted"); \ - return KL_VEC_TYPE (); \ - } \ - if (l == 0) \ - return KL_VEC_TYPE (); \ - DO_VV_OP (OP); \ - return KL_VEC_TYPE (result, l); \ - } - -KL_VV_OP(operator +, +, "addition") -KL_VV_OP(operator -, -, "subtraction") -KL_VV_OP(product, *, "product") -KL_VV_OP(quotient, /, "quotient") - -// Unary MArray ops. - -KL_VEC_TYPE -operator - (const KL_VEC_TYPE& a) -{ - NEG_V; - return KL_VEC_TYPE (result, l); -} - -#endif - -#ifdef KLUDGE_MATRICES - -// Like type operations for matrices. - -// Element by element matrix by scalar ops. - -#define KL_MS_OP(OP) \ - KL_MAT_TYPE \ - operator OP (const KL_MAT_TYPE& a, const TYPE& s) \ - { \ - DO_VS_OP (OP); \ - return KL_MAT_TYPE (result, a.rows (), a.cols ()); \ - } - -KL_MS_OP(+) -KL_MS_OP(-) -KL_MS_OP(*) -KL_MS_OP(/) - -// Element by element scalar by matrix ops. - -#define KL_SM_OP(OP) \ - KL_MAT_TYPE \ - operator OP (const TYPE& s, const KL_MAT_TYPE& a) \ - { \ - DO_SV_OP (OP); \ - return KL_MAT_TYPE (result, a.rows (), a.cols ()); \ - } - -KL_SM_OP(+) -KL_SM_OP(-) -KL_SM_OP(*) -KL_SM_OP(/) - -// Element by element matrix by matrix ops. - -#define KL_MM_OP(FCN, OP, OP_STR) \ - KL_MAT_TYPE \ - FCN (const KL_MAT_TYPE& a, const KL_MAT_TYPE& b) \ - { \ - int r = a.rows (); \ - int c = a.cols (); \ - if (r != b.rows () || c != b.cols ()) \ - { \ - (*current_liboctave_error_handler) \ - ("nonconformant array " OP_STR " attempted"); \ - return KL_MAT_TYPE (); \ - } \ - if (r == 0 || c == 0) \ - return KL_MAT_TYPE (r, c); \ - int l = a.length (); \ - DO_VV_OP (+); \ - return KL_MAT_TYPE (result, r, c); \ - } - -KL_MM_OP (operator +, +, "addition") -KL_MM_OP (operator -, -, "subtraction") -KL_MM_OP (product, *, "product") -KL_MM_OP (quotient, /, "quotient") - -// Unary matrix ops. - -KL_MAT_TYPE -operator - (const KL_MAT_TYPE& a) -{ - NEG_V; - return KL_MAT_TYPE (result, a.rows (), a.cols ()); -} - -#endif - -#ifdef KLUDGE_DIAG_MATRICES - -// Like type operations for diagonal matrices. - -// Element by element MDiagArray by scalar ops. - -#define KL_DMS_OP(OP) \ - KL_DMAT_TYPE \ - operator OP (const KL_DMAT_TYPE& a, const TYPE& s) \ - { \ - DO_VS_OP (OP); \ - return KL_DMAT_TYPE (result, a.rows (), a.cols ()); \ - } - -KL_DMS_OP (*) -KL_DMS_OP (/) - -// Element by element scalar by MDiagArray ops. - -#define KL_SDM_OP(OP) \ - KL_DMAT_TYPE \ - operator OP (const TYPE& s, const KL_DMAT_TYPE& a) \ - { \ - DO_SV_OP (OP); \ - return KL_DMAT_TYPE (result, a.rows (), a.cols ()); \ - } - -KL_SDM_OP (*) - -// Element by element MDiagArray by MDiagArray ops. - -#define KL_DMDM_OP(FCN, OP, OP_STR) \ - KL_DMAT_TYPE \ - FCN (const KL_DMAT_TYPE& a, const KL_DMAT_TYPE& b) \ - { \ - int r = a.rows (); \ - int c = a.cols (); \ - if (r != b.rows () || c != b.cols ()) \ - { \ - (*current_liboctave_error_handler) \ - ("nonconformant diagonal array " OP_STR " attempted"); \ - return KL_DMAT_TYPE (); \ - } \ - if (c == 0 || r == 0) \ - return KL_DMAT_TYPE (); \ - int l = a.length (); \ - DO_VV_OP (OP); \ - return KL_DMAT_TYPE (result, r, c); \ - } - -KL_DMDM_OP (operator +, +, "addition") -KL_DMDM_OP (operator -, -, "subtraction") -KL_DMDM_OP (product, *, "product") - -// Unary MDiagArray ops. - -KL_DMAT_TYPE -operator - (const KL_DMAT_TYPE& a) -{ - NEG_V; - return KL_DMAT_TYPE (result, a.rows (), a.cols ()); -} - -#endif - -#undef DO_VS_OP -#undef DO_SV_OP -#undef DO_VV_OP -#undef NEG_V -#undef KL_VS_OP -#undef KL_SV_OP -#undef KL_VV_OP -#undef KL_MS_OP -#undef KL_SM_OP -#undef KL_MM_OP -#undef KL_DMS_OP -#undef KL_SDM_OP -#undef KL_DMDM_OP - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; page-delimiter: "^/\\*" *** -;;; End: *** -*/ diff --git a/liboctave/mx-kludge.h b/liboctave/mx-kludge.h deleted file mode 100644 --- a/liboctave/mx-kludge.h +++ /dev/null @@ -1,128 +0,0 @@ -// kludge.h -*- C++ -*- -/* - -Copyright (C) 1996 John W. Eaton - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 2, or (at your option) any -later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, write to the Free -Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - -*/ - -// Nothing like a little CPP abuse to brighten everyone's day. Would -// have been nice to do this with template functions but as of 2.5.x, -// g++ seems to fail in various ways, either not resolving general -// template functions, or not instatiating non-member template -// functions. -// -// When templates work more reliably in g++, this will be replaced by -// the MArray class. - -#ifdef KLUDGE_VECTORS - -// Like type operations for vectors. - -// Element by element vector by scalar ops. - -friend KL_VEC_TYPE operator + (const KL_VEC_TYPE& a, const TYPE& s); -friend KL_VEC_TYPE operator - (const KL_VEC_TYPE& a, const TYPE& s); -friend KL_VEC_TYPE operator * (const KL_VEC_TYPE& a, const TYPE& s); -friend KL_VEC_TYPE operator / (const KL_VEC_TYPE& a, const TYPE& s); - -// Element by element scalar by vector ops. - -friend KL_VEC_TYPE operator + (const TYPE& s, const KL_VEC_TYPE& a); -friend KL_VEC_TYPE operator - (const TYPE& s, const KL_VEC_TYPE& a); -friend KL_VEC_TYPE operator * (const TYPE& s, const KL_VEC_TYPE& a); -friend KL_VEC_TYPE operator / (const TYPE& s, const KL_VEC_TYPE& a); - -// Element by element vector by vector ops. - -friend KL_VEC_TYPE operator + (const KL_VEC_TYPE& a, const KL_VEC_TYPE& b); -friend KL_VEC_TYPE operator - (const KL_VEC_TYPE& a, const KL_VEC_TYPE& b); - -friend KL_VEC_TYPE product (const KL_VEC_TYPE& a, const KL_VEC_TYPE& b); -friend KL_VEC_TYPE quotient (const KL_VEC_TYPE& a, const KL_VEC_TYPE& b); - -// Unary MArray ops. - -friend KL_VEC_TYPE operator - (const KL_VEC_TYPE& a); - -#endif - -#ifdef KLUDGE_MATRICES - -// Like type operations for matrices - -// Element by element matrix by scalar ops. - -friend KL_MAT_TYPE operator + (const KL_MAT_TYPE& a, const TYPE& s); -friend KL_MAT_TYPE operator - (const KL_MAT_TYPE& a, const TYPE& s); -friend KL_MAT_TYPE operator * (const KL_MAT_TYPE& a, const TYPE& s); -friend KL_MAT_TYPE operator / (const KL_MAT_TYPE& a, const TYPE& s); - -// Element by element scalar by matrix ops. - -friend KL_MAT_TYPE operator + (const TYPE& s, const KL_MAT_TYPE& a); -friend KL_MAT_TYPE operator - (const TYPE& s, const KL_MAT_TYPE& a); -friend KL_MAT_TYPE operator * (const TYPE& s, const KL_MAT_TYPE& a); -friend KL_MAT_TYPE operator / (const TYPE& s, const KL_MAT_TYPE& a); - -// Element by element matrix by matrix ops. - -friend KL_MAT_TYPE operator + (const KL_MAT_TYPE& a, const KL_MAT_TYPE& b); -friend KL_MAT_TYPE operator - (const KL_MAT_TYPE& a, const KL_MAT_TYPE& b); - -friend KL_MAT_TYPE product (const KL_MAT_TYPE& a, const KL_MAT_TYPE& b); -friend KL_MAT_TYPE quotient (const KL_MAT_TYPE& a, const KL_MAT_TYPE& b); - -// Unary matrix ops. - -friend KL_MAT_TYPE operator - (const KL_MAT_TYPE& a); - -#endif - -#ifdef KLUDGE_DIAG_MATRICES - -// Like type operations for diagonal matrices. - -// Element by element MDiagArray by scalar ops. - -friend KL_DMAT_TYPE operator * (const KL_DMAT_TYPE& a, const TYPE& s); -friend KL_DMAT_TYPE operator / (const KL_DMAT_TYPE& a, const TYPE& s); - -// Element by element scalar by MDiagArray ops. - -friend KL_DMAT_TYPE operator * (const TYPE& s, const KL_DMAT_TYPE& a); - -// Element by element MDiagArray by MDiagArray ops. - -friend KL_DMAT_TYPE operator + (const KL_DMAT_TYPE& a, const KL_DMAT_TYPE& b); -friend KL_DMAT_TYPE operator - (const KL_DMAT_TYPE& a, const KL_DMAT_TYPE& b); - -friend KL_DMAT_TYPE product (const KL_DMAT_TYPE& a, const KL_DMAT_TYPE& b); - -// Unary MDiagArray ops. - -friend KL_DMAT_TYPE operator - (const KL_DMAT_TYPE& a); - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; page-delimiter: "^/\\*" *** -;;; End: *** -*/