Mercurial > hg > octave-nkf
diff liboctave/dColVector.cc @ 458:38cb88095913
[project @ 1994-06-06 00:41:10 by jwe]
Initial revision
author | jwe |
---|---|
date | Mon, 06 Jun 1994 00:41:10 +0000 |
parents | |
children | 2ca256b77602 |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/liboctave/dColVector.cc @@ -0,0 +1,491 @@ +// ColumnVector manipulations. -*- C++ -*- +/* + +Copyright (C) 1992, 1993, 1994 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#if defined (__GNUG__) +#pragma implementation +#endif + +#include <iostream.h> + +#include <Complex.h> + +#include "mx-base.h" +#include "mx-inlines.cc" +#include "f77-uscore.h" +#include "lo-error.h" + +// Fortran functions we call. + +extern "C" +{ + int F77_FCN (dgemm) (const char*, const char*, const int*, + const int*, const int*, const double*, + const double*, const int*, const double*, + const int*, const double*, double*, const int*, + long, long); +} + +/* + * Column Vector class. + */ + +#define KLUDGE_VECTORS +#define TYPE double +#define KL_VEC_TYPE ColumnVector +#include "mx-kludge.cc" +#undef KLUDGE_VECTORS +#undef TYPE +#undef KL_VEC_TYPE + +#if 0 +ColumnVector& +ColumnVector::resize (int n) +{ + if (n < 0) + { + (*current_liboctave_error_handler) + ("can't resize to negative dimension"); + return *this; + } + + double *new_data = (double *) NULL; + if (n > 0) + { + new_data = new double [n]; + int min_len = len < n ? len : n; + + for (int i = 0; i < min_len; i++) + new_data[i] = data[i]; + } + + delete [] data; + len = n; + data = new_data; + + return *this; +} + +ColumnVector& +ColumnVector::resize (int n, double val) +{ + int old_len = len; + resize (n); + for (int i = old_len; i < len; i++) + data[i] = val; + + return *this; +} +#endif + +int +ColumnVector::operator == (const ColumnVector& a) const +{ + int len = length (); + if (len != a.length ()) + return 0; + return equal (data (), a.data (), len); +} + +int +ColumnVector::operator != (const ColumnVector& a) const +{ + return !(*this == a); +} + +ColumnVector& +ColumnVector::insert (const ColumnVector& a, int r) +{ + int a_len = a.length (); + if (r < 0 || r + a_len - 1 > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + for (int i = 0; i < a_len; i++) + elem (r+i) = a.elem (i); + + return *this; +} + +ColumnVector& +ColumnVector::fill (double val) +{ + int len = length (); + if (len > 0) + for (int i = 0; i < len; i++) + elem (i) = val; + return *this; +} + +ColumnVector& +ColumnVector::fill (double val, int r1, int r2) +{ + int len = length (); + if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + for (int i = r1; i <= r2; i++) + elem (i) = val; + + return *this; +} + +ColumnVector +ColumnVector::stack (const ColumnVector& a) const +{ + int len = length (); + int nr_insert = len; + ColumnVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nr_insert); + return retval; +} + +RowVector +ColumnVector::transpose (void) const +{ + int len = length (); + return RowVector (dup (data (), len), len); +} + +// resize is the destructive equivalent for this one + +ColumnVector +ColumnVector::extract (int r1, int r2) const +{ + if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; } + + int new_r = r2 - r1 + 1; + + ColumnVector result (new_r); + + for (int i = 0; i < new_r; i++) + result.elem (i) = elem (r1+i); + + return result; +} + +// column vector by column vector -> column vector operations + +ColumnVector& +ColumnVector::operator += (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector += operation attempted"); + return ColumnVector (); + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + add2 (d, a.data (), len); + return *this; +} + +ColumnVector& +ColumnVector::operator -= (const ColumnVector& a) +{ + int len = length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector -= operation attempted"); + return ColumnVector (); + } + + if (len == 0) + return *this; + + double *d = fortran_vec (); // Ensures only one reference to my privates! + + subtract2 (d, a.data (), len); + return *this; +} + +// scalar by column vector -> column vector operations + +ComplexColumnVector +operator + (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (add (a.data (), len, s), len); +} + +ComplexColumnVector +operator - (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (subtract (a.data (), len, s), len); +} + +ComplexColumnVector +operator * (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (multiply (a.data (), len, s), len); +} + +ComplexColumnVector +operator / (const ColumnVector& a, const Complex& s) +{ + int len = a.length (); + return ComplexColumnVector (divide (a.data (), len, s), len); +} + +// scalar by column vector -> column vector operations + +ComplexColumnVector +operator + (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (add (a.data (), a_len, s), a_len); +} + +ComplexColumnVector +operator - (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (subtract (s, a.data (), a_len), a_len); +} + +ComplexColumnVector +operator * (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (multiply (a.data (), a_len, s), a_len); +} + +ComplexColumnVector +operator / (const Complex& s, const ColumnVector& a) +{ + int a_len = a.length (); + return ComplexColumnVector (divide (s, a.data (), a_len), a_len); +} + +// column vector by row vector -> matrix operations + +Matrix +operator * (const ColumnVector& v, const RowVector& a) +{ + int len = v.length (); + int a_len = a.length (); + if (len != a_len) + { + (*current_liboctave_error_handler) + ("nonconformant vector multiplication attempted"); + return Matrix (); + } + + if (len == 0) + return Matrix (len, len, 0.0); + + char transa = 'N'; + char transb = 'N'; + double alpha = 1.0; + double beta = 0.0; + int anr = 1; + + double *c = new double [len * a_len]; + + F77_FCN (dgemm) (&transa, &transb, &len, &a_len, &anr, &alpha, + v.data (), &len, a.data (), &anr, &beta, c, &len, + 1L, 1L); + + return Matrix (c, len, a_len); +} + +ComplexMatrix +operator * (const ColumnVector& v, const ComplexRowVector& a) +{ + ComplexColumnVector tmp (v); + return tmp * a; +} + +ComplexColumnVector +operator + (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (add (v.data (), a.data (), len), len); +} + +ComplexColumnVector +operator - (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector subtraction attempted"); + return ComplexColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (subtract (v.data (), a.data (), len), len); +} + +ComplexColumnVector +product (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector product attempted"); + return ColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (multiply (v.data (), a.data (), len), len); +} + +ComplexColumnVector +quotient (const ColumnVector& v, const ComplexColumnVector& a) +{ + int len = v.length (); + if (len != a.length ()) + { + (*current_liboctave_error_handler) + ("nonconformant vector quotient attempted"); + return ColumnVector (); + } + + if (len == 0) + return ComplexColumnVector (0); + + return ComplexColumnVector (divide (v.data (), a.data (), len), len); +} + +// other operations + +ColumnVector +map (d_d_Mapper f, const ColumnVector& a) +{ + ColumnVector b (a); + b.map (f); + return b; +} + +void +ColumnVector::map (d_d_Mapper f) +{ + for (int i = 0; i < length (); i++) + elem (i) = f (elem (i)); +} + +double +ColumnVector::min (void) const +{ + int len = length (); + if (len == 0) + return 0.0; + + double res = elem (0); + + for (int i = 1; i < len; i++) + if (elem (i) < res) + res = elem (i); + + return res; +} + +double +ColumnVector::max (void) const +{ + int len = length (); + if (len == 0) + return 0.0; + + double res = elem (0); + + for (int i = 1; i < len; i++) + if (elem (i) > res) + res = elem (i); + + return res; +} + +ostream& +operator << (ostream& os, const ColumnVector& a) +{ +// int field_width = os.precision () + 7; + for (int i = 0; i < a.length (); i++) + os << /* setw (field_width) << */ a.elem (i) << "\n"; + return os; +} + +istream& +operator >> (istream& is, ColumnVector& a) +{ + int len = a.length(); + + if (len < 1) + is.clear (ios::badbit); + else + { + double tmp; + for (int i = 0; i < len; i++) + { + is >> tmp; + if (is) + a.elem (i) = tmp; + else + break; + } + } +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/