Mercurial > hg > octave-avbm
changeset 9653:e087d7c77ff9
improve linspace in liboctave
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Fri, 18 Sep 2009 15:27:09 +0200 |
parents | ecdb275bd41b |
children | a307a6f77fb3 |
files | liboctave/Array.h liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/CRowVector.cc liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/dRowVector.cc liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fCRowVector.cc liboctave/fMatrix.cc liboctave/fMatrix.h liboctave/fRowVector.cc |
diffstat | 14 files changed, 231 insertions(+), 60 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/Array.h +++ b/liboctave/Array.h @@ -679,6 +679,39 @@ static void instantiation_guard (); }; +// This is a simple wrapper template that will subclass an Array<T> type or any +// later type derived from it and override the default non-const operator() to +// not check for the array's uniqueness. It is, however, the user's +// responsibility to ensure the array is actually unaliased whenever elements +// are accessed. + +template<class ArrayClass> +class NoAlias : public ArrayClass +{ + typedef typename ArrayClass::element_type T; +public: + NoAlias () : ArrayClass () { } + + // FIXME: this would be simpler once C++0x is available + template <class X> + explicit NoAlias (X x) : ArrayClass (x) { } + + template <class X, class Y> + explicit NoAlias (X x, Y y) : ArrayClass (x, y) { } + + template <class X, class Y, class Z> + explicit NoAlias (X x, Y y, Z z) : ArrayClass (x, y, z) { } + + T& operator () (octave_idx_type n) + { return ArrayClass::xelem (n); } + T& operator () (octave_idx_type i, octave_idx_type j) + { return ArrayClass::xelem (i, j); } + T& operator () (octave_idx_type i, octave_idx_type j, octave_idx_type k) + { return ArrayClass::xelem (i, j, k); } + T& operator () (const Array<octave_idx_type>& ra_idx) + { return ArrayClass::xelem (ra_idx); } +}; + #endif /*
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -4066,6 +4066,39 @@ return result; } +ComplexMatrix linspace (const ComplexColumnVector& x1, + const ComplexColumnVector& x2, + octave_idx_type n) + +{ + if (n < 1) n = 1; + + octave_idx_type m = x1.length (); + + if (x2.length () != m) + (*current_liboctave_error_handler) ("linspace: vectors must be of equal length"); + + NoAlias<ComplexMatrix> retval; + + retval.clear (m, n); + for (octave_idx_type i = 0; i < m; i++) + retval(i, 0) = x1(i); + + // The last column is not needed while using delta. + Complex *delta = &retval(0, 1); + for (octave_idx_type i = 0; i < m; i++) + delta[i] = (x2(i) - x1(i)) / (n - 1.0); + + for (octave_idx_type j = 1; j < n-1; j++) + for (octave_idx_type i = 0; i < m; i++) + retval(i, j) = retval(i, j-1) + delta[i]; + + for (octave_idx_type i = 0; i < m; i++) + retval(i, n-1) = x2(i); + + return retval; +} + MS_CMP_OPS (ComplexMatrix, Complex) MS_BOOL_OPS (ComplexMatrix, Complex)
--- a/liboctave/CMatrix.h +++ b/liboctave/CMatrix.h @@ -408,6 +408,11 @@ extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& m, const Complex& c); extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& a, const ComplexMatrix& b); +extern OCTAVE_API ComplexMatrix linspace (const ComplexColumnVector& x1, + const ComplexColumnVector& x2, + octave_idx_type n); + + MS_CMP_OP_DECLS (ComplexMatrix, Complex, OCTAVE_API) MS_BOOL_OP_DECLS (ComplexMatrix, Complex, OCTAVE_API)
--- a/liboctave/CRowVector.cc +++ b/liboctave/CRowVector.cc @@ -482,22 +482,15 @@ ComplexRowVector linspace (const Complex& x1, const Complex& x2, octave_idx_type n) { - ComplexRowVector retval; + if (n < 1) n = 1; + + NoAlias<ComplexRowVector> retval (n); - if (n > 0) - { - retval.resize (n); - Complex delta = (x2 - x1) / (n - 1.0); - retval.elem (0) = x1; - for (octave_idx_type i = 1; i < n-1; i++) - retval.elem (i) = x1 + 1.0 * i * delta; - retval.elem (n-1) = x2; - } - else - { - retval.resize (1); - retval.elem (0) = x2; - } + Complex delta = (x2 - x1) / (n - 1.0); + Complex y = retval(0) = x1; + for (octave_idx_type i = 1; i < n-1; i++) + retval(i) = y += delta; + retval(n-1) = x2; return retval; }
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,19 @@ +2009-09-18 Jaroslav Hajek <highegg@gmail.com> + + * Array.h (NoAlias): New template class. + * dRowVector.cc (linspace): Rewrite. + * fRowVector.cc (linspace): Rewrite. + * CRowVector.cc (linspace): Rewrite. + * fCRowVector.cc (linspace): Rewrite. + * dMatrix.cc (linspace): New method. + * dMatrix.h (linspace): Declare it. + * fMatrix.cc (linspace): New method. + * fMatrix.h (linspace): Declare it. + * CMatrix.cc (linspace): New method. + * CMatrix.h (linspace): Declare it. + * fCMatrix.cc (linspace): New method. + * fCMatrix.h (linspace): Declare it. + 2009-09-17 John W. Eaton <jwe@octave.org> * oct-types.h.in: Delete.
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -3395,6 +3395,39 @@ return result; } +Matrix linspace (const ColumnVector& x1, + const ColumnVector& x2, + octave_idx_type n) + +{ + if (n < 1) n = 1; + + octave_idx_type m = x1.length (); + + if (x2.length () != m) + (*current_liboctave_error_handler) ("linspace: vectors must be of equal length"); + + NoAlias<Matrix> retval; + + retval.clear (m, n); + for (octave_idx_type i = 0; i < m; i++) + retval(i, 0) = x1(i); + + // The last column is not needed while using delta. + double *delta = &retval(0, 1); + for (octave_idx_type i = 0; i < m; i++) + delta[i] = (x2(i) - x1(i)) / (n - 1); + + for (octave_idx_type j = 1; j < n-1; j++) + for (octave_idx_type i = 0; i < m; i++) + retval(i, j) = retval(i, j-1) + delta[i]; + + for (octave_idx_type i = 0; i < m; i++) + retval(i, n-1) = x2(i); + + return retval; +} + MS_CMP_OPS (Matrix, double) MS_BOOL_OPS (Matrix, double)
--- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -357,6 +357,10 @@ extern OCTAVE_API Matrix max (const Matrix& m, double d); extern OCTAVE_API Matrix max (const Matrix& a, const Matrix& b); +extern OCTAVE_API Matrix linspace (const ColumnVector& x1, + const ColumnVector& x2, + octave_idx_type n); + MS_CMP_OP_DECLS (Matrix, double, OCTAVE_API) MS_BOOL_OP_DECLS (Matrix, double, OCTAVE_API)
--- a/liboctave/dRowVector.cc +++ b/liboctave/dRowVector.cc @@ -311,22 +311,15 @@ RowVector linspace (double x1, double x2, octave_idx_type n) { - RowVector retval; + if (n < 1) n = 1; + + NoAlias<RowVector> retval (n); - if (n > 1) - { - retval.resize (n); - double delta = (x2 - x1) / (n - 1); - retval.elem (0) = x1; - for (octave_idx_type i = 1; i < n-1; i++) - retval.elem (i) = x1 + i * delta; - retval.elem (n-1) = x2; - } - else - { - retval.resize (1); - retval.elem (0) = x2; - } + double delta = (x2 - x1) / (n - 1); + double y = retval(0) = x1; + for (octave_idx_type i = 1; i < n-1; i++) + retval(i) = y += delta; + retval(n-1) = x2; return retval; }
--- a/liboctave/fCMatrix.cc +++ b/liboctave/fCMatrix.cc @@ -4059,6 +4059,39 @@ return result; } +FloatComplexMatrix linspace (const FloatComplexColumnVector& x1, + const FloatComplexColumnVector& x2, + octave_idx_type n) + +{ + if (n < 1) n = 1; + + octave_idx_type m = x1.length (); + + if (x2.length () != m) + (*current_liboctave_error_handler) ("linspace: vectors must be of equal length"); + + NoAlias<FloatComplexMatrix> retval; + + retval.clear (m, n); + for (octave_idx_type i = 0; i < m; i++) + retval(i, 0) = x1(i); + + // The last column is not needed while using delta. + FloatComplex *delta = &retval(0, 1); + for (octave_idx_type i = 0; i < m; i++) + delta[i] = (x2(i) - x1(i)) / (n - 1.0f); + + for (octave_idx_type j = 1; j < n-1; j++) + for (octave_idx_type i = 0; i < m; i++) + retval(i, j) = retval(i, j-1) + delta[i]; + + for (octave_idx_type i = 0; i < m; i++) + retval(i, n-1) = x2(i); + + return retval; +} + MS_CMP_OPS (FloatComplexMatrix, FloatComplex) MS_BOOL_OPS (FloatComplexMatrix, FloatComplex)
--- a/liboctave/fCMatrix.h +++ b/liboctave/fCMatrix.h @@ -408,6 +408,10 @@ extern OCTAVE_API FloatComplexMatrix max (const FloatComplexMatrix& m, const FloatComplex& c); extern OCTAVE_API FloatComplexMatrix max (const FloatComplexMatrix& a, const FloatComplexMatrix& b); +extern OCTAVE_API FloatComplexMatrix linspace (const FloatComplexColumnVector& x1, + const FloatComplexColumnVector& x2, + octave_idx_type n); + MS_CMP_OP_DECLS (FloatComplexMatrix, FloatComplex, OCTAVE_API) MS_BOOL_OP_DECLS (FloatComplexMatrix, FloatComplex, OCTAVE_API)
--- a/liboctave/fCRowVector.cc +++ b/liboctave/fCRowVector.cc @@ -482,22 +482,15 @@ FloatComplexRowVector linspace (const FloatComplex& x1, const FloatComplex& x2, octave_idx_type n) { - FloatComplexRowVector retval; + if (n < 1) n = 1; + + NoAlias<FloatComplexRowVector> retval (n); - if (n > 0) - { - retval.resize (n); - FloatComplex delta = (x2 - x1) / static_cast<float> (n - 1.0); - retval.elem (0) = x1; - for (octave_idx_type i = 1; i < n-1; i++) - retval.elem (i) = x1 + static_cast<float> (1.0) * i * delta; - retval.elem (n-1) = x2; - } - else - { - retval.resize (1); - retval.elem (0) = x2; - } + FloatComplex delta = (x2 - x1) / (n - 1.0f); + FloatComplex y = retval(0) = x1; + for (octave_idx_type i = 1; i < n-1; i++) + retval(i) = y += delta; + retval(n-1) = x2; return retval; }
--- a/liboctave/fMatrix.cc +++ b/liboctave/fMatrix.cc @@ -3394,6 +3394,39 @@ return result; } +FloatMatrix linspace (const FloatColumnVector& x1, + const FloatColumnVector& x2, + octave_idx_type n) + +{ + if (n < 1) n = 1; + + octave_idx_type m = x1.length (); + + if (x2.length () != m) + (*current_liboctave_error_handler) ("linspace: vectors must be of equal length"); + + NoAlias<FloatMatrix> retval; + + retval.clear (m, n); + for (octave_idx_type i = 0; i < m; i++) + retval(i, 0) = x1(i); + + // The last column is not needed while using delta. + float *delta = &retval(0, 1); + for (octave_idx_type i = 0; i < m; i++) + delta[i] = (x2(i) - x1(i)) / (n - 1); + + for (octave_idx_type j = 1; j < n-1; j++) + for (octave_idx_type i = 0; i < m; i++) + retval(i, j) = retval(i, j-1) + delta[i]; + + for (octave_idx_type i = 0; i < m; i++) + retval(i, n-1) = x2(i); + + return retval; +} + MS_CMP_OPS (FloatMatrix, float) MS_BOOL_OPS (FloatMatrix, float)
--- a/liboctave/fMatrix.h +++ b/liboctave/fMatrix.h @@ -358,6 +358,11 @@ extern OCTAVE_API FloatMatrix max (const FloatMatrix& m, float d); extern OCTAVE_API FloatMatrix max (const FloatMatrix& a, const FloatMatrix& b); +extern OCTAVE_API FloatMatrix linspace (const FloatColumnVector& x1, + const FloatColumnVector& x2, + octave_idx_type n); + + MS_CMP_OP_DECLS (FloatMatrix, float, OCTAVE_API) MS_BOOL_OP_DECLS (FloatMatrix, float, OCTAVE_API)
--- a/liboctave/fRowVector.cc +++ b/liboctave/fRowVector.cc @@ -311,22 +311,15 @@ FloatRowVector linspace (float x1, float x2, octave_idx_type n) { - FloatRowVector retval; + if (n < 1) n = 1; + + NoAlias<FloatRowVector> retval (n); - if (n > 1) - { - retval.resize (n); - float delta = (x2 - x1) / (n - 1); - retval.elem (0) = x1; - for (octave_idx_type i = 1; i < n-1; i++) - retval.elem (i) = x1 + i * delta; - retval.elem (n-1) = x2; - } - else - { - retval.resize (1); - retval.elem (0) = x2; - } + float delta = (x2 - x1) / (n - 1); + float y = retval(0) = x1; + for (octave_idx_type i = 1; i < n-1; i++) + retval(i) = y += delta; + retval(n-1) = x2; return retval; }