Mercurial > octave
changeset 1205:8302fab9fe24
[project @ 1995-04-04 02:05:01 by jwe]
| author | jwe |
|---|---|
| date | Tue, 04 Apr 1995 02:05:01 +0000 |
| parents | 68d147abe7ca |
| children | 83b59118ab1a |
| files | liboctave/CColVector.cc liboctave/CColVector.h liboctave/CDiagMatrix.cc liboctave/CDiagMatrix.h liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/CRowVector.cc liboctave/CRowVector.h liboctave/dColVector.cc liboctave/dColVector.h liboctave/dDiagMatrix.cc liboctave/dDiagMatrix.h liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/dRowVector.cc liboctave/dRowVector.h |
| diffstat | 16 files changed, 1653 insertions(+), 2241 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CColVector.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CColVector.cc Tue Apr 04 02:05:01 1995 +0000 @@ -38,11 +38,10 @@ extern "C" { - int F77_FCN (zgemm) (const char*, const char*, const int*, - const int*, const int*, const Complex*, + int F77_FCN (zgemv) (const char*, const int*, const int*, + const Complex*, const Complex*, const int*, const Complex*, const int*, const Complex*, - const int*, const Complex*, Complex*, const int*, - long, long); + Complex*, const int*, long); } /* @@ -64,57 +63,6 @@ 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 = 0; - 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 { @@ -256,26 +204,6 @@ 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) { @@ -416,6 +344,34 @@ return ComplexColumnVector (divide (v.data (), len, s), len); } +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 @@ -446,36 +402,70 @@ return ComplexColumnVector (divide (s, a.data (), a_len), a_len); } -// column vector by row vector -> matrix 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); +} -ComplexMatrix -operator * (const ComplexColumnVector& v, const ComplexRowVector& a) +ComplexColumnVector +operator / (const Complex& s, const ColumnVector& a) { - int len = v.length (); int a_len = a.length (); - if (len != a_len) + return ComplexColumnVector (divide (s, a.data (), a_len), a_len); +} + +// 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 vector multiplication attempted"); - return ComplexMatrix (); + ("nonconformant matrix multiplication attempted"); + return ComplexColumnVector (); } - if (len == 0) - return ComplexMatrix (len, len, 0.0); + if (nc == 0 || nr == 0) + return ComplexColumnVector (0); - char transa = 'N'; - char transb = 'N'; + char trans = 'N'; + int ld = nr; Complex alpha (1.0); Complex beta (0.0); - int anr = 1; + int i_one = 1; - Complex *c = new Complex [len * a_len]; + Complex *y = new Complex [nr]; - F77_FCN (zgemm) (&transa, &transb, &len, &a_len, &anr, &alpha, - v.data (), &len, a.data (), &anr, &beta, c, &len, - 1L, 1L); + F77_FCN (zgemv) (&trans, &nr, &nc, &alpha, m.data (), &ld, a.data (), + &i_one, &beta, y, &i_one, 1L); - return ComplexMatrix (c, len, a_len); + return ComplexColumnVector (y, nr); } // column vector by column vector -> column vector operations @@ -515,6 +505,40 @@ } 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 ComplexColumnVector& v, const ColumnVector& a) { int len = v.length (); @@ -548,6 +572,132 @@ return ComplexColumnVector (divide (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); +} + +// matrix by column vector -> column vector operations + +ComplexColumnVector +operator * (const Matrix& m, const ComplexColumnVector& a) +{ + ComplexMatrix tmp (m); + return tmp * a; +} + +// diagonal matrix by column vector -> column vector operations + +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; +} + +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; +} + // other operations ComplexColumnVector @@ -558,16 +708,6 @@ 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) {
--- a/liboctave/CColVector.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CColVector.h Tue Apr 04 02:05:01 1995 +0000 @@ -32,21 +32,17 @@ class ComplexColumnVector : public Array<Complex> { -friend class ColumnVector; +friend class ComplexMatrix; friend class ComplexRowVector; -friend class ComplexMatrix; public: ComplexColumnVector (void) : Array<Complex> () { } ComplexColumnVector (int n) : Array<Complex> (n) { } - ComplexColumnVector (int n, const Complex& val) - : Array<Complex> (n, val) { } + ComplexColumnVector (int n, const Complex& val) : Array<Complex> (n, val) { } ComplexColumnVector (const ColumnVector& a); ComplexColumnVector (const Array<Complex>& a) : Array<Complex> (a) { } ComplexColumnVector (const ComplexColumnVector& a) : Array<Complex> (a) { } -// ComplexColumnVector (double a) : Array<Complex> (1, a) { } -// ComplexColumnVector (const Complex& a) : Array<Complex> (1, a) { } ComplexColumnVector& operator = (const ComplexColumnVector& a) { @@ -75,8 +71,6 @@ ComplexRowVector hermitian (void) const; // complex conjugate transpose. ComplexRowVector transpose (void) const; - friend ColumnVector real (const ComplexColumnVector& a); - friend ColumnVector imag (const ComplexColumnVector& a); friend ComplexColumnVector conj (const ComplexColumnVector& a); // resize is the destructive equivalent for this one @@ -102,6 +96,15 @@ friend ComplexColumnVector operator / (const ComplexColumnVector& a, double s); + friend ComplexColumnVector operator + (const ColumnVector& a, + const Complex& s); + friend ComplexColumnVector operator - (const ColumnVector& a, + const Complex& s); + friend ComplexColumnVector operator * (const ColumnVector& a, + const Complex& s); + friend ComplexColumnVector operator / (const ColumnVector& a, + const Complex& s); + // scalar by column vector -> column vector operations friend ComplexColumnVector operator + (double s, @@ -113,10 +116,22 @@ friend ComplexColumnVector operator / (double s, const ComplexColumnVector& a); -// column vector by row vector -> matrix operations + friend ComplexColumnVector operator + (const Complex& s, + const ColumnVector& a); + friend ComplexColumnVector operator - (const Complex& s, + const ColumnVector& a); + friend ComplexColumnVector operator * (const Complex& s, + const ColumnVector& a); + friend ComplexColumnVector operator / (const Complex& s, + const ColumnVector& a); - friend ComplexMatrix operator * (const ComplexColumnVector& a, - const ComplexRowVector& b); +// matrix by column vector -> column vector operations + + friend ComplexColumnVector operator * (const ComplexMatrix& a, + const ColumnVector& b); + + friend ComplexColumnVector operator * (const ComplexMatrix& a, + const ComplexColumnVector& b); // column vector by column vector -> column vector operations @@ -125,15 +140,40 @@ friend ComplexColumnVector operator - (const ComplexColumnVector& a, const ColumnVector& b); + friend ComplexColumnVector operator + (const ColumnVector& a, + const ComplexColumnVector& b); + friend ComplexColumnVector operator - (const ColumnVector& a, + const ComplexColumnVector& b); + friend ComplexColumnVector product (const ComplexColumnVector& a, const ColumnVector& b); friend ComplexColumnVector quotient (const ComplexColumnVector& a, const ColumnVector& b); + friend ComplexColumnVector product (const ColumnVector& a, + const ComplexColumnVector& b); + friend ComplexColumnVector quotient (const ColumnVector& a, + const ComplexColumnVector& b); + +// matrix by column vector -> column vector operations + + friend ComplexColumnVector operator * (const Matrix& a, + const ComplexColumnVector& b); + +// diagonal matrix by column vector -> column vector operations + + friend ComplexColumnVector operator * (const DiagMatrix& a, + const ComplexColumnVector& b); + + friend ComplexColumnVector operator * (const ComplexDiagMatrix& a, + const ColumnVector& b); + + friend ComplexColumnVector operator * (const ComplexDiagMatrix& a, + const ComplexColumnVector& b); + // other operations friend ComplexColumnVector map (c_c_Mapper f, const ComplexColumnVector& a); - friend ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a); void map (c_c_Mapper f); Complex min (void) const; @@ -152,10 +192,8 @@ #undef TYPE #undef KL_VEC_TYPE -// private: -// XXX FIXME XXX -- why does it not work to make this private, with -// ColumnVector declared as a friend of ComplexColumnVector? It seems -// to work for the similar case with Matrix/ComplexMatrix. Hmm... +private: + ComplexColumnVector (Complex *d, int l) : Array<Complex> (d, l) { } };
--- a/liboctave/CDiagMatrix.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CDiagMatrix.cc Tue Apr 04 02:05:01 1995 +0000 @@ -66,107 +66,6 @@ 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 = 0; - 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 = 0; - 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 = 0; - 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 { @@ -182,12 +81,6 @@ return !(*this == a); } -ComplexDiagMatrix -ComplexDiagMatrix::hermitian (void) const -{ - return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ()); -} - ComplexDiagMatrix& ComplexDiagMatrix::fill (double val) { @@ -363,33 +256,17 @@ } ComplexDiagMatrix +ComplexDiagMatrix::hermitian (void) const +{ + return ComplexDiagMatrix (conj_dup (data (), length ()), cols (), rows ()); +} + +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) { @@ -622,36 +499,6 @@ 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 @@ -668,34 +515,18 @@ a.rows (), a.cols ()); } -// scalar by diagonal matrix -> matrix operations - -ComplexMatrix -operator + (double s, const ComplexDiagMatrix& a) +ComplexDiagMatrix +operator * (const DiagMatrix& a, const Complex& s) { - ComplexMatrix tmp (a.rows (), a.cols (), s); - return tmp + a; + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); } -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) +ComplexDiagMatrix +operator / (const DiagMatrix& a, const Complex& s) { - 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; + return ComplexDiagMatrix (divide (a.data (), a.length (), s), + a.rows (), a.cols ()); } // scalar by diagonal matrix -> diagonal matrix operations @@ -707,60 +538,11 @@ a.rows (), a.cols ()); } -// diagonal matrix by column vector -> column vector operations - -ComplexColumnVector -operator * (const ComplexDiagMatrix& m, const ColumnVector& a) +ComplexDiagMatrix +operator * (const Complex& s, const DiagMatrix& 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; + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); } // diagonal matrix by diagonal matrix -> diagonal matrix operations @@ -881,6 +663,82 @@ } 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 ComplexDiagMatrix& m, const DiagMatrix& a) { int nr = m.rows (); @@ -899,190 +757,23 @@ 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) +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 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 (); + ("nonconformant matrix product attempted"); + return ComplexDiagMatrix (); } - 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; -} + if (nc == 0 || nr == 0) + return ComplexDiagMatrix (nr, nc); -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; + return ComplexDiagMatrix (multiply (m.data (), a.data (), m.length ()), + nr, nc); } // other operations
--- a/liboctave/CDiagMatrix.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CDiagMatrix.h Tue Apr 04 02:05:01 1995 +0000 @@ -37,8 +37,6 @@ class ComplexDiagMatrix : public DiagArray<Complex> { - friend DiagMatrix; - public: ComplexDiagMatrix (void) : DiagArray<Complex> () { } @@ -57,7 +55,6 @@ ComplexDiagMatrix (const DiagArray<Complex>& a) : DiagArray<Complex> (a) { } ComplexDiagMatrix (const ComplexDiagMatrix& a) : DiagArray<Complex> (a) { } -// ComplexDiagMatrix (const Complex& a) : DiagArray<Complex> (1, a) { } ComplexDiagMatrix& operator = (const ComplexDiagMatrix& a) { @@ -65,8 +62,6 @@ return *this; } -// operator DiagArray<Complex>& () const { return *this; } - int operator == (const ComplexDiagMatrix& a) const; int operator != (const ComplexDiagMatrix& a) const; @@ -86,8 +81,6 @@ ComplexDiagMatrix hermitian (void) const; // complex conjugate transpose ComplexDiagMatrix transpose (void) const; - friend DiagMatrix real (const ComplexDiagMatrix& a); - friend DiagMatrix imag (const ComplexDiagMatrix& a); friend ComplexDiagMatrix conj (const ComplexDiagMatrix& a); // resize is the destructive analog for this one @@ -113,42 +106,19 @@ ComplexDiagMatrix& operator += (const ComplexDiagMatrix& a); ComplexDiagMatrix& operator -= (const ComplexDiagMatrix& a); -// diagonal matrix by scalar -> matrix operations - - friend ComplexMatrix operator + (const ComplexDiagMatrix& a, double s); - friend ComplexMatrix operator - (const ComplexDiagMatrix& a, double s); - - friend ComplexMatrix operator + (const ComplexDiagMatrix& a, - const Complex& s); - friend ComplexMatrix operator - (const ComplexDiagMatrix& a, - const Complex& s); - // diagonal matrix by scalar -> diagonal matrix operations friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, double s); friend ComplexDiagMatrix operator / (const ComplexDiagMatrix& a, double s); -// scalar by diagonal matrix -> matrix operations - - friend ComplexMatrix operator + (double s, const ComplexDiagMatrix& a); - friend ComplexMatrix operator - (double s, const ComplexDiagMatrix& a); - - friend ComplexMatrix operator + (const Complex& s, - const ComplexDiagMatrix& a); - friend ComplexMatrix operator - (const Complex& s, - const ComplexDiagMatrix& a); + friend ComplexDiagMatrix operator * (const DiagMatrix& a, const Complex& s); + friend ComplexDiagMatrix operator / (const DiagMatrix& a, const Complex& s); // scalar by diagonal matrix -> diagonal matrix operations friend ComplexDiagMatrix operator * (double s, const ComplexDiagMatrix& a); -// diagonal matrix by column vector -> column vector operations - - friend ComplexColumnVector operator * (const ComplexDiagMatrix& a, - const ColumnVector& b); - - friend ComplexColumnVector operator * (const ComplexDiagMatrix& a, - const ComplexColumnVector& b); + friend ComplexDiagMatrix operator * (const Complex& s, const DiagMatrix& a); // diagonal matrix by diagonal matrix -> diagonal matrix operations @@ -162,24 +132,18 @@ friend ComplexDiagMatrix operator * (const ComplexDiagMatrix& a, const DiagMatrix& b); + friend ComplexDiagMatrix operator + (const DiagMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexDiagMatrix operator - (const DiagMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexDiagMatrix operator * (const DiagMatrix& a, + const ComplexDiagMatrix& b); + friend ComplexDiagMatrix product (const ComplexDiagMatrix& a, const DiagMatrix& b); -// diagonal matrix by matrix -> matrix operations - - friend ComplexMatrix operator + (const ComplexDiagMatrix& a, - const Matrix& b); - friend ComplexMatrix operator - (const ComplexDiagMatrix& a, - const Matrix& b); - friend ComplexMatrix operator * (const ComplexDiagMatrix& a, - const Matrix& b); - - friend ComplexMatrix operator + (const ComplexDiagMatrix& a, - const ComplexMatrix& b); - friend ComplexMatrix operator - (const ComplexDiagMatrix& a, - const ComplexMatrix& b); - friend ComplexMatrix operator * (const ComplexDiagMatrix& a, - const ComplexMatrix& b); + friend ComplexDiagMatrix product (const DiagMatrix& a, + const ComplexDiagMatrix& b); // other operations
--- a/liboctave/CMatrix.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CMatrix.cc Tue Apr 04 02:05:01 1995 +0000 @@ -48,11 +48,6 @@ 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*); @@ -113,117 +108,6 @@ 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 = 0; - 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 = 0; - 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 = 0; - 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 { @@ -770,26 +654,6 @@ 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) { @@ -1518,6 +1382,140 @@ return retval; } +// column vector by row vector -> matrix operations + +ComplexMatrix +operator * (const ColumnVector& v, const ComplexRowVector& a) +{ + ComplexColumnVector tmp (v); + return tmp * a; +} + +ComplexMatrix +operator * (const ComplexColumnVector& a, const RowVector& b) +{ + ComplexRowVector tmp (b); + return a * tmp; +} + +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); +} + +// diagonal matrix by scalar -> matrix operations + +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; +} + +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; +} + +// scalar by diagonal matrix -> matrix operations + +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; +} + +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; +} + // matrix by diagonal matrix -> matrix operations ComplexMatrix& @@ -1592,6 +1590,378 @@ return *this; } +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 = 0; + + 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); +} + +// diagonal matrix by matrix -> matrix operations + +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; +} + +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; +} + // matrix by matrix -> matrix operations ComplexMatrix& @@ -1689,6 +2059,34 @@ // 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 ()); +} + +ComplexMatrix operator + (const ComplexMatrix& a, double s) { return ComplexMatrix (add (a.data (), a.length (), s), @@ -1746,42 +2144,32 @@ a.rows (), a.cols ()); } -// matrix by column vector -> column vector operations - -ComplexColumnVector -operator * (const ComplexMatrix& m, const ColumnVector& a) +ComplexMatrix +operator + (const Complex& s, const Matrix& a) { - ComplexColumnVector tmp (a); - return m * tmp; + return ComplexMatrix (add (s, a.data (), a.length ()), + a.rows (), a.cols ()); } -ComplexColumnVector -operator * (const ComplexMatrix& m, const ComplexColumnVector& a) +ComplexMatrix +operator - (const Complex& s, const Matrix& 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); + 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 diagonal matrix -> matrix operations @@ -2011,6 +2399,39 @@ } 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 operator * (const ComplexMatrix& m, const Matrix& a) { ComplexMatrix tmp (a); @@ -2018,6 +2439,13 @@ } ComplexMatrix +operator * (const Matrix& m, const ComplexMatrix& a) +{ + ComplexMatrix tmp (m); + return tmp * a; +} + +ComplexMatrix operator * (const ComplexMatrix& m, const ComplexMatrix& a) { int nr = m.rows (); @@ -2086,6 +2514,42 @@ return ComplexMatrix (divide (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 ComplexMatrix @@ -2096,18 +2560,6 @@ 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) {
--- a/liboctave/CMatrix.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CMatrix.h Tue Apr 04 02:05:01 1995 +0000 @@ -34,6 +34,7 @@ class ComplexMatrix : public Array2<Complex> { +friend class Matrix; friend class ComplexCHOL; friend class ComplexHESS; friend class ComplexLU; @@ -41,8 +42,6 @@ friend class ComplexQRP; friend class ComplexSCHUR; friend class ComplexSVD; -friend class ComplexColumnVector; -friend class Matrix; public: @@ -56,8 +55,6 @@ ComplexMatrix (const DiagMatrix& a); ComplexMatrix (const DiagArray<Complex>& a) : Array2<Complex> (a) { } ComplexMatrix (const ComplexDiagMatrix& a); -// ComplexMatrix (double a) : Array2<Complex> (1, 1, a) { } -// ComplexMatrix (const Complex& a) : Array2<Complex> (1, 1, a) { } ComplexMatrix& operator = (const ComplexMatrix& a) { @@ -65,8 +62,6 @@ return *this; } -// operator Array2<Complex>& () const { return *this; } - int operator == (const ComplexMatrix& a) const; int operator != (const ComplexMatrix& a) const; @@ -110,8 +105,6 @@ ComplexMatrix hermitian (void) const; // complex conjugate transpose ComplexMatrix transpose (void) const; - friend Matrix real (const ComplexMatrix& a); - friend Matrix imag (const ComplexMatrix& a); friend ComplexMatrix conj (const ComplexMatrix& a); // resize is the destructive equivalent for this one @@ -165,6 +158,43 @@ ComplexColumnVector lssolve (const ComplexColumnVector& b, int& info, int& rank) const; +// column vector by row vector -> matrix operations + + friend ComplexMatrix operator * (const ColumnVector& a, + const ComplexRowVector& b); + + friend ComplexMatrix operator * (const ComplexColumnVector& a, + const RowVector& b); + + friend ComplexMatrix operator * (const ComplexColumnVector& a, + const ComplexRowVector& b); + +// diagonal matrix by scalar -> matrix operations + + friend ComplexMatrix operator + (const DiagMatrix& a, const Complex& s); + friend ComplexMatrix operator - (const DiagMatrix& a, const Complex& s); + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, double s); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, double s); + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, + const Complex& s); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, + const Complex& s); + +// scalar by diagonal matrix -> matrix operations + + friend ComplexMatrix operator + (const Complex& s, const DiagMatrix& a); + friend ComplexMatrix operator - (const Complex& s, const DiagMatrix& a); + + friend ComplexMatrix operator + (double s, const ComplexDiagMatrix& a); + friend ComplexMatrix operator - (double s, const ComplexDiagMatrix& a); + + friend ComplexMatrix operator + (const Complex& s, + const ComplexDiagMatrix& a); + friend ComplexMatrix operator - (const Complex& s, + const ComplexDiagMatrix& a); + // matrix by diagonal matrix -> matrix operations ComplexMatrix& operator += (const DiagMatrix& a); @@ -173,6 +203,36 @@ ComplexMatrix& operator += (const ComplexDiagMatrix& a); ComplexMatrix& operator -= (const ComplexDiagMatrix& a); + friend ComplexMatrix operator + (const Matrix& a, + const ComplexDiagMatrix& b); + friend ComplexMatrix operator - (const Matrix& a, + const ComplexDiagMatrix& b); + friend ComplexMatrix operator * (const Matrix& a, + const ComplexDiagMatrix& b); + +// diagonal matrix by matrix -> matrix operations + + friend ComplexMatrix operator + (const DiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator - (const DiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator * (const DiagMatrix& a, + const ComplexMatrix& b); + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, + const Matrix& b); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, + const Matrix& b); + friend ComplexMatrix operator * (const ComplexDiagMatrix& a, + const Matrix& b); + + friend ComplexMatrix operator + (const ComplexDiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator - (const ComplexDiagMatrix& a, + const ComplexMatrix& b); + friend ComplexMatrix operator * (const ComplexDiagMatrix& a, + const ComplexMatrix& b); + // matrix by matrix -> matrix operations ComplexMatrix& operator += (const Matrix& a); @@ -187,6 +247,11 @@ // matrix by scalar -> matrix operations + friend ComplexMatrix operator + (const Matrix& a, const Complex& s); + friend ComplexMatrix operator - (const Matrix& a, const Complex& s); + friend ComplexMatrix operator * (const Matrix& a, const Complex& s); + friend ComplexMatrix operator / (const Matrix& a, const Complex& s); + friend ComplexMatrix operator + (const ComplexMatrix& a, double s); friend ComplexMatrix operator - (const ComplexMatrix& a, double s); friend ComplexMatrix operator * (const ComplexMatrix& a, double s); @@ -199,13 +264,10 @@ friend ComplexMatrix operator * (double s, const ComplexMatrix& a); friend ComplexMatrix operator / (double s, const ComplexMatrix& a); -// matrix by column vector -> column vector operations - - friend ComplexColumnVector operator * (const ComplexMatrix& a, - const ColumnVector& b); - - friend ComplexColumnVector operator * (const ComplexMatrix& a, - const ComplexColumnVector& b); + friend ComplexMatrix operator + (const Complex& s, const Matrix& a); + friend ComplexMatrix operator - (const Complex& s, const Matrix& a); + friend ComplexMatrix operator * (const Complex& s, const Matrix& a); + friend ComplexMatrix operator / (const Complex& s, const Matrix& a); // matrix by diagonal matrix -> matrix operations @@ -228,17 +290,25 @@ friend ComplexMatrix operator + (const ComplexMatrix& a, const Matrix& b); friend ComplexMatrix operator - (const ComplexMatrix& a, const Matrix& b); + friend ComplexMatrix operator + (const Matrix& a, const ComplexMatrix& b); + friend ComplexMatrix operator - (const Matrix& a, const ComplexMatrix& b); + friend ComplexMatrix operator * (const ComplexMatrix& a, const Matrix& b); + + friend ComplexMatrix operator * (const Matrix& a, const ComplexMatrix& b); + friend ComplexMatrix operator * (const ComplexMatrix& a, const ComplexMatrix& b); friend ComplexMatrix product (const ComplexMatrix& a, const Matrix& b); friend ComplexMatrix quotient (const ComplexMatrix& a, const Matrix& b); + friend ComplexMatrix product (const Matrix& a, const ComplexMatrix& b); + friend ComplexMatrix quotient (const Matrix& a, const ComplexMatrix& b); + // other operations friend ComplexMatrix map (c_c_Mapper f, const ComplexMatrix& a); - friend Matrix map (d_c_Mapper f, const ComplexMatrix& a); void map (c_c_Mapper f); Matrix all (void) const;
--- a/liboctave/CRowVector.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CRowVector.cc Tue Apr 04 02:05:01 1995 +0000 @@ -63,57 +63,6 @@ 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 = 0; - 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 { @@ -255,26 +204,6 @@ 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) { @@ -414,6 +343,34 @@ return ComplexRowVector (divide (v.data (), len, s), len); } +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 @@ -444,32 +401,28 @@ return ComplexRowVector (divide (s, a.data (), a_len), a_len); } -// row vector by column vector -> scalar - -Complex -operator * (const ComplexRowVector& v, const ColumnVector& a) +ComplexRowVector +operator + (const Complex& s, const RowVector& a) { - ComplexColumnVector tmp (a); - return v * tmp; + return ComplexRowVector (); } -Complex -operator * (const ComplexRowVector& v, const ComplexColumnVector& a) +ComplexRowVector +operator - (const Complex& s, const RowVector& a) { - int len = v.length (); - if (len != a.length ()) - { - (*current_liboctave_error_handler) - ("nonconformant vector multiplication attempted"); - return 0.0; - } + return ComplexRowVector (); +} - Complex retval (0.0, 0.0); +ComplexRowVector +operator * (const Complex& s, const RowVector& a) +{ + return ComplexRowVector (); +} - for (int i = 0; i < len; i++) - retval += v.elem (i) * a.elem (i); - - return retval; +ComplexRowVector +operator / (const Complex& s, const RowVector& a) +{ + return ComplexRowVector (); } // row vector by matrix -> row vector @@ -507,6 +460,13 @@ return ComplexRowVector (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 @@ -544,6 +504,40 @@ } 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 ComplexRowVector& v, const RowVector& a) { int len = v.length (); @@ -577,6 +571,40 @@ return ComplexRowVector (divide (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 ComplexRowVector @@ -587,16 +615,6 @@ 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) { @@ -604,24 +622,6 @@ elem (i) = f (elem (i)); } -ComplexRowVector -linspace (const Complex& x1, const Complex& x2, int n) -{ - ComplexRowVector retval; - - if (n > 0) - { - retval.resize (n); - Complex delta = (x2 - x1) / (n - 1); - retval.elem (0) = x1; - for (int i = 1; i < n-1; i++) - retval.elem (i) = x1 + i * delta; - retval.elem (n-1) = x2; - } - - return retval; -} - Complex ComplexRowVector::min (void) const { @@ -695,6 +695,56 @@ return is; } +// row vector by column vector -> scalar + +// 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; +} + +// other operations + +ComplexRowVector +linspace (const Complex& x1, const Complex& x2, int n) +{ + ComplexRowVector retval; + + if (n > 0) + { + retval.resize (n); + Complex delta = (x2 - x1) / (n - 1); + retval.elem (0) = x1; + for (int i = 1; i < n-1; i++) + retval.elem (i) = x1 + i * delta; + retval.elem (n-1) = x2; + } + + return retval; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/CRowVector.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/CRowVector.h Tue Apr 04 02:05:01 1995 +0000 @@ -33,7 +33,6 @@ class ComplexRowVector : public Array<Complex> { friend class ComplexColumnVector; -friend class RowVector; public: @@ -43,8 +42,6 @@ ComplexRowVector (const RowVector& a); ComplexRowVector (const Array<Complex>& a) : Array<Complex> (a) { } ComplexRowVector (const ComplexRowVector& a) : Array<Complex> (a) { } -// ComplexRowVector (double a) : Array<Complex> (1, a) { } -// ComplexRowVector (const Complex& a) : Array<Complex> (1, a) { } ComplexRowVector& operator = (const ComplexRowVector& a) { @@ -52,8 +49,6 @@ return *this; } -// operator Array<Complex>& () const { return *this; } - int operator == (const ComplexRowVector& a) const; int operator != (const ComplexRowVector& a) const; @@ -73,8 +68,6 @@ ComplexColumnVector hermitian (void) const; // complex conjugate transpose. ComplexColumnVector transpose (void) const; - friend RowVector real (const ComplexRowVector& a); - friend RowVector imag (const ComplexRowVector& a); friend ComplexRowVector conj (const ComplexRowVector& a); // resize is the destructive equivalent for this one @@ -96,6 +89,11 @@ friend ComplexRowVector operator * (const ComplexRowVector& a, double s); friend ComplexRowVector operator / (const ComplexRowVector& a, double s); + friend ComplexRowVector operator + (const RowVector& a, const Complex& s); + friend ComplexRowVector operator - (const RowVector& a, const Complex& s); + friend ComplexRowVector operator * (const RowVector& a, const Complex& s); + friend ComplexRowVector operator / (const RowVector& a, const Complex& s); + // scalar by row vector -> row vector operations friend ComplexRowVector operator + (double s, const ComplexRowVector& a); @@ -103,18 +101,19 @@ friend ComplexRowVector operator * (double s, const ComplexRowVector& a); friend ComplexRowVector operator / (double s, const ComplexRowVector& a); -// row vector by column vector -> scalar - - friend Complex operator * (const ComplexRowVector& a, const ColumnVector& b); - - friend Complex operator * (const ComplexRowVector& a, - const ComplexColumnVector& b); + friend ComplexRowVector operator + (const Complex& s, const RowVector& a); + friend ComplexRowVector operator - (const Complex& s, const RowVector& a); + friend ComplexRowVector operator * (const Complex& s, const RowVector& a); + friend ComplexRowVector operator / (const Complex& s, const RowVector& a); // row vector by matrix -> row vector friend ComplexRowVector operator * (const ComplexRowVector& a, const ComplexMatrix& b); + friend ComplexRowVector operator * (const RowVector& a, + const ComplexMatrix& b); + // row vector by row vector -> row vector operations friend ComplexRowVector operator + (const ComplexRowVector& a, @@ -122,15 +121,24 @@ friend ComplexRowVector operator - (const ComplexRowVector& a, const RowVector& b); + friend ComplexRowVector operator + (const RowVector& a, + const ComplexRowVector& b); + friend ComplexRowVector operator - (const RowVector& a, + const ComplexRowVector& b); + friend ComplexRowVector product (const ComplexRowVector& a, const RowVector& b); friend ComplexRowVector quotient (const ComplexRowVector& a, const RowVector& b); + friend ComplexRowVector product (const RowVector& a, + const ComplexRowVector& b); + friend ComplexRowVector quotient (const RowVector& a, + const ComplexRowVector& b); + // other operations friend ComplexRowVector map (c_c_Mapper f, const ComplexRowVector& a); - friend RowVector map (d_c_Mapper f, const ComplexRowVector& a); void map (c_c_Mapper f); Complex min (void) const; @@ -154,6 +162,14 @@ ComplexRowVector (Complex *d, int l) : Array<Complex> (d, l) { } }; +// row vector by column vector -> scalar + +Complex operator * (const ComplexRowVector& a, const ColumnVector& b); + +Complex operator * (const ComplexRowVector& a, const ComplexColumnVector& b); + +// other operations + ComplexRowVector linspace (const Complex& x1, const Complex& x2, int n); } // extern "C++"
--- a/liboctave/dColVector.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dColVector.cc Tue Apr 04 02:05:01 1995 +0000 @@ -38,11 +38,10 @@ extern "C" { - int F77_FCN (dgemm) (const char*, const char*, const int*, - const int*, const int*, const double*, + int F77_FCN (dgemv) (const char*, const int*, const int*, + const double*, const double*, const int*, const double*, const int*, const double*, - const int*, const double*, double*, const int*, - long, long); + double*, const int*, long); } /* @@ -57,46 +56,6 @@ #undef TYPE #undef KL_VEC_TYPE -#if 0 -ColumnVector& -ColumnVector::resize (int n) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return *this; - } - - double *new_data = 0; - 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 { @@ -174,6 +133,26 @@ return RowVector (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; +} + // resize is the destructive equivalent for this one ColumnVector @@ -233,171 +212,64 @@ 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 +// matrix 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) +ColumnVector +operator * (const Matrix& m, 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) + int nr = m.rows (); + int nc = m.cols (); + if (nc != a.length ()) { (*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"); + ("nonconformant matrix multiplication attempted"); return ColumnVector (); } - if (len == 0) - return ComplexColumnVector (0); + if (nr == 0 || nc == 0) + return ColumnVector (0); - return ComplexColumnVector (multiply (v.data (), a.data (), len), len); + 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 -quotient (const ColumnVector& v, const ComplexColumnVector& a) +// diagonal matrix by column vector -> column vector operations + +ColumnVector +operator * (const DiagMatrix& m, const ColumnVector& a) { - int len = v.length (); - if (len != a.length ()) + int nr = m.rows (); + int nc = m.cols (); + int a_len = a.length (); + if (nc != a_len) { (*current_liboctave_error_handler) - ("nonconformant vector quotient attempted"); + ("nonconformant matrix multiplication attempted"); return ColumnVector (); } - if (len == 0) - return ComplexColumnVector (0); + if (nc == 0 || nr == 0) + return ColumnVector (0); + + ColumnVector result (nr); - return ComplexColumnVector (divide (v.data (), a.data (), len), len); + 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; } // other operations @@ -410,6 +282,16 @@ 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 ColumnVector::map (d_d_Mapper f) {
--- a/liboctave/dColVector.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dColVector.h Tue Apr 04 02:05:01 1995 +0000 @@ -34,7 +34,6 @@ { friend class Matrix; friend class RowVector; -friend class ComplexColumnVector; public: @@ -43,7 +42,6 @@ ColumnVector (int n, double val) : Array<double> (n, val) { } ColumnVector (const Array<double>& a) : Array<double> (a) { } ColumnVector (const ColumnVector& a) : Array<double> (a) { } -// ColumnVector (double a) : Array<double> (1, a) { } ColumnVector& operator = (const ColumnVector& a) { @@ -51,8 +49,6 @@ return *this; } -// operator Array<double>& () const { return *this; } - int operator == (const ColumnVector& a) const; int operator != (const ColumnVector& a) const; @@ -67,6 +63,9 @@ RowVector transpose (void) const; + friend ColumnVector real (const ComplexColumnVector& a); + friend ColumnVector imag (const ComplexColumnVector& a); + // resize is the destructive equivalent for this one ColumnVector extract (int r1, int r2) const; @@ -76,52 +75,18 @@ ColumnVector& operator += (const ColumnVector& a); ColumnVector& operator -= (const ColumnVector& a); -// column vector by scalar -> column vector operations +// matrix by column vector -> column vector operations - friend ComplexColumnVector operator + (const ColumnVector& a, - const Complex& s); - friend ComplexColumnVector operator - (const ColumnVector& a, - const Complex& s); - friend ComplexColumnVector operator * (const ColumnVector& a, - const Complex& s); - friend ComplexColumnVector operator / (const ColumnVector& a, - const Complex& s); - -// scalar by column vector -> column vector operations + friend ColumnVector operator * (const Matrix& a, const ColumnVector& b); - friend ComplexColumnVector operator + (const Complex& s, - const ColumnVector& a); - friend ComplexColumnVector operator - (const Complex& s, - const ColumnVector& a); - friend ComplexColumnVector operator * (const Complex& s, - const ColumnVector& a); - friend ComplexColumnVector operator / (const Complex& s, - const ColumnVector& a); - -// column vector by row vector -> matrix operations - - friend Matrix operator * (const ColumnVector& a, const RowVector& a); +// diagonal matrix by column vector -> column vector operations - friend ComplexMatrix operator * (const ColumnVector& a, - const ComplexRowVector& b); - -// column vector by column vector -> column vector operations - - friend ComplexColumnVector operator + (const ComplexColumnVector& a, - const ComplexColumnVector& b); - - friend ComplexColumnVector operator - (const ComplexColumnVector& a, - const ComplexColumnVector& b); - - friend ComplexColumnVector product (const ComplexColumnVector& a, - const ComplexColumnVector& b); - - friend ComplexColumnVector quotient (const ComplexColumnVector& a, - const ComplexColumnVector& b); + friend ColumnVector operator * (const DiagMatrix& a, const ColumnVector& b); // other operations friend ColumnVector map (d_d_Mapper f, const ColumnVector& a); + friend ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a); void map (d_d_Mapper f); double min (void) const;
--- a/liboctave/dDiagMatrix.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dDiagMatrix.cc Tue Apr 04 02:05:01 1995 +0000 @@ -45,73 +45,6 @@ #undef TYPE #undef KL_DMAT_TYPE -#if 0 -DiagMatrix& -DiagMatrix::resize (int r, int c) -{ - if (r < 0 || c < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimensions"); - return *this; - } - - int new_len = r < c ? r : c; - double *new_data = 0; - 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 = 0; - 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 { @@ -220,6 +153,28 @@ return DiagMatrix (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; +} + Matrix DiagMatrix::extract (int r1, int c1, int r2, int c2) const { @@ -400,147 +355,6 @@ 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 @@ -582,287 +396,6 @@ 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
--- a/liboctave/dDiagMatrix.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dDiagMatrix.h Tue Apr 04 02:05:01 1995 +0000 @@ -37,7 +37,6 @@ { friend class SVD; friend class ComplexSVD; -friend class ComplexDiagMatrix; public: @@ -72,6 +71,9 @@ DiagMatrix transpose (void) const; + friend DiagMatrix real (const ComplexDiagMatrix& a); + friend DiagMatrix imag (const ComplexDiagMatrix& a); + // resize is the destructive analog for this one Matrix extract (int r1, int c1, int r2, int c2) const; @@ -92,66 +94,11 @@ DiagMatrix& operator += (const DiagMatrix& a); DiagMatrix& operator -= (const DiagMatrix& a); -// diagonal matrix by scalar -> matrix operations - - friend Matrix operator + (const DiagMatrix& a, double s); - friend Matrix operator - (const DiagMatrix& a, double s); - - friend ComplexMatrix operator + (const DiagMatrix& a, const Complex& s); - friend ComplexMatrix operator - (const DiagMatrix& a, const Complex& s); - -// diagonal matrix by scalar -> diagonal matrix operations - - friend ComplexDiagMatrix operator * (const DiagMatrix& a, const Complex& s); - friend ComplexDiagMatrix operator / (const DiagMatrix& a, const Complex& s); - -// scalar by diagonal matrix -> matrix operations - - friend Matrix operator + (double s, const DiagMatrix& a); - friend Matrix operator - (double s, const DiagMatrix& a); - - friend ComplexMatrix operator + (const Complex& s, const DiagMatrix& a); - friend ComplexMatrix operator - (const Complex& s, const DiagMatrix& a); - -// scalar by diagonal matrix -> diagonal matrix operations - - friend ComplexDiagMatrix operator * (const Complex& s, const DiagMatrix& a); - -// diagonal matrix by column vector -> column vector operations - - friend ColumnVector operator * (const DiagMatrix& a, const ColumnVector& b); - - friend ComplexColumnVector operator * (const DiagMatrix& a, - const ComplexColumnVector& b); - // diagonal matrix by diagonal matrix -> diagonal matrix operations friend DiagMatrix operator * (const DiagMatrix& a, const DiagMatrix& b); - friend ComplexDiagMatrix operator + (const DiagMatrix& a, - const ComplexDiagMatrix& b); - friend ComplexDiagMatrix operator - (const DiagMatrix& a, - const ComplexDiagMatrix& b); - friend ComplexDiagMatrix operator * (const DiagMatrix& a, - const ComplexDiagMatrix& b); - - friend ComplexDiagMatrix product (const DiagMatrix& a, - const ComplexDiagMatrix& b); - -// diagonal matrix by matrix -> matrix operations - - friend Matrix operator + (const DiagMatrix& a, const Matrix& b); - friend Matrix operator - (const DiagMatrix& a, const Matrix& b); - friend Matrix operator * (const DiagMatrix& a, const Matrix& b); - - friend ComplexMatrix operator + (const DiagMatrix& a, - const ComplexMatrix& b); - friend ComplexMatrix operator - (const DiagMatrix& a, - const ComplexMatrix& b); - friend ComplexMatrix operator * (const DiagMatrix& a, - const ComplexMatrix& b); - // other operations ColumnVector diag (void) const;
--- a/liboctave/dMatrix.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dMatrix.cc Tue Apr 04 02:05:01 1995 +0000 @@ -49,11 +49,6 @@ 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*); @@ -99,79 +94,6 @@ 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 = 0; - 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 = 0; - 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 { @@ -453,6 +375,26 @@ } 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; +} + +Matrix Matrix::extract (int r1, int c1, int r2, int c2) const { if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } @@ -1309,102 +1251,68 @@ 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 ()); -} +// column vector by row vector -> matrix operations -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) +Matrix +operator * (const ColumnVector& v, const RowVector& a) { - return ComplexMatrix (multiply (a.data (), a.length (), s), - a.rows (), a.cols ()); -} + int len = v.length (); + int a_len = a.length (); + if (len != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return Matrix (); + } -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. + if (len == 0) + return Matrix (len, len, 0.0); -ComplexMatrix -operator + (const Complex& s, const Matrix& a) -{ - return ComplexMatrix (add (s, a.data (), a.length ()), - a.rows (), a.cols ()); -} + char transa = 'N'; + char transb = 'N'; + double alpha = 1.0; + double beta = 0.0; + int anr = 1; -ComplexMatrix -operator - (const Complex& s, const Matrix& a) -{ - return ComplexMatrix (subtract (s, a.data (), a.length ()), - a.rows (), a.cols ()); + 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 Complex& s, const Matrix& a) +// diagonal matrix by scalar -> matrix operations + +Matrix +operator + (const DiagMatrix& a, double s) { - 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 tmp (a.rows (), a.cols (), s); + return a + tmp; } -// matrix by column vector -> column vector operations - -ColumnVector -operator * (const Matrix& m, const ColumnVector& a) +Matrix +operator - (const DiagMatrix& a, double s) { - 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); + Matrix tmp (a.rows (), a.cols (), -s); + return a + tmp; } -ComplexColumnVector -operator * (const Matrix& m, const ComplexColumnVector& a) +// scalar by diagonal matrix -> matrix operations + +Matrix +operator + (double s, const DiagMatrix& a) { - ComplexMatrix tmp (m); - return tmp * 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; } // matrix by diagonal matrix -> matrix operations @@ -1506,8 +1414,10 @@ return Matrix (c, nr, a_nc); } -ComplexMatrix -operator + (const Matrix& m, const ComplexDiagMatrix& a) +// diagonal matrix by matrix -> matrix operations + +Matrix +operator + (const DiagMatrix& m, const Matrix& a) { int nr = m.rows (); int nc = m.cols (); @@ -1515,21 +1425,21 @@ { (*current_liboctave_error_handler) ("nonconformant matrix addition attempted"); - return ComplexMatrix (); + return Matrix (); } if (nr == 0 || nc == 0) - return ComplexMatrix (nr, nc); + return Matrix (nr, nc); - ComplexMatrix result (m); - for (int i = 0; i < a.length (); i++) - result.elem (i, i) += a.elem (i, i); + Matrix result (a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); return result; } -ComplexMatrix -operator - (const Matrix& m, const ComplexDiagMatrix& a) +Matrix +operator - (const DiagMatrix& m, const Matrix& a) { int nr = m.rows (); int nc = m.cols (); @@ -1537,21 +1447,21 @@ { (*current_liboctave_error_handler) ("nonconformant matrix subtraction attempted"); - return ComplexMatrix (); + return Matrix (); } if (nr == 0 || nc == 0) - return ComplexMatrix (nr, nc); + return Matrix (nr, nc); - ComplexMatrix result (m); - for (int i = 0; i < a.length (); i++) - result.elem (i, i) -= a.elem (i, i); + Matrix result (-a); + for (int i = 0; i < m.length (); i++) + result.elem (i, i) += m.elem (i, i); return result; } -ComplexMatrix -operator * (const Matrix& m, const ComplexDiagMatrix& a) +Matrix +operator * (const DiagMatrix& m, const Matrix& a) { int nr = m.rows (); int nc = m.cols (); @@ -1561,43 +1471,41 @@ { (*current_liboctave_error_handler) ("nonconformant matrix multiplication attempted"); - return ComplexMatrix (); + return Matrix (); } if (nr == 0 || nc == 0 || a_nc == 0) - return ComplexMatrix (nr, a_nc, 0.0); + return Matrix (nr, a_nc, 0.0); - Complex *c = new Complex [nr*a_nc]; - Complex *ctmp = 0; + Matrix c (nr, a_nc); - for (int j = 0; j < a.length (); j++) + for (int i = 0; i < m.length (); i++) { - int idx = j * nr; - ctmp = c + idx; - if (a.elem (j, j) == 1.0) + if (m.elem (i, i) == 1.0) { - for (int i = 0; i < nr; i++) - ctmp[i] = m.elem (i, j); + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = a.elem (i, j); } - else if (a.elem (j, j) == 0.0) + else if (m.elem (i, i) == 0.0) { - for (int i = 0; i < nr; i++) - ctmp[i] = 0.0; + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = 0.0; } else { - for (int i = 0; i < nr; i++) - ctmp[i] = a.elem (j, j) * m.elem (i, j); + for (int j = 0; j < a_nc; j++) + c.elem (i, j) = m.elem (i, i) * a.elem (i, j); } } - if (a_nr < a_nc) + if (nr > nc) { - for (int i = nr * nc; i < nr * a_nc; i++) - ctmp[i] = 0.0; + for (int j = 0; j < a_nc; j++) + for (int i = a_nr; i < nr; i++) + c.elem (i, j) = 0.0; } - return ComplexMatrix (c, nr, a_nc); + return c; } // matrix by matrix -> matrix operations @@ -1636,82 +1544,6 @@ 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 @@ -1722,6 +1554,18 @@ 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 Matrix::map (d_d_Mapper f) {
--- a/liboctave/dMatrix.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dMatrix.h Tue Apr 04 02:05:01 1995 +0000 @@ -35,7 +35,7 @@ class Matrix : public Array2<double> { -friend class ColumnVector; +friend class ComplexMatrix; friend class AEPBAL; friend class CHOL; friend class GEPBAL; @@ -45,7 +45,6 @@ friend class QRP; friend class SCHUR; friend class SVD; -friend class ComplexMatrix; public: @@ -56,7 +55,6 @@ Matrix (const Matrix& a) : Array2<double> (a) { } Matrix (const DiagArray<double>& a) : Array2<double> (a) { } Matrix (const DiagMatrix& a); -// Matrix (double a) : Array2<double> (1, 1, a) { } Matrix& operator = (const Matrix& a) { @@ -89,6 +87,9 @@ Matrix transpose (void) const; + friend Matrix real (const ComplexMatrix& a); + friend Matrix imag (const ComplexMatrix& a); + // resize is the destructive equivalent for this one Matrix extract (int r1, int c1, int r2, int c2) const; @@ -162,25 +163,19 @@ Matrix operator ! (void) const; -// matrix by scalar -> matrix operations +// column vector by row vector -> matrix operations - friend ComplexMatrix operator + (const Matrix& a, const Complex& s); - friend ComplexMatrix operator - (const Matrix& a, const Complex& s); - friend ComplexMatrix operator * (const Matrix& a, const Complex& s); - friend ComplexMatrix operator / (const Matrix& a, const Complex& s); + friend Matrix operator * (const ColumnVector& a, const RowVector& a); -// scalar by matrix -> matrix operations +// diagonal matrix by scalar -> matrix operations - friend ComplexMatrix operator + (const Complex& s, const Matrix& a); - friend ComplexMatrix operator - (const Complex& s, const Matrix& a); - friend ComplexMatrix operator * (const Complex& s, const Matrix& a); - friend ComplexMatrix operator / (const Complex& s, const Matrix& a); + friend Matrix operator + (const DiagMatrix& a, double s); + friend Matrix operator - (const DiagMatrix& a, double s); -// matrix by column vector -> column vector operations +// scalar by diagonal matrix -> matrix operations - friend ColumnVector operator * (const Matrix& a, const ColumnVector& b); - friend ComplexColumnVector operator * (const Matrix& a, - const ComplexColumnVector& b); + friend Matrix operator + (double s, const DiagMatrix& a); + friend Matrix operator - (double s, const DiagMatrix& a); // matrix by diagonal matrix -> matrix operations @@ -188,27 +183,20 @@ friend Matrix operator - (const Matrix& a, const DiagMatrix& b); friend Matrix operator * (const Matrix& a, const DiagMatrix& b); - friend ComplexMatrix operator + (const Matrix& a, - const ComplexDiagMatrix& b); - friend ComplexMatrix operator - (const Matrix& a, - const ComplexDiagMatrix& b); - friend ComplexMatrix operator * (const Matrix& a, - const ComplexDiagMatrix& b); +// diagonal matrix by matrix -> matrix operations + + friend Matrix operator + (const DiagMatrix& a, const Matrix& b); + friend Matrix operator - (const DiagMatrix& a, const Matrix& b); + friend Matrix operator * (const DiagMatrix& a, const Matrix& b); // matrix by matrix -> matrix operations friend Matrix operator * (const Matrix& a, const Matrix& b); - friend ComplexMatrix operator * (const Matrix& a, const ComplexMatrix& b); - - friend ComplexMatrix operator + (const Matrix& a, const ComplexMatrix& b); - friend ComplexMatrix operator - (const Matrix& a, const ComplexMatrix& b); - - friend ComplexMatrix product (const Matrix& a, const ComplexMatrix& b); - friend ComplexMatrix quotient (const Matrix& a, const ComplexMatrix& b); // other operations friend Matrix map (d_d_Mapper f, const Matrix& a); + friend Matrix map (d_c_Mapper f, const ComplexMatrix& a); void map (d_d_Mapper f); Matrix all (void) const;
--- a/liboctave/dRowVector.cc Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dRowVector.cc Tue Apr 04 02:05:01 1995 +0000 @@ -59,46 +59,6 @@ #undef TYPE #undef KL_VEC_TYPE -#if 0 -RowVector& -RowVector::resize (int n) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't resize to negative dimension"); - return *this; - } - - double *new_data = 0; - 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 { @@ -177,6 +137,26 @@ } 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; +} + +RowVector RowVector::extract (int c1, int c2) const { if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; } @@ -233,86 +213,6 @@ 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 @@ -348,83 +248,6 @@ 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 @@ -435,6 +258,16 @@ 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 RowVector::map (d_d_Mapper f) { @@ -442,24 +275,6 @@ elem (i) = f (elem (i)); } -RowVector -linspace (double x1, double x2, int n) -{ - RowVector retval; - - if (n > 0) - { - retval.resize (n); - double delta = (x2 - x1) / (n - 1); - retval.elem (0) = x1; - for (int i = 1; i < n-1; i++) - retval.elem (i) = x1 + i * delta; - retval.elem (n-1) = x2; - } - - return retval; -} - double RowVector::min (void) const { @@ -523,6 +338,50 @@ return is; } +// other operations + +RowVector +linspace (double x1, double x2, int n) +{ + RowVector retval; + + if (n > 0) + { + retval.resize (n); + double delta = (x2 - x1) / (n - 1); + retval.elem (0) = x1; + for (int i = 1; i < n-1; i++) + retval.elem (i) = x1 + i * delta; + retval.elem (n-1) = x2; + } + + return retval; +} + +// 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; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/dRowVector.h Tue Apr 04 01:42:14 1995 +0000 +++ b/liboctave/dRowVector.h Tue Apr 04 02:05:01 1995 +0000 @@ -33,7 +33,6 @@ class RowVector : public Array<double> { friend class ColumnVector; -friend class ComplexRowVector; public: @@ -42,7 +41,6 @@ RowVector (int n, double val) : Array<double> (n, val) { } RowVector (const Array<double>& a) : Array<double> (a) { } RowVector (const RowVector& a) : Array<double> (a) { } -// RowVector (double a) : Array<double> (1, a) { } RowVector& operator = (const RowVector& a) { @@ -50,8 +48,6 @@ return *this; } -// operator Array<double>& () const { return *this; } - int operator == (const RowVector& a) const; int operator != (const RowVector& a) const; @@ -66,6 +62,9 @@ ColumnVector transpose (void) const; + friend RowVector real (const ComplexRowVector& a); + friend RowVector imag (const ComplexRowVector& a); + // resize is the destructive equivalent for this one RowVector extract (int c1, int c2) const; @@ -75,48 +74,14 @@ RowVector& operator += (const RowVector& a); RowVector& operator -= (const RowVector& a); -// row vector by scalar -> row vector operations - - friend ComplexRowVector operator + (const RowVector& a, const Complex& s); - friend ComplexRowVector operator - (const RowVector& a, const Complex& s); - friend ComplexRowVector operator * (const RowVector& a, const Complex& s); - friend ComplexRowVector operator / (const RowVector& a, const Complex& s); - -// scalar by row vector -> row vector operations - - friend ComplexRowVector operator + (const Complex& s, const RowVector& a); - friend ComplexRowVector operator - (const Complex& s, const RowVector& a); - friend ComplexRowVector operator * (const Complex& s, const RowVector& a); - friend ComplexRowVector operator / (const Complex& s, const RowVector& a); - -// row vector by column vector -> scalar - - friend double operator * (const RowVector& a, const ColumnVector& b); - - friend Complex operator * (const RowVector& a, const ComplexColumnVector& b); - // row vector by matrix -> row vector friend RowVector operator * (const RowVector& a, const Matrix& b); - friend ComplexRowVector operator * (const RowVector& a, - const ComplexMatrix& b); - -// row vector by row vector -> row vector operations - - friend ComplexRowVector operator + (const RowVector& a, - const ComplexRowVector& b); - friend ComplexRowVector operator - (const RowVector& a, - const ComplexRowVector& b); - - friend ComplexRowVector product (const RowVector& a, - const ComplexRowVector& b); - friend ComplexRowVector quotient (const RowVector& a, - const ComplexRowVector& b); - // other operations friend RowVector map (d_d_Mapper f, const RowVector& a); + friend RowVector map (d_c_Mapper f, const ComplexRowVector& a); void map (d_d_Mapper f); double min (void) const; @@ -140,6 +105,14 @@ RowVector (double *d, int l) : Array<double> (d, l) { } }; +// row vector by column vector -> scalar + +double operator * (const RowVector& a, const ColumnVector& b); + +Complex operator * (const RowVector& a, const ComplexColumnVector& b); + +// other operations + RowVector linspace (double x1, double x2, int n); } // extern "C++"