Mercurial > hg > octave-nkf
changeset 238:780cbbc57b7c
[project @ 1993-11-30 20:23:04 by jwe]
line wrap: on
line diff
--- a/libcruft/misc/lo-error.cc +++ b/libcruft/misc/lo-error.cc @@ -21,8 +21,8 @@ */ -#ifdef __GNUG__ -#pragma implementation +#ifdef HAVE_CONFIG_H +#include "config.h" #endif #include <stdio.h>
--- a/liboctave/Array.h +++ b/liboctave/Array.h @@ -21,112 +21,173 @@ */ -// Written by John C. Campbell <jcc@che.utexas.edu>. - #if !defined (_Array_h) #define _Array_h 1 -#include <iostream.h> -#include <assert.h> +#if defined (__GNUG__) && defined (USE_EXTERNAL_TEMPLATES) +#pragma interface +#endif + +// Classes we declare. -template <class T> class Array; +template <class T> class ArrayRep; +template <class T> class Array; +template <class T> class Array2; +template <class T> class Array3; +template <class T> class DiagArray; + +/* + * The real representation of all arrays. + */ template <class T> class ArrayRep { +// Rethink resize()? friend class Array<T>; + friend class Array2<T>; + friend class Array3<T>; + friend class DiagArray<T>; + +protected: + + ArrayRep (T *d, int l); public: - ArrayRep (void); - ArrayRep (int); - ArrayRep (const ArrayRep<T>& a); + ArrayRep (void); + ArrayRep (int n); + ArrayRep (const ArrayRep<T>& a); ~ArrayRep (void); - + int length (void) const; - + T& elem (int n); - T& checkelem (int n); - T& operator () (int n); - + T elem (int n) const; - T checkelem (int n) const; - T operator () (int n) const; - + + void resize (int n); + private: - + T *data; int len; int count; }; +/* + * One dimensional array class. Handles the reference counting for + * all the derived classes. + */ + template <class T> class Array { +protected: + + ArrayRep<T> *rep; + + Array (T *d, int l); + public: - + Array (void); - Array (int); - Array (int n, T val); + Array (int n); + Array (int n, const T& val); + Array (const Array<T>& a); ~Array (void); Array<T>& operator = (const Array<T>& a); - + + int capacity (void) const; int length (void) const; T& elem (int n); T& checkelem (int n); T& operator () (int n); +// No checking. + T& xelem (int n); + T elem (int n) const; T checkelem (int n) const; T operator () (int n) const; -protected: + void resize (int n); + void resize (int n, const T& val); + + const T *data (void) const; - ArrayRep<T> *rep; + T *fortran_vec (void); }; +/* + * Two dimensional array class. + */ + template <class T> class Array2 : public Array<T> { +protected: + + int d1; + int d2; + + Array2 (T *d, int n, int m); + public: Array2 (void); Array2 (int n, int m); - Array2 (int n, int m, T val); + Array2 (int n, int m, const T& val); Array2 (const Array2<T>& a); + Array2 (const DiagArray<T>& a); Array2<T>& operator = (const Array2<T>& a); int dim1 (void) const; int dim2 (void) const; - + + int rows (void) const; + int cols (void) const; + int columns (void) const; + T& elem (int i, int j); - T& checkelem (int i, int j); + T& checkelem (int i, int j); T& operator () (int i, int j); - + +// No checking. + T& xelem (int i, int j); + T elem (int i, int j) const; T checkelem (int i, int j) const; T operator () (int i, int j) const; -protected: - - int d1; - int d2; + void resize (int n, int m); + void resize (int n, int m, const T& val); }; +/* + * Three dimensional array class. + */ + template <class T> class Array3 : public Array2<T> { +protected: + + int d3; + + Array3 (T *d, int n, int m, int k); + public: Array3 (void); Array3 (int n, int m, int k); - Array3 (int n, int m, int k, T val); + Array3 (int n, int m, int k, const T& val); Array3 (const Array3<T>& a); Array3<T>& operator = (const Array3<T>& a); @@ -134,32 +195,49 @@ int dim3 (void) const; T& elem (int i, int j, int k); - T& checkelem (int i, int j, int k); - T& operator()(int i,int j,int k); - + T& checkelem (int i, int j, int k); + T& operator () (int i, int j, int k); + +// No checking. + T& xelem (int i, int j, int k); + T elem (int i, int j, int k) const; - T checkelem(int i,int j,int k)const; - T operator()(int i,int j,int k) const; + T checkelem (int i, int j, int k) const; + T operator () (int i, int j, int k) const; -protected: - - int d3; + void resize (int n, int m, int k); + void resize (int n, int m, int k, const T& val); }; +/* + * A two-dimensional array with diagonal elements only. + */ + template <class T> class DiagArray : public Array<T> { +protected: + + int nr; + int nc; + + DiagArray (T *d, int r, int c); + public: - + DiagArray (void); - DiagArray (int n): Array<T> (n) {} + DiagArray (int n); + DiagArray (int n, const T& val); DiagArray (int r, int c); - DiagArray (int r, int c, T val); + DiagArray (int r, int c, const T& val); DiagArray (const Array<T>& a); DiagArray (const DiagArray<T>& a); DiagArray<T>& operator = (const DiagArray<T>& a); + int dim1 (void) const; + int dim2 (void) const; + int rows (void) const; int cols (void) const; int columns (void) const; @@ -168,17 +246,18 @@ T& checkelem (int r, int c); T& operator () (int r, int c); +// No checking. + T& xelem (int r, int c); + T elem (int r, int c) const; T checkelem (int r, int c) const; T operator () (int r, int c) const; -protected: - - int nr; - int nc; + void resize (int n, int m); + void resize (int n, int m, const T& val); }; -#ifdef __GNUG__ +#if defined (__GNUG__) && ! defined (USE_EXTERNAL_TEMPLATES) #include "Array.cc" #endif
--- a/liboctave/Bounds.cc +++ b/liboctave/Bounds.cc @@ -21,11 +21,12 @@ */ -#ifdef __GNUG__ -#pragma implementation +#ifdef HAVE_CONFIG_H +#include "config.h" #endif #include <iostream.h> + #include "Bounds.h" #include "lo-error.h"
--- a/liboctave/Bounds.h +++ b/liboctave/Bounds.h @@ -24,11 +24,8 @@ #if !defined (_Bounds_h) #define _Bounds_h 1 -#ifdef __GNUG__ -#pragma interface -#endif +class ostream; -#include <iostream.h> #include "Matrix.h" #ifndef Vector
--- a/liboctave/ColVector.cc +++ b/liboctave/ColVector.cc @@ -21,17 +21,16 @@ */ -// I\'m not sure how this is supposed to work if the .h file declares -// several classes, each of which is defined in a separate file... -// -// #ifdef __GNUG__ -// #pragma implementation "Matrix.h" -// #endif +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include <iostream.h> #include "Matrix.h" #include "mx-inlines.cc" +#include "f77-uscore.h" #include "lo-error.h" -#include "f77-uscore.h" // Fortran functions we call. @@ -69,111 +68,7 @@ * Column Vector class. */ -ColumnVector::ColumnVector (int n) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create column vector with negative dimension"); - len = 0; - data = (double *) NULL; - return; - } - - len = n; - if (n > 0) - data = new double [len]; - else - data = (double *) NULL; -} - -ColumnVector::ColumnVector (int n, double val) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create column vector with negative dimension"); - len = 0; - data = (double *) NULL; - return; - } - - len = n; - if (n > 0) - { - data = new double [len]; - copy (data, len, val); - } - else - data = (double *) NULL; -} - -ColumnVector::ColumnVector (const ColumnVector& a) -{ - len = a.len; - if (len > 0) - { - data = new double [len]; - copy (data, a.data, len); - } - else - data = (double *) NULL; -} - -ColumnVector::ColumnVector (double a) -{ - len = 1; - data = new double [1]; - data[0] = a; -} - -ColumnVector& -ColumnVector::operator = (const ColumnVector& a) -{ - if (this != &a) - { - delete [] data; - len = a.len; - if (len > 0) - { - data = new double [len]; - copy (data, a.data, len); - } - else - data = (double *) NULL; - } - return *this; -} - -double& -ColumnVector::checkelem (int n) -{ -#ifndef NO_RANGE_CHECK - if (n < 0 || n >= len) - { - (*current_liboctave_error_handler) ("range error"); - static double foo = 0.0; - return foo; - } -#endif - - return elem (n); -} - -double -ColumnVector::checkelem (int n) const -{ -#ifndef NO_RANGE_CHECK - if (n < 0 || n >= len) - { - (*current_liboctave_error_handler) ("range error"); - return 0.0; - } -#endif - - return elem (n); -} - +#if 0 ColumnVector& ColumnVector::resize (int n) { @@ -211,34 +106,35 @@ return *this; } +#endif int ColumnVector::operator == (const ColumnVector& a) const { - if (len != a.len) + int len = length (); + if (len != a.length ()) return 0; - return equal (data, a.data, len); + return equal (data (), a.data (), len); } int ColumnVector::operator != (const ColumnVector& a) const { - if (len != a.len) - return 1; - return !equal (data, a.data, len); + return !(*this == a); } ColumnVector& ColumnVector::insert (const ColumnVector& a, int r) { - if (r < 0 || r + a.len - 1 > len) + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) { (*current_liboctave_error_handler) ("range error for insert"); return *this; } - for (int i = 0; i < a.len; i++) - data[r+i] = a.data[i]; + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); return *this; } @@ -246,14 +142,17 @@ ColumnVector& ColumnVector::fill (double val) { + int len = length (); if (len > 0) - copy (data, len, val); + for (int i = 0; i < len; i++) + elem (i) = val; return *this; } ColumnVector& ColumnVector::fill (double val, int r1, int r2) { + int len = length (); if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) { (*current_liboctave_error_handler) ("range error for fill"); @@ -263,7 +162,7 @@ if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } for (int i = r1; i <= r2; i++) - data[i] = val; + elem (i) = val; return *this; } @@ -271,8 +170,9 @@ ColumnVector ColumnVector::stack (const ColumnVector& a) const { + int len = length (); int nr_insert = len; - ColumnVector retval (len + a.len); + ColumnVector retval (len + a.length ()); retval.insert (*this, 0); retval.insert (a, nr_insert); return retval; @@ -281,7 +181,8 @@ RowVector ColumnVector::transpose (void) const { - return RowVector (dup (data, len), len); + int len = length (); + return RowVector (dup (data (), len), len); } // resize is the destructive equivalent for this one @@ -296,93 +197,121 @@ ColumnVector result (new_r); for (int i = 0; i < new_r; i++) - result.data[i] = elem (r1+i); + result.elem (i) = elem (r1+i); return result; } -// column vector by scalar -> column vector operations +// column vector by column vector -> column vector operations -ColumnVector -ColumnVector::operator + (double s) const +ColumnVector& +ColumnVector::operator += (const ColumnVector& a) { - return ColumnVector (add (data, len, s), len); + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return ColumnVector (); + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; } -ColumnVector -ColumnVector::operator - (double s) const +ColumnVector& +ColumnVector::operator -= (const ColumnVector& a) { - return ColumnVector (subtract (data, len, s), len); -} + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return ColumnVector (); + } -ColumnVector -ColumnVector::operator * (double s) const -{ - return ColumnVector (multiply (data, len, s), len); -} + if (len == 0) + return *this; -ColumnVector -ColumnVector::operator / (double s) const -{ - return ColumnVector (divide (data, len, s), len); + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; } // scalar by column vector -> column vector operations -ColumnVector -operator + (double s, const ColumnVector& a) +ComplexColumnVector +operator + (const ColumnVector& a, const Complex& s) { - return ColumnVector (add (a.data, a.len, s), a.len); + int len = a.length (); + return ComplexColumnVector (add (a.data (), len, s), len); } -ColumnVector -operator - (double s, const ColumnVector& a) +ComplexColumnVector +operator - (const ColumnVector& a, const Complex& s) { - return ColumnVector (subtract (s, a.data, a.len), a.len); + int len = a.length (); + return ComplexColumnVector (subtract (a.data (), len, s), len); } -ColumnVector -operator * (double s, const ColumnVector& a) +ComplexColumnVector +operator * (const ColumnVector& a, const Complex& s) { - return ColumnVector (multiply (a.data, a.len, s), a.len); -} - -ColumnVector -operator / (double s, const ColumnVector& a) -{ - return ColumnVector (divide (s, a.data, a.len), a.len); + int len = a.length (); + return ComplexColumnVector (multiply (a.data (), len, s), len); } ComplexColumnVector -ColumnVector::operator + (const Complex& s) const +operator / (const ColumnVector& a, const Complex& s) { - return ComplexColumnVector (add (data, len, s), len); + int len = a.length (); + return ComplexColumnVector (divide (a.data (), len, s), len); +} + +// scalar by column vector -> column vector operations + +ComplexColumnVector +operator + (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (add (a.data (), a_len, s), a_len); } ComplexColumnVector -ColumnVector::operator - (const Complex& s) const +operator - (const Complex& s, const ColumnVector& a) { - return ComplexColumnVector (subtract (data, len, s), len); + int a_len = a.length (); + return ComplexColumnVector (subtract (s, a.data (), a_len), a_len); } ComplexColumnVector -ColumnVector::operator * (const Complex& s) const +operator * (const Complex& s, const ColumnVector& a) { - return ComplexColumnVector (multiply (data, len, s), len); + int a_len = a.length (); + return ComplexColumnVector (multiply (a.data (), a_len, s), a_len); } ComplexColumnVector -ColumnVector::operator / (const Complex& s) const +operator / (const Complex& s, const ColumnVector& a) { - return ComplexColumnVector (divide (data, len, s), len); + int a_len = a.length (); + return ComplexColumnVector (divide (s, a.data (), a_len), a_len); } // column vector by row vector -> matrix operations Matrix -ColumnVector::operator * (const RowVector& a) const +operator * (const ColumnVector& v, const RowVector& a) { - if (len != a.len) + int len = v.length (); + int a_len = a.length (); + if (len != a_len) { (*current_liboctave_error_handler) ("nonconformant vector multiplication attempted"); @@ -397,77 +326,28 @@ double alpha = 1.0; double beta = 0.0; int anr = 1; - int anc = a.len; - double *c = new double [len * a.len]; + double *c = new double [len * a_len]; - F77_FCN (dgemm) (&transa, &transb, &len, &anc, &anr, &alpha, data, - &len, a.data, &anr, &beta, c, &len, 1L, 1L); + 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); + return Matrix (c, len, a_len); } ComplexMatrix -ColumnVector::operator * (const ComplexRowVector& a) const +operator * (const ColumnVector& v, const ComplexRowVector& a) { - ComplexColumnVector tmp (*this); + ComplexColumnVector tmp (v); return tmp * a; } -// column vector by column vector -> column vector operations - -ColumnVector -ColumnVector::operator + (const ColumnVector& a) const -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector addition attempted"); - return ColumnVector (); - } - - if (len == 0) - return ColumnVector (0); - - return ColumnVector (add (data, a.data, len), len); -} - -ColumnVector -ColumnVector::operator - (const ColumnVector& a) const +ComplexColumnVector +operator + (const ColumnVector& v, const ComplexColumnVector& a) { - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector subtraction attempted"); - return ColumnVector (); - } - - if (len == 0) - return ColumnVector (0); - - return ColumnVector (subtract (data, a.data, len), len); -} - -ComplexColumnVector -ColumnVector::operator + (const ComplexColumnVector& a) const -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector addition attempted"); - return ComplexColumnVector (); - } - - if (len == 0) - return ComplexColumnVector (0); - - return ComplexColumnVector (add (data, a.data, len), len); -} - -ComplexColumnVector -ColumnVector::operator - (const ComplexColumnVector& a) const -{ - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector subtraction attempted"); @@ -477,45 +357,31 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (subtract (data, a.data, len), len); + return ComplexColumnVector (add (v.data (), a.data (), len), len); } -ColumnVector -ColumnVector::product (const ColumnVector& a) const +ComplexColumnVector +operator - (const ColumnVector& v, const ComplexColumnVector& a) { - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) - ("nonconformant vector product attempted"); - return ColumnVector (); + ("nonconformant vector subtraction attempted"); + return ComplexColumnVector (); } if (len == 0) - return ColumnVector (0); - - return ColumnVector (multiply (data, a.data, len), len); -} + return ComplexColumnVector (0); -ColumnVector -ColumnVector::quotient (const ColumnVector& a) const -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector quotient attempted"); - return ColumnVector (); - } - - if (len == 0) - return ColumnVector (0); - - return ColumnVector (divide (data, a.data, len), len); + return ComplexColumnVector (subtract (v.data (), a.data (), len), len); } ComplexColumnVector -ColumnVector::product (const ComplexColumnVector& a) const +product (const ColumnVector& v, const ComplexColumnVector& a) { - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector product attempted"); @@ -525,13 +391,14 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (multiply (data, a.data, len), len); + return ComplexColumnVector (multiply (v.data (), a.data (), len), len); } ComplexColumnVector -ColumnVector::quotient (const ComplexColumnVector& a) const +quotient (const ColumnVector& v, const ComplexColumnVector& a) { - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector quotient attempted"); @@ -541,53 +408,10 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (divide (data, a.data, len), len); -} - -ColumnVector& -ColumnVector::operator += (const ColumnVector& a) -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector += operation attempted"); - return ColumnVector (); - } - - if (len == 0) - return *this; - - add2 (data, a.data, len); - return *this; + return ComplexColumnVector (divide (v.data (), a.data (), len), len); } -ColumnVector& -ColumnVector::operator -= (const ColumnVector& a) -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector -= operation attempted"); - return ColumnVector (); - } - - if (len == 0) - return *this; - - subtract2 (data, a.data, len); - return *this; -} - -// unary operations - -ColumnVector -ColumnVector::operator - (void) const -{ - if (len == 0) - return ColumnVector (0); - - return ColumnVector (negate (data, len), len); -} +// other operations ColumnVector map (d_d_Mapper f, const ColumnVector& a) @@ -600,21 +424,22 @@ void ColumnVector::map (d_d_Mapper f) { - for (int i = 0; i < len; i++) - data[i] = f (data[i]); + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); } double ColumnVector::min (void) const { + int len = length (); if (len == 0) return 0.0; - double res = data[0]; + double res = elem (0); for (int i = 1; i < len; i++) - if (data[i] < res) - res = data[i]; + if (elem (i) < res) + res = elem (i); return res; } @@ -622,14 +447,15 @@ double ColumnVector::max (void) const { + int len = length (); if (len == 0) return 0.0; - double res = data[0]; + double res = elem (0); for (int i = 1; i < len; i++) - if (data[i] > res) - res = data[i]; + if (elem (i) > res) + res = elem (i); return res; } @@ -638,8 +464,8 @@ operator << (ostream& os, const ColumnVector& a) { // int field_width = os.precision () + 7; - for (int i = 0; i < a.len; i++) - os << /* setw (field_width) << */ a.data[i] << "\n"; + for (int i = 0; i < a.length (); i++) + os << /* setw (field_width) << */ a.elem (i) << "\n"; return os; } @@ -647,167 +473,14 @@ * Complex Column Vector class */ -ComplexColumnVector::ComplexColumnVector (int n) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create column vector with negative dimension"); - len = 0; - data = (Complex *) NULL; - return; - } - - len = n; - if (n > 0) - data = new Complex [len]; - else - data = (Complex *) NULL; -} - -ComplexColumnVector::ComplexColumnVector (int n, double val) +ComplexColumnVector::ComplexColumnVector (const ColumnVector& a) + : Array<Complex> (a.length ()) { - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create column vector with negative dimension"); - len = 0; - data = (Complex *) NULL; - return; - } - - len = n; - if (n > 0) - { - data = new Complex [len]; - copy (data, len, val); - } - else - data = (Complex *) NULL; -} - -ComplexColumnVector::ComplexColumnVector (int n, const Complex& val) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create column vector with negative dimension"); - len = 0; - data = (Complex *) NULL; - return; - } - - len = n; - if (n > 0) - { - data = new Complex [len]; - copy (data, len, val); - } - else - data = (Complex *) NULL; -} - -ComplexColumnVector::ComplexColumnVector (const ColumnVector& a) -{ - len = a.len; - if (len > 0) - { - data = new Complex [len]; - copy (data, a.data, len); - } - else - data = (Complex *) NULL; + for (int i = 0; i < length (); i++) + elem (i) = a.elem (i); } -ComplexColumnVector::ComplexColumnVector (const ComplexColumnVector& a) -{ - len = a.len; - if (len > 0) - { - data = new Complex [len]; - copy (data, a.data, len); - } - else - data = (Complex *) NULL; -} - -ComplexColumnVector::ComplexColumnVector (double a) -{ - len = 1; - data = new Complex [1]; - data[0] = a; -} - -ComplexColumnVector::ComplexColumnVector (const Complex& a) -{ - len = 1; - data = new Complex [1]; - data[0] = Complex (a); -} - -ComplexColumnVector& -ComplexColumnVector::operator = (const ColumnVector& a) -{ - delete [] data; - len = a.len; - if (len > 0) - { - data = new Complex [len]; - copy (data, a.data, len); - } - else - data = (Complex *) NULL; - - return *this; -} - -ComplexColumnVector& -ComplexColumnVector::operator = (const ComplexColumnVector& a) -{ - if (this != &a) - { - delete [] data; - len = a.len; - if (len > 0) - { - data = new Complex [len]; - copy (data, a.data, len); - } - else - data = (Complex *) NULL; - } - return *this; -} - -Complex& -ComplexColumnVector::checkelem (int n) -{ -#ifndef NO_RANGE_CHECK - if (n < 0 || n >= len) - { - (*current_liboctave_error_handler) ("range error"); - static Complex foo (0.0); - return foo; - } -#endif - - return elem (n); -} - -Complex -ComplexColumnVector::checkelem (int n) const -{ -#ifndef NO_RANGE_CHECK - if (n < 0 || n >= len) - { - (*current_liboctave_error_handler) ("range error"); - return Complex (0.0); - } -#endif - - return elem (n); -} - +#if 0 ComplexColumnVector& ComplexColumnVector::resize (int n) { @@ -856,21 +529,21 @@ return *this; } +#endif int ComplexColumnVector::operator == (const ComplexColumnVector& a) const { - if (len != a.len) + int len = length (); + if (len != a.length ()) return 0; - return equal (data, a.data, len); + return equal (data (), a.data (), len); } int ComplexColumnVector::operator != (const ComplexColumnVector& a) const { - if (len != a.len) - return 0; - return !equal (data, a.data, len); + return !(*this == a); } // destructive insert/delete/reorder operations @@ -878,14 +551,15 @@ ComplexColumnVector& ComplexColumnVector::insert (const ColumnVector& a, int r) { - if (r < 0 || r + a.len - 1 > len) + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) { (*current_liboctave_error_handler) ("range error for insert"); return *this; } - for (int i = 0; i < a.len; i++) - data[r+i] = a.data[i]; + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); return *this; } @@ -893,14 +567,15 @@ ComplexColumnVector& ComplexColumnVector::insert (const ComplexColumnVector& a, int r) { - if (r < 0 || r + a.len - 1 > len) + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) { (*current_liboctave_error_handler) ("range error for insert"); return *this; } - for (int i = 0; i < a.len; i++) - data[r+i] = a.data[i]; + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); return *this; } @@ -908,22 +583,27 @@ ComplexColumnVector& ComplexColumnVector::fill (double val) { + int len = length (); if (len > 0) - copy (data, len, val); + for (int i = 0; i < len; i++) + elem (i) = val; return *this; } ComplexColumnVector& ComplexColumnVector::fill (const Complex& val) { + int len = length (); if (len > 0) - copy (data, len, val); + for (int i = 0; i < len; i++) + elem (i) = val; return *this; } ComplexColumnVector& ComplexColumnVector::fill (double val, int r1, int r2) { + int len = length (); if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) { (*current_liboctave_error_handler) ("range error for fill"); @@ -933,7 +613,7 @@ if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } for (int i = r1; i <= r2; i++) - data[i] = val; + elem (i) = val; return *this; } @@ -941,6 +621,7 @@ ComplexColumnVector& ComplexColumnVector::fill (const Complex& val, int r1, int r2) { + int len = length (); if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) { (*current_liboctave_error_handler) ("range error for fill"); @@ -950,7 +631,7 @@ if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } for (int i = r1; i <= r2; i++) - data[i] = val; + elem (i) = val; return *this; } @@ -958,8 +639,9 @@ ComplexColumnVector ComplexColumnVector::stack (const ColumnVector& a) const { + int len = length (); int nr_insert = len; - ComplexColumnVector retval (len + a.len); + ComplexColumnVector retval (len + a.length ()); retval.insert (*this, 0); retval.insert (a, nr_insert); return retval; @@ -968,8 +650,9 @@ ComplexColumnVector ComplexColumnVector::stack (const ComplexColumnVector& a) const { + int len = length (); int nr_insert = len; - ComplexColumnVector retval (len + a.len); + ComplexColumnVector retval (len + a.length ()); retval.insert (*this, 0); retval.insert (a, nr_insert); return retval; @@ -978,39 +661,44 @@ ComplexRowVector ComplexColumnVector::hermitian (void) const { - return ComplexRowVector (conj_dup (data, len), len); + int len = length (); + return ComplexRowVector (conj_dup (data (), len), len); } ComplexRowVector ComplexColumnVector::transpose (void) const { - return ComplexRowVector (dup (data, len), len); + int len = length (); + return ComplexRowVector (dup (data (), len), len); } ColumnVector real (const ComplexColumnVector& a) { + int a_len = a.length (); ColumnVector retval; - if (a.len > 0) - retval = ColumnVector (real_dup (a.data, a.len), a.len); + 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); + if (a_len > 0) + retval = ColumnVector (imag_dup (a.data (), a_len), a_len); return retval; } ComplexColumnVector conj (const ComplexColumnVector& a) { + int a_len = a.length (); ComplexColumnVector retval; - if (a.len > 0) - retval = ComplexColumnVector (conj_dup (a.data, a.len), a.len); + if (a_len > 0) + retval = ComplexColumnVector (conj_dup (a.data (), a_len), a_len); return retval; } @@ -1026,59 +714,122 @@ ComplexColumnVector result (new_r); for (int i = 0; i < new_r; i++) - result.data[i] = elem (r1+i); + result.elem (i) = elem (r1+i); return result; } +// column vector by column vector -> column vector operations + +ComplexColumnVector& +ComplexColumnVector::operator += (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::operator -= (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::operator += (const ComplexColumnVector& a) +{ + int len = length (); + + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ComplexColumnVector& +ComplexColumnVector::operator -= (const ComplexColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return *this; + } + + if (len == 0) + return *this; + + Complex *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + // column vector by scalar -> column vector operations ComplexColumnVector -ComplexColumnVector::operator + (double s) const -{ - return ComplexColumnVector (add (data, len, s), len); -} - -ComplexColumnVector -ComplexColumnVector::operator - (double s) const +operator + (const ComplexColumnVector& v, double s) { - return ComplexColumnVector (subtract (data, len, s), len); -} - -ComplexColumnVector -ComplexColumnVector::operator * (double s) const -{ - return ComplexColumnVector (multiply (data, len, s), len); + int len = v.length (); + return ComplexColumnVector (add (v.data (), len, s), len); } ComplexColumnVector -ComplexColumnVector::operator / (double s) const +operator - (const ComplexColumnVector& v, double s) { - return ComplexColumnVector (divide (data, len, s), len); -} - -ComplexColumnVector -ComplexColumnVector::operator + (const Complex& s) const -{ - return ComplexColumnVector (add (data, len, s), len); + int len = v.length (); + return ComplexColumnVector (subtract (v.data (), len, s), len); } ComplexColumnVector -ComplexColumnVector::operator - (const Complex& s) const +operator * (const ComplexColumnVector& v, double s) { - return ComplexColumnVector (subtract (data, len, s), len); + int len = v.length (); + return ComplexColumnVector (multiply (v.data (), len, s), len); } ComplexColumnVector -ComplexColumnVector::operator * (const Complex& s) const +operator / (const ComplexColumnVector& v, double s) { - return ComplexColumnVector (multiply (data, len, s), len); -} - -ComplexColumnVector -ComplexColumnVector::operator / (const Complex& s) const -{ - return ComplexColumnVector (divide (data, len, s), len); + int len = v.length (); + return ComplexColumnVector (divide (v.data (), len, s), len); } // scalar by column vector -> column vector operations @@ -1086,64 +837,39 @@ ComplexColumnVector operator + (double s, const ComplexColumnVector& a) { - return ComplexColumnVector (add (a.data, a.len, s), a.len); + int a_len = a.length (); + return ComplexColumnVector (add (a.data (), a_len, s), a_len); } ComplexColumnVector operator - (double s, const ComplexColumnVector& a) { - return ComplexColumnVector (subtract (s, a.data, a.len), a.len); + int a_len = a.length (); + return ComplexColumnVector (subtract (s, a.data (), a_len), a_len); } ComplexColumnVector operator * (double s, const ComplexColumnVector& a) { - return ComplexColumnVector (multiply (a.data, a.len, s), a.len); + int a_len = a.length (); + return ComplexColumnVector (multiply (a.data (), a_len, s), a_len); } ComplexColumnVector operator / (double s, const ComplexColumnVector& a) { - return ComplexColumnVector (divide (s, a.data, a.len), a.len); -} - -ComplexColumnVector -operator + (const Complex& s, const ComplexColumnVector& a) -{ - return ComplexColumnVector (add (a.data, a.len, s), a.len); -} - -ComplexColumnVector -operator - (const Complex& s, const ComplexColumnVector& a) -{ - return ComplexColumnVector (subtract (s, a.data, a.len), a.len); -} - -ComplexColumnVector -operator * (const Complex& s, const ComplexColumnVector& a) -{ - return ComplexColumnVector (multiply (a.data, a.len, s), a.len); -} - -ComplexColumnVector -operator / (const Complex& s, const ComplexColumnVector& a) -{ - return ComplexColumnVector (divide (s, a.data, a.len), a.len); + int a_len = a.length (); + return ComplexColumnVector (divide (s, a.data (), a_len), a_len); } // column vector by row vector -> matrix operations ComplexMatrix -ComplexColumnVector::operator * (const RowVector& a) const +operator * (const ComplexColumnVector& v, const ComplexRowVector& a) { - ComplexRowVector tmp (a); - return *this * tmp; -} - -ComplexMatrix -ComplexColumnVector::operator * (const ComplexRowVector& a) const -{ - if (len != a.len) + int len = v.length (); + int a_len = a.length (); + if (len != a_len) { (*current_liboctave_error_handler) ("nonconformant vector multiplication attempted"); @@ -1158,22 +884,23 @@ Complex alpha (1.0); Complex beta (0.0); int anr = 1; - int anc = a.len; - Complex *c = new Complex [len * a.len]; + Complex *c = new Complex [len * a_len]; - F77_FCN (zgemm) (&transa, &transb, &len, &anc, &anr, &alpha, data, - &len, a.data, &anr, &beta, c, &len, 1L, 1L); + 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); + return ComplexMatrix (c, len, a_len); } // column vector by column vector -> column vector operations ComplexColumnVector -ComplexColumnVector::operator + (const ColumnVector& a) const +operator + (const ComplexColumnVector& v, const ColumnVector& a) { - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector addition attempted"); @@ -1183,13 +910,14 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (add (data, a.data, len), len); + return ComplexColumnVector (add (v.data (), a.data (), len), len); } ComplexColumnVector -ComplexColumnVector::operator - (const ColumnVector& a) const +operator - (const ComplexColumnVector& v, const ColumnVector& a) { - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector subtraction attempted"); @@ -1199,45 +927,14 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (subtract (data, a.data, len), len); -} - -ComplexColumnVector -ComplexColumnVector::operator + (const ComplexColumnVector& a) const -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector addition attempted"); - return ComplexColumnVector (); - } - - if (len == 0) - return ComplexColumnVector (0); - - return ComplexColumnVector (add (data, a.data, len), len); + return ComplexColumnVector (subtract (v.data (), a.data (), len), len); } ComplexColumnVector -ComplexColumnVector::operator - (const ComplexColumnVector& a) const +product (const ComplexColumnVector& v, const ColumnVector& a) { - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector subtraction attempted"); - return ComplexColumnVector (); - } - - if (len == 0) - return ComplexColumnVector (0); - - return ComplexColumnVector (subtract (data, a.data, len), len); -} - -ComplexColumnVector -ComplexColumnVector::product (const ColumnVector& a) const -{ - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector product attempted"); @@ -1247,13 +944,14 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (multiply (data, a.data, len), len); + return ComplexColumnVector (multiply (v.data (), a.data (), len), len); } ComplexColumnVector -ComplexColumnVector::quotient (const ColumnVector& a) const +quotient (const ComplexColumnVector& v, const ColumnVector& a) { - if (len != a.len) + int len = v.length (); + if (len != a.length ()) { (*current_liboctave_error_handler) ("nonconformant vector quotient attempted"); @@ -1263,119 +961,10 @@ if (len == 0) return ComplexColumnVector (0); - return ComplexColumnVector (divide (data, a.data, len), len); -} - -ComplexColumnVector -ComplexColumnVector::product (const ComplexColumnVector& a) const -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector product attempted"); - return ComplexColumnVector (); - } - - if (len == 0) - return ComplexColumnVector (0); - - return ComplexColumnVector (multiply (data, a.data, len), len); -} - -ComplexColumnVector -ComplexColumnVector::quotient (const ComplexColumnVector& a) const -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector quotient attempted"); - return ComplexColumnVector (); - } - - if (len == 0) - return ComplexColumnVector (0); - - return ComplexColumnVector (divide (data, a.data, len), len); -} - -ComplexColumnVector& -ComplexColumnVector::operator += (const ColumnVector& a) -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector += operation attempted"); - return *this; - } - - if (len == 0) - return *this; - - add2 (data, a.data, len); - return *this; + return ComplexColumnVector (divide (v.data (), a.data (), len), len); } -ComplexColumnVector& -ComplexColumnVector::operator -= (const ColumnVector& a) -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector -= operation attempted"); - return *this; - } - - if (len == 0) - return *this; - - subtract2 (data, a.data, len); - return *this; -} - -ComplexColumnVector& -ComplexColumnVector::operator += (const ComplexColumnVector& a) -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector += operation attempted"); - return *this; - } - - if (len == 0) - return *this; - - add2 (data, a.data, len); - return *this; -} - -ComplexColumnVector& -ComplexColumnVector::operator -= (const ComplexColumnVector& a) -{ - if (len != a.len) - { - (*current_liboctave_error_handler) - ("nonconformant vector -= operation attempted"); - return *this; - } - - if (len == 0) - return *this; - - subtract2 (data, a.data, len); - return *this; -} - -// unary operations - -ComplexColumnVector -ComplexColumnVector::operator - (void) const -{ - if (len == 0) - return ComplexColumnVector (0); - - return ComplexColumnVector (negate (data, len), len); -} +// other operations ComplexColumnVector map (c_c_Mapper f, const ComplexColumnVector& a) @@ -1388,8 +977,9 @@ ColumnVector map (d_c_Mapper f, const ComplexColumnVector& a) { - ColumnVector b (a.len); - for (int i = 0; i < a.len; i++) + 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; } @@ -1397,23 +987,24 @@ void ComplexColumnVector::map (c_c_Mapper f) { - for (int i = 0; i < len; i++) - data[i] = f (data[i]); + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); } Complex ComplexColumnVector::min (void) const { + int len = length (); if (len == 0) return 0.0; - Complex res = data[0]; + Complex res = elem (0); double absres = abs (res); for (int i = 1; i < len; i++) - if (abs (data[i]) < absres) + if (abs (elem (i)) < absres) { - res = data[i]; + res = elem (i); absres = abs (res); } @@ -1423,16 +1014,17 @@ Complex ComplexColumnVector::max (void) const { + int len = length (); if (len == 0) return 0.0; - Complex res = data[0]; + Complex res = elem (0); double absres = abs (res); for (int i = 1; i < len; i++) - if (abs (data[i]) > absres) + if (abs (elem (i)) > absres) { - res = data[i]; + res = elem (i); absres = abs (res); } @@ -1445,8 +1037,8 @@ operator << (ostream& os, const ComplexColumnVector& a) { // int field_width = os.precision () + 7; - for (int i = 0; i < a.len; i++) - os << /* setw (field_width) << */ a.data[i] << "\n"; + for (int i = 0; i < a.length (); i++) + os << /* setw (field_width) << */ a.elem (i) << "\n"; return os; }
--- a/liboctave/CollocWt.cc +++ b/liboctave/CollocWt.cc @@ -21,11 +21,12 @@ */ -#ifdef __GNUG__ -#pragma implementation +#ifdef HAVE_CONFIG_H +#include "config.h" #endif #include <iostream.h> + #include "CollocWt.h" #include "f77-uscore.h" #include "lo-error.h"
--- a/liboctave/CollocWt.h +++ b/liboctave/CollocWt.h @@ -24,11 +24,8 @@ #if !defined (_CollocWt_h) #define _CollocWt_h 1 -#ifdef __GNUG__ -#pragma interface -#endif +class ostream; -#include <iostream.h> #include "Matrix.h" #ifndef Vector
--- a/liboctave/DAE.h +++ b/liboctave/DAE.h @@ -24,15 +24,9 @@ #if !defined (_DAE_h) #define _DAE_h 1 -#ifdef __GNUG__ -#pragma interface -#endif - -#include <iostream.h> #include "ODE.h" #include "DAEFunc.h" #include "Matrix.h" -#include "f77-uscore.h" #ifndef Vector #define Vector ColumnVector
--- a/liboctave/DAEFunc.cc +++ b/liboctave/DAEFunc.cc @@ -21,11 +21,10 @@ */ -#ifdef __GNUG__ -#pragma implementation +#ifdef HAVE_CONFIG_H +#include "config.h" #endif -#include <iostream.h> #include "DAEFunc.h" DAEFunc::DAEFunc (void)
--- a/liboctave/DAEFunc.h +++ b/liboctave/DAEFunc.h @@ -24,11 +24,6 @@ #if !defined (_DAEFunc_h) #define _DAEFunc_h 1 -#ifdef __GNUG__ -#pragma interface -#endif - -#include <iostream.h> #include "Matrix.h" #ifndef Vector
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -21,12 +21,12 @@ */ -#ifdef __GNUG__ -#pragma implementation +#ifdef HAVE_CONFIG_H +#include "config.h" #endif -#include <iostream.h> #include "DAE.h" +#include "f77-uscore.h" #include "lo-error.h" extern "C" @@ -232,7 +232,7 @@ // Fix up the matrix of partial derivatives for dassl. - tmp_dfdx = tmp_dfdx + (*cj * tmp_dfdxdot); + tmp_dfdx = tmp_dfdx + (tmp_dfdxdot * (*cj)); for (int j = 0; j < nn; j++) for (int i = 0; i < nn; i++)
--- a/liboctave/DiagMatrix.cc +++ b/liboctave/DiagMatrix.cc @@ -21,12 +21,11 @@ */ -// I\'m not sure how this is supposed to work if the .h file declares -// several classes, each of which is defined in a separate file... -// -// #ifdef __GNUG__ -// #pragma implementation "Matrix.h" -// #endif +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include <iostream.h> #include "Matrix.h" #include "mx-inlines.cc" @@ -36,200 +35,7 @@ * Diagonal Matrix class. */ -DiagMatrix::DiagMatrix (int n) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create matrix with negative dimensions"); - nr = 0; - nc = 0; - len = 0; - data = (double *) NULL; - return; - } - - nr = n; - nc = n; - len = n; - if (len > 0) - data = new double [len]; - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (int n, double val) -{ - if (n < 0) - { - (*current_liboctave_error_handler) - ("can't create matrix with negative dimensions"); - nr = 0; - nc = 0; - len = 0; - data = (double *) NULL; - return; - } - - nr = n; - nc = n; - len = n; - if (len > 0) - { - data = new double [len]; - copy (data, len, val); - } - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (int r, int c) -{ - if (r < 0 || c < 0) - { - (*current_liboctave_error_handler) - ("can't create matrix with negative dimensions"); - nr = 0; - nc = 0; - len = 0; - data = (double *) NULL; - return; - } - - nr = r; - nc = c; - len = r < c ? r : c; - if (len > 0) - data = new double [len]; - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (int r, int c, double val) -{ - if (r < 0 || c < 0) - { - (*current_liboctave_error_handler) - ("can't create matrix with negative dimensions"); - nr = 0; - nc = 0; - len = 0; - data = (double *) NULL; - return; - } - - nr = r; - nc = c; - len = r < c ? r : c; - if (len > 0) - { - data = new double [len]; - copy (data, len, val); - } - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (const RowVector& a) -{ - nr = a.len; - nc = nr; - len = nr; - if (len > 0) - { - data = new double [len]; - copy (data, a.data, len); - } - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (const ColumnVector& a) -{ - nr = a.len; - nc = nr; - len = nr; - if (len > 0) - { - data = new double [len]; - copy (data, a.data, len); - } - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (const DiagMatrix& a) -{ - nr = a.nr; - nc = a.nc; - len = a.len; - if (len > 0) - { - data = new double [len]; - copy (data, a.data, len); - } - else - data = (double *) NULL; -} - -DiagMatrix::DiagMatrix (double a) -{ - nr = 1; - nc = 1; - len = 1; - data = new double [1]; - data[0] = a; -} - -DiagMatrix& -DiagMatrix::operator = (const DiagMatrix& a) -{ - if (this != &a) - { - delete [] data; - nr = a.nr; - nc = a.nc; - len = a.len; - if (len > 0) - { - data = new double [len]; - copy (data, a.data, len); - } - else - data = (double *) NULL; - } - return *this; -} - -double& -DiagMatrix::checkelem (int r, int c) -{ -#ifndef NO_RANGE_CHECK - if (r < 0 || r >= nr || c < 0 || c >= nc) - { - (*current_liboctave_error_handler) ("range error"); - static double foo = 0.0; - return foo; - } -#endif - - return elem (r, c); -} - -double -DiagMatrix::checkelem (int r, int c) const -{ -#ifndef NO_RANGE_CHECK - if (r < 0 || r >= nr || c < 0 || c >= nc) - { - (*current_liboctave_error_handler) ("range error"); - return 0.0; - } -#endif - - return elem (r, c); -} - +#if 0 DiagMatrix& DiagMatrix::resize (int r, int c) { @@ -294,102 +100,114 @@ return *this; } +#endif int DiagMatrix::operator == (const DiagMatrix& a) const { - if (nr != a.nr || nc != a.nc) + if (rows () != a.rows () || cols () != a.cols ()) return 0; - return equal (data, a.data, len); + return equal (data (), a.data (), length ()); } int DiagMatrix::operator != (const DiagMatrix& a) const { - if (nr != a.nr || nc != a.nc) - return 1; - - return !equal (data, a.data, len); + return !(*this == a); } DiagMatrix& DiagMatrix::fill (double val) { - copy (data, len, val); + for (int i = 0; i < length (); i++) + elem (i, i) = val; return *this; } DiagMatrix& DiagMatrix::fill (double val, int beg, int end) { - if (beg < 0 || end >= len || end < beg) + if (beg < 0 || end >= length () || end < beg) { (*current_liboctave_error_handler) ("range error for fill"); return *this; } - if (end > beg) - copy (data+beg, beg-end, val); + for (int i = beg; i < end; i++) + elem (i, i) = val; + return *this; } DiagMatrix& DiagMatrix::fill (const ColumnVector& a) { - if (a.len != len) + int len = length (); + if (a.length () != len) { (*current_liboctave_error_handler) ("range error for fill"); return *this; } - copy (data, a.data, len); + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + return *this; } DiagMatrix& DiagMatrix::fill (const RowVector& a) { - if (a.len != len) + int len = length (); + if (a.length () != len) { (*current_liboctave_error_handler) ("range error for fill"); return *this; } - copy (data, a.data, len); + for (int i = 0; i < len; i++) + elem (i, i) = a.elem (i); + return *this; } DiagMatrix& DiagMatrix::fill (const ColumnVector& a, int beg) { - if (beg < 0 || beg + a.len >= len) + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) { (*current_liboctave_error_handler) ("range error for fill"); return *this; } - copy (data+beg, a.data, a.len); + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + return *this; } DiagMatrix& DiagMatrix::fill (const RowVector& a, int beg) { - if (beg < 0 || beg + a.len >= len) + int a_len = a.length (); + if (beg < 0 || beg + a_len >= length ()) { (*current_liboctave_error_handler) ("range error for fill"); return *this; } - copy (data+beg, a.data, a.len); + for (int i = 0; i < a_len; i++) + elem (i+beg, i+beg) = a.elem (i); + return *this; } DiagMatrix DiagMatrix::transpose (void) const { - return DiagMatrix (dup (data, len), nc, nr); + return DiagMatrix (dup (data (), length ()), cols (), rows ()); } Matrix @@ -405,7 +223,7 @@ for (int j = 0; j < new_c; j++) for (int i = 0; i < new_r; i++) - result.data[new_r*j+i] = elem (r1+i, c1+j); + result.elem (i, j) = elem (r1+i, c1+j); return result; } @@ -415,6 +233,8 @@ RowVector DiagMatrix::row (int i) const { + int nr = rows (); + int nc = cols (); if (i < 0 || i >= nr) { (*current_liboctave_error_handler) ("invalid row selection"); @@ -423,7 +243,7 @@ RowVector retval (nc, 0.0); if (nr <= nc || (nr > nc && i < nc)) - retval.data [i] = data[i]; + retval.elem (i) = elem (i, i); return retval; } @@ -441,7 +261,7 @@ if (c == 'f' || c == 'F') return row (0); else if (c == 'l' || c == 'L') - return row (nr - 1); + return row (rows () - 1); else { (*current_liboctave_error_handler) ("invalid row selection"); @@ -452,6 +272,8 @@ ColumnVector DiagMatrix::column (int i) const { + int nr = rows (); + int nc = cols (); if (i < 0 || i >= nc) { (*current_liboctave_error_handler) ("invalid column selection"); @@ -460,7 +282,7 @@ ColumnVector retval (nr, 0.0); if (nr >= nc || (nr < nc && i < nr)) - retval.data [i] = data[i]; + retval.elem (i) = elem (i, i); return retval; } @@ -478,7 +300,7 @@ if (c == 'f' || c == 'F') return column (0); else if (c == 'l' || c == 'L') - return column (nc - 1); + return column (cols () - 1); else { (*current_liboctave_error_handler) ("invalid column selection"); @@ -487,94 +309,131 @@ } DiagMatrix -DiagMatrix::inverse (int &info) const -{ - if (nr != nc) - { - (*current_liboctave_error_handler) ("inverse requires square matrix"); - return DiagMatrix (); - } - - info = 0; - double *tmp_data = dup (data, len); - for (int i = 0; i < len; i++) - { - if (data[i] == 0.0) - { - info = -1; - copy (tmp_data, data, len); // Restore contents. - break; - } - else - { - tmp_data[i] = 1.0 / data[i]; - } - } - - return DiagMatrix (tmp_data, nr, nc); -} - -DiagMatrix DiagMatrix::inverse (void) const { int info; return inverse (info); } +DiagMatrix +DiagMatrix::inverse (int &info) const +{ + int nr = rows (); + int nc = cols (); + int len = length (); + if (nr != nc) + { + (*current_liboctave_error_handler) ("inverse requires square matrix"); + return DiagMatrix (); + } + + info = 0; + double *tmp_data = dup (data (), len); + for (int i = 0; i < len; i++) + { + if (elem (i, i) == 0.0) + { + info = -1; + copy (tmp_data, data (), len); // Restore contents. + break; + } + else + { + tmp_data[i] = 1.0 / elem (i, i); + } + } + + return DiagMatrix (tmp_data, nr, nc); +} + +// diagonal matrix by diagonal matrix -> diagonal matrix operations + +DiagMatrix& +DiagMatrix::operator += (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix += operation attempted"); + return *this; + } + + if (nc == 0 || nr == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), length ()); + return *this; +} + +DiagMatrix& +DiagMatrix::operator -= (const DiagMatrix& a) +{ + int nr = rows (); + int nc = cols (); + if (nr != a.rows () || nc != a.cols ()) + { + (*current_liboctave_error_handler) + ("nonconformant matrix -= operation attempted"); + return *this; + } + + if (nr == 0 || nc == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), length ()); + return *this; +} + // diagonal matrix by scalar -> matrix operations Matrix -DiagMatrix::operator + (double s) const +operator + (const DiagMatrix& a, double s) { - Matrix tmp (nr, nc, s); - return *this + tmp; + Matrix tmp (a.rows (), a.cols (), s); + return a + tmp; } Matrix -DiagMatrix::operator - (double s) const +operator - (const DiagMatrix& a, double s) { - Matrix tmp (nr, nc, -s); - return *this + tmp; + Matrix tmp (a.rows (), a.cols (), -s); + return a + tmp; } ComplexMatrix -DiagMatrix::operator + (const Complex& s) const +operator + (const DiagMatrix& a, const Complex& s) { - ComplexMatrix tmp (nr, nc, s); - return *this + tmp; + ComplexMatrix tmp (a.rows (), a.cols (), s); + return a + tmp; } ComplexMatrix -DiagMatrix::operator - (const Complex& s) const +operator - (const DiagMatrix& a, const Complex& s) { - ComplexMatrix tmp (nr, nc, -s); - return *this + tmp; + ComplexMatrix tmp (a.rows (), a.cols (), -s); + return a + tmp; } // diagonal matrix by scalar -> diagonal matrix operations -DiagMatrix -DiagMatrix::operator * (double s) const +ComplexDiagMatrix +operator * (const DiagMatrix& a, const Complex& s) { - return DiagMatrix (multiply (data, len, s), nr, nc); -} - -DiagMatrix -DiagMatrix::operator / (double s) const -{ - return DiagMatrix (divide (data, len, s), nr, nc); + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); } ComplexDiagMatrix -DiagMatrix::operator * (const Complex& s) const +operator / (const DiagMatrix& a, const Complex& s) { - return ComplexDiagMatrix (multiply (data, len, s), nr, nc); -} - -ComplexDiagMatrix -DiagMatrix::operator / (const Complex& s) const -{ - return ComplexDiagMatrix (divide (data, len, s), nr, nc); + return ComplexDiagMatrix (divide (a.data (), a.length (), s), + a.rows (), a.cols ()); } // scalar by diagonal matrix -> matrix operations @@ -582,35 +441,49 @@ Matrix operator + (double s, const DiagMatrix& a) { - return a + s; + Matrix tmp (a.rows (), a.cols (), s); + return tmp + a; } Matrix operator - (double s, const DiagMatrix& a) { - return -a + s; + 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 -DiagMatrix -operator * (double s, const DiagMatrix& a) +ComplexDiagMatrix +operator * (const Complex& s, const DiagMatrix& a) { - return DiagMatrix (multiply (a.data, a.len, s), a.nr, a.nc); -} - -DiagMatrix -operator / (double s, const DiagMatrix& a) -{ - return DiagMatrix (divide (s, a.data, a.len), a.nr, a.nc); + return ComplexDiagMatrix (multiply (a.data (), a.length (), s), + a.rows (), a.cols ()); } // diagonal matrix by column vector -> column vector operations ColumnVector -DiagMatrix::operator * (const ColumnVector& a) const +operator * (const DiagMatrix& m, const ColumnVector& a) { - if (nc != a.len) + 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"); @@ -622,19 +495,22 @@ ColumnVector result (nr); - for (int i = 0; i < a.len; i++) - result.data[i] = a.data[i] * data[i]; + 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.data[i] = 0.0; + for (i = a_len; i < nr; i++) + result.elem (i) = 0.0; return result; } ComplexColumnVector -DiagMatrix::operator * (const ComplexColumnVector& a) const +operator * (const DiagMatrix& m, const ComplexColumnVector& a) { - if (nc != a.len) + 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"); @@ -646,69 +522,23 @@ ComplexColumnVector result (nr); - for (int i = 0; i < a.len; i++) - result.data[i] = a.data[i] * data[i]; + 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.data[i] = 0.0; + for (i = a_len; i < nr; i++) + result.elem (i) = 0.0; return result; } // diagonal matrix by diagonal matrix -> diagonal matrix operations -DiagMatrix -DiagMatrix::operator + (const DiagMatrix& a) const -{ - if (nr != a.nr || nc != a.nc) - { - (*current_liboctave_error_handler) - ("nonconformant matrix addition attempted"); - return DiagMatrix (); - } - - if (nc == 0 || nr == 0) - return DiagMatrix (nr, nc); - - return DiagMatrix (add (data, a.data, len), nr , nc); -} - -DiagMatrix -DiagMatrix::operator - (const DiagMatrix& a) const +ComplexDiagMatrix +operator + (const DiagMatrix& m, const ComplexDiagMatrix& a) { - if (nr != a.nr || nc != a.nc) - { - (*current_liboctave_error_handler) - ("nonconformant matrix subtraction attempted"); - return DiagMatrix (); - } - - if (nc == 0 || nr == 0) - return DiagMatrix (nr, nc); - - return DiagMatrix (subtract (data, a.data, len), nr, nc); -} - -DiagMatrix -DiagMatrix::operator * (const DiagMatrix& a) const -{ - if (nr != a.nr || nc != a.nc) - { - (*current_liboctave_error_handler) - ("nonconformant matrix multiplication attempted"); - return DiagMatrix (); - } - - if (nc == 0 || nr == 0) - return DiagMatrix (nr, nc); - - return DiagMatrix (multiply (data, a.data, len), nr, nc); -} - -ComplexDiagMatrix -DiagMatrix::operator + (const ComplexDiagMatrix& a) const -{ - if (nr != a.nr || nc != a.nc)