Mercurial > hg > octave-nkf
diff liboctave/fColVector.cc @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | |
children | eb63fbe60fab |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/liboctave/fColVector.cc @@ -0,0 +1,345 @@ +// ColumnVector manipulations. +/* + +Copyright (C) 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2004, + 2005, 2007 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 3 of the License, 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, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <iostream> + +#include "Array-util.h" +#include "f77-fcn.h" +#include "functor.h" +#include "lo-error.h" +#include "mx-base.h" +#include "mx-inlines.cc" +#include "oct-cmplx.h" + +// Fortran functions we call. + +extern "C" +{ + F77_RET_T + F77_FUNC (sgemv, SGEMV) (F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const float&, + const float*, const octave_idx_type&, const float*, + const octave_idx_type&, const float&, float*, + const octave_idx_type& + F77_CHAR_ARG_LEN_DECL); +} + +// Column Vector class. + +bool +FloatColumnVector::operator == (const FloatColumnVector& a) const +{ + octave_idx_type len = length (); + if (len != a.length ()) + return 0; + return mx_inline_equal (data (), a.data (), len); +} + +bool +FloatColumnVector::operator != (const FloatColumnVector& a) const +{ + return !(*this == a); +} + +FloatColumnVector& +FloatColumnVector::insert (const FloatColumnVector& a, octave_idx_type r) +{ + octave_idx_type a_len = a.length (); + + if (r < 0 || r + a_len > length ()) + { + (*current_liboctave_error_handler) ("range error for insert"); + return *this; + } + + if (a_len > 0) + { + make_unique (); + + for (octave_idx_type i = 0; i < a_len; i++) + xelem (r+i) = a.elem (i); + } + + return *this; +} + +FloatColumnVector& +FloatColumnVector::fill (float val) +{ + octave_idx_type len = length (); + + if (len > 0) + { + make_unique (); + + for (octave_idx_type i = 0; i < len; i++) + xelem (i) = val; + } + + return *this; +} + +FloatColumnVector& +FloatColumnVector::fill (float val, octave_idx_type r1, octave_idx_type r2) +{ + octave_idx_type len = length (); + + if (r1 < 0 || r2 < 0 || r1 >= len || r2 >= len) + { + (*current_liboctave_error_handler) ("range error for fill"); + return *this; + } + + if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } + + if (r2 >= r1) + { + make_unique (); + + for (octave_idx_type i = r1; i <= r2; i++) + xelem (i) = val; + } + + return *this; +} + +FloatColumnVector +FloatColumnVector::stack (const FloatColumnVector& a) const +{ + octave_idx_type len = length (); + octave_idx_type nr_insert = len; + FloatColumnVector retval (len + a.length ()); + retval.insert (*this, 0); + retval.insert (a, nr_insert); + return retval; +} + +FloatRowVector +FloatColumnVector::transpose (void) const +{ + return MArray<float>::transpose(); +} + +FloatColumnVector +real (const FloatComplexColumnVector& a) +{ + octave_idx_type a_len = a.length (); + FloatColumnVector retval; + if (a_len > 0) + retval = FloatColumnVector (mx_inline_real_dup (a.data (), a_len), a_len); + return retval; +} + +FloatColumnVector +imag (const FloatComplexColumnVector& a) +{ + octave_idx_type a_len = a.length (); + FloatColumnVector retval; + if (a_len > 0) + retval = FloatColumnVector (mx_inline_imag_dup (a.data (), a_len), a_len); + return retval; +} + +// resize is the destructive equivalent for this one + +FloatColumnVector +FloatColumnVector::extract (octave_idx_type r1, octave_idx_type r2) const +{ + if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } + + octave_idx_type new_r = r2 - r1 + 1; + + FloatColumnVector result (new_r); + + for (octave_idx_type i = 0; i < new_r; i++) + result.xelem (i) = elem (r1+i); + + return result; +} + +FloatColumnVector +FloatColumnVector::extract_n (octave_idx_type r1, octave_idx_type n) const +{ + FloatColumnVector result (n); + + for (octave_idx_type i = 0; i < n; i++) + result.xelem (i) = elem (r1+i); + + return result; +} + +// matrix by column vector -> column vector operations + +FloatColumnVector +operator * (const FloatMatrix& m, const FloatColumnVector& a) +{ + FloatColumnVector retval; + + octave_idx_type nr = m.rows (); + octave_idx_type nc = m.cols (); + + octave_idx_type a_len = a.length (); + + if (nc != a_len) + gripe_nonconformant ("operator *", nr, nc, a_len, 1); + else + { + if (nr == 0 || nc == 0) + retval.resize (nr, 0.0); + else + { + octave_idx_type ld = nr; + + retval.resize (nr); + float *y = retval.fortran_vec (); + + F77_XFCN (sgemv, SGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), + nr, nc, 1.0, m.data (), ld, + a.data (), 1, 0.0, y, 1 + F77_CHAR_ARG_LEN (1))); + } + } + + return retval; +} + +// diagonal matrix by column vector -> column vector operations + +FloatColumnVector +operator * (const FloatDiagMatrix& m, const FloatColumnVector& a) +{ + FloatColumnVector retval; + + octave_idx_type nr = m.rows (); + octave_idx_type nc = m.cols (); + + octave_idx_type a_len = a.length (); + + if (nc != a_len) + gripe_nonconformant ("operator *", nr, nc, a_len, 1); + else + { + if (nr == 0 || nc == 0) + retval.resize (nr, 0.0); + else + { + retval.resize (nr); + + for (octave_idx_type i = 0; i < a_len; i++) + retval.elem (i) = a.elem (i) * m.elem (i, i); + + for (octave_idx_type i = a_len; i < nr; i++) + retval.elem (i) = 0.0; + } + } + + return retval; +} + +// other operations + +FloatColumnVector +FloatColumnVector::map (dmapper fcn) const +{ + return MArray<float>::map<float> (func_ptr (fcn)); +} + +FloatComplexColumnVector +FloatColumnVector::map (cmapper fcn) const +{ + return MArray<float>::map<FloatComplex> (func_ptr (fcn)); +} + +float +FloatColumnVector::min (void) const +{ + octave_idx_type len = length (); + if (len == 0) + return 0.0; + + float res = elem (0); + + for (octave_idx_type i = 1; i < len; i++) + if (elem (i) < res) + res = elem (i); + + return res; +} + +float +FloatColumnVector::max (void) const +{ + octave_idx_type len = length (); + if (len == 0) + return 0.0; + + float res = elem (0); + + for (octave_idx_type i = 1; i < len; i++) + if (elem (i) > res) + res = elem (i); + + return res; +} + +std::ostream& +operator << (std::ostream& os, const FloatColumnVector& a) +{ +// int field_width = os.precision () + 7; + for (octave_idx_type i = 0; i < a.length (); i++) + os << /* setw (field_width) << */ a.elem (i) << "\n"; + return os; +} + +std::istream& +operator >> (std::istream& is, FloatColumnVector& a) +{ + octave_idx_type len = a.length(); + + if (len < 1) + is.clear (std::ios::badbit); + else + { + float tmp; + for (octave_idx_type i = 0; i < len; i++) + { + is >> tmp; + if (is) + a.elem (i) = tmp; + else + break; + } + } + return is; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/