Mercurial > hg > octave-terminal
changeset 8736:53b4fdeacc2e
improve reduction functions
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Fri, 13 Feb 2009 21:04:50 +0100 |
parents | afbfd7f4fd93 |
children | ae51dc447bab |
files | liboctave/CMatrix.cc liboctave/CNDArray.cc liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dNDArray.cc liboctave/fCMatrix.cc liboctave/fCNDArray.cc liboctave/fMatrix.cc liboctave/fNDArray.cc liboctave/mx-inlines.cc |
diffstat | 10 files changed, 295 insertions(+), 109 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -3264,42 +3264,31 @@ ComplexMatrix ComplexMatrix::cumprod (int dim) const { - MX_CUMULATIVE_OP (ComplexMatrix, Complex, *=); + return do_mx_cum_op<ComplexMatrix> (*this, dim, mx_inline_cumprod); } ComplexMatrix ComplexMatrix::cumsum (int dim) const { - MX_CUMULATIVE_OP (ComplexMatrix, Complex, +=); + return do_mx_cum_op<ComplexMatrix> (*this, dim, mx_inline_cumsum); } ComplexMatrix ComplexMatrix::prod (int dim) const { - MX_REDUCTION_OP (ComplexMatrix, *=, 1.0, 1.0); + return do_mx_red_op<ComplexMatrix> (*this, dim, mx_inline_prod); } ComplexMatrix ComplexMatrix::sum (int dim) const { - MX_REDUCTION_OP (ComplexMatrix, +=, 0.0, 0.0); + return do_mx_red_op<ComplexMatrix> (*this, dim, mx_inline_sum); } ComplexMatrix ComplexMatrix::sumsq (int dim) const { -#define ROW_EXPR \ - Complex d = elem (i, j); \ - retval.elem (i, 0) += d * conj (d) - -#define COL_EXPR \ - Complex d = elem (i, j); \ - retval.elem (0, j) += d * conj (d) - - MX_BASE_REDUCTION_OP (ComplexMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0); - -#undef ROW_EXPR -#undef COL_EXPR + return do_mx_red_op<Matrix> (*this, dim, mx_inline_sumsq); } Matrix ComplexMatrix::abs (void) const
--- a/liboctave/CNDArray.cc +++ b/liboctave/CNDArray.cc @@ -639,34 +639,31 @@ ComplexNDArray ComplexNDArray::cumprod (int dim) const { - MX_ND_CUMULATIVE_OP (ComplexNDArray, Complex, Complex (1, 0), *); + return do_mx_cum_op<ComplexNDArray> (*this, dim, mx_inline_cumprod); } ComplexNDArray ComplexNDArray::cumsum (int dim) const { - MX_ND_CUMULATIVE_OP (ComplexNDArray, Complex, Complex (0, 0), +); + return do_mx_cum_op<ComplexNDArray> (*this, dim, mx_inline_cumsum); } ComplexNDArray ComplexNDArray::prod (int dim) const { - MX_ND_COMPLEX_OP_REDUCTION (*= elem (iter_idx), Complex (1, 0)); + return do_mx_red_op<ComplexNDArray> (*this, dim, mx_inline_prod); +} + +ComplexNDArray +ComplexNDArray::sum (int dim) const +{ + return do_mx_red_op<ComplexNDArray> (*this, dim, mx_inline_sum); } ComplexNDArray ComplexNDArray::sumsq (int dim) const { - MX_ND_COMPLEX_OP_REDUCTION - (+= std::imag (elem (iter_idx)) - ? elem (iter_idx) * conj (elem (iter_idx)) - : std::pow (elem (iter_idx), 2), Complex (0, 0)); -} - -ComplexNDArray -ComplexNDArray::sum (int dim) const -{ - MX_ND_COMPLEX_OP_REDUCTION (+= elem (iter_idx), Complex (0, 0)); + return do_mx_red_op<NDArray> (*this, dim, mx_inline_sumsq); } ComplexNDArray
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,24 @@ +2009-02-13 Jaroslav Hajek <highegg@gmail.com> + + * mx-inlines.cc (OP_RED_SUM, OP_RED_PROD, OP_RED_SUMSQ, OP_RED_SUMSQC, + OP_RED_FCN, OP_RED_FCN2, OP_RED_FCNN, OP_CUM_FCN, OP_CUM_FCN2, + OP_CUM_FCNN): New macros. + (mx_inline_sum, mx_inline_prod, mx_inline_sumsq, mx_inline_cumsum, + mx_inline_cumprod, get_extent_triplet, do_mx_red_op, do_mx_cum_op): + New template functions. + * dMatrix.cc (Matrix::cumprod, Matrix::cumsum, Matrix::prod, + Matrix::sum, Matrix::sumsq): Use do_mx_red_op and do_mx_cum_op. + * fMatrix.cc (FloatMatrix::cumprod, FloatMatrix::cumsum, + FloatMatrix::prod, FloatMatrix::sum, FloatMatrix::sumsq): Use + do_mx_red_op and do_mx_cum_op. + * CMatrix.cc (ComplexMatrix::cumprod, ComplexMatrix::cumsum, + ComplexMatrix::prod, ComplexMatrix::sum, ComplexMatrix::sumsq): Use + do_mx_red_op and do_mx_cum_op. + * fCMatrix.cc (FloatComplexMatrix::cumprod, + FloatComplexMatrix::cumsum, FloatComplexMatrix::prod, + FloatComplexMatrix::sum, FloatComplexMatrix::sumsq): Use do_mx_red_op + and do_mx_cum_op. + 2009-02-12 Jaroslav Hajek <highegg@gmail.com> * oct-inttypes.h (if_else_type): Remove
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -2798,42 +2798,31 @@ Matrix Matrix::cumprod (int dim) const { - MX_CUMULATIVE_OP (Matrix, double, *=); + return do_mx_cum_op<Matrix> (*this, dim, mx_inline_cumprod); } Matrix Matrix::cumsum (int dim) const { - MX_CUMULATIVE_OP (Matrix, double, +=); + return do_mx_cum_op<Matrix> (*this, dim, mx_inline_cumsum); } Matrix Matrix::prod (int dim) const { - MX_REDUCTION_OP (Matrix, *=, 1.0, 1.0); + return do_mx_red_op<Matrix> (*this, dim, mx_inline_prod); } Matrix Matrix::sum (int dim) const { - MX_REDUCTION_OP (Matrix, +=, 0.0, 0.0); + return do_mx_red_op<Matrix> (*this, dim, mx_inline_sum); } Matrix Matrix::sumsq (int dim) const { -#define ROW_EXPR \ - double d = elem (i, j); \ - retval.elem (i, 0) += d * d - -#define COL_EXPR \ - double d = elem (i, j); \ - retval.elem (0, j) += d * d - - MX_BASE_REDUCTION_OP (Matrix, ROW_EXPR, COL_EXPR, 0.0, 0.0); - -#undef ROW_EXPR -#undef COL_EXPR + return do_mx_red_op<Matrix> (*this, dim, mx_inline_sumsq); } Matrix
--- a/liboctave/dNDArray.cc +++ b/liboctave/dNDArray.cc @@ -703,31 +703,31 @@ NDArray NDArray::cumprod (int dim) const { - MX_ND_CUMULATIVE_OP (NDArray, double, 1, *); + return do_mx_cum_op<NDArray> (*this, dim, mx_inline_cumprod); } NDArray NDArray::cumsum (int dim) const { - MX_ND_CUMULATIVE_OP (NDArray, double, 0, +); + return do_mx_cum_op<NDArray> (*this, dim, mx_inline_cumsum); } NDArray NDArray::prod (int dim) const { - MX_ND_REAL_OP_REDUCTION (*= elem (iter_idx), 1); + return do_mx_red_op<NDArray> (*this, dim, mx_inline_prod); +} + +NDArray +NDArray::sum (int dim) const +{ + return do_mx_red_op<NDArray> (*this, dim, mx_inline_sum); } NDArray NDArray::sumsq (int dim) const { - MX_ND_REAL_OP_REDUCTION (+= std::pow (elem (iter_idx), 2), 0); -} - -NDArray -NDArray::sum (int dim) const -{ - MX_ND_REAL_OP_REDUCTION (+= elem (iter_idx), 0); + return do_mx_red_op<NDArray> (*this, dim, mx_inline_sumsq); } NDArray
--- a/liboctave/fCMatrix.cc +++ b/liboctave/fCMatrix.cc @@ -3298,42 +3298,31 @@ FloatComplexMatrix FloatComplexMatrix::cumprod (int dim) const { - MX_CUMULATIVE_OP (FloatComplexMatrix, FloatComplex, *=); + return do_mx_cum_op<FloatComplexMatrix> (*this, dim, mx_inline_cumprod); } FloatComplexMatrix FloatComplexMatrix::cumsum (int dim) const { - MX_CUMULATIVE_OP (FloatComplexMatrix, FloatComplex, +=); + return do_mx_cum_op<FloatComplexMatrix> (*this, dim, mx_inline_cumsum); } FloatComplexMatrix FloatComplexMatrix::prod (int dim) const { - MX_REDUCTION_OP (FloatComplexMatrix, *=, 1.0, 1.0); + return do_mx_red_op<FloatComplexMatrix> (*this, dim, mx_inline_prod); } FloatComplexMatrix FloatComplexMatrix::sum (int dim) const { - MX_REDUCTION_OP (FloatComplexMatrix, +=, 0.0, 0.0); + return do_mx_red_op<FloatComplexMatrix> (*this, dim, mx_inline_sum); } FloatComplexMatrix FloatComplexMatrix::sumsq (int dim) const { -#define ROW_EXPR \ - FloatComplex d = elem (i, j); \ - retval.elem (i, 0) += d * conj (d) - -#define COL_EXPR \ - FloatComplex d = elem (i, j); \ - retval.elem (0, j) += d * conj (d) - - MX_BASE_REDUCTION_OP (FloatComplexMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0); - -#undef ROW_EXPR -#undef COL_EXPR + return do_mx_red_op<FloatMatrix> (*this, dim, mx_inline_sumsq); } FloatMatrix FloatComplexMatrix::abs (void) const
--- a/liboctave/fCNDArray.cc +++ b/liboctave/fCNDArray.cc @@ -634,33 +634,31 @@ FloatComplexNDArray FloatComplexNDArray::cumprod (int dim) const { - MX_ND_CUMULATIVE_OP (FloatComplexNDArray, FloatComplex, FloatComplex (1, 0), *); + return do_mx_cum_op<FloatComplexNDArray> (*this, dim, mx_inline_cumprod); } FloatComplexNDArray FloatComplexNDArray::cumsum (int dim) const { - MX_ND_CUMULATIVE_OP (FloatComplexNDArray, FloatComplex, FloatComplex (0, 0), +); + return do_mx_cum_op<FloatComplexNDArray> (*this, dim, mx_inline_cumsum); } FloatComplexNDArray FloatComplexNDArray::prod (int dim) const { - MX_ND_REDUCTION (retval(result_idx) *= elem (iter_idx), FloatComplex (1, 0), FloatComplexNDArray); + return do_mx_red_op<FloatComplexNDArray> (*this, dim, mx_inline_prod); +} + +FloatComplexNDArray +FloatComplexNDArray::sum (int dim) const +{ + return do_mx_red_op<FloatComplexNDArray> (*this, dim, mx_inline_sum); } FloatComplexNDArray FloatComplexNDArray::sumsq (int dim) const { - MX_ND_REDUCTION (retval(result_idx) += std::imag (elem (iter_idx)) - ? elem (iter_idx) * conj (elem (iter_idx)) - : std::pow (elem (iter_idx), 2), FloatComplex (0, 0), FloatComplexNDArray); -} - -FloatComplexNDArray -FloatComplexNDArray::sum (int dim) const -{ - MX_ND_REDUCTION (retval(result_idx) += elem (iter_idx), FloatComplex (0, 0), FloatComplexNDArray); + return do_mx_red_op<FloatNDArray> (*this, dim, mx_inline_sumsq); } FloatComplexNDArray
--- a/liboctave/fMatrix.cc +++ b/liboctave/fMatrix.cc @@ -2797,42 +2797,31 @@ FloatMatrix FloatMatrix::cumprod (int dim) const { - MX_CUMULATIVE_OP (FloatMatrix, float, *=); + return do_mx_cum_op<FloatMatrix> (*this, dim, mx_inline_cumprod); } FloatMatrix FloatMatrix::cumsum (int dim) const { - MX_CUMULATIVE_OP (FloatMatrix, float, +=); + return do_mx_cum_op<FloatMatrix> (*this, dim, mx_inline_cumsum); } FloatMatrix FloatMatrix::prod (int dim) const { - MX_REDUCTION_OP (FloatMatrix, *=, 1.0, 1.0); + return do_mx_red_op<FloatMatrix> (*this, dim, mx_inline_prod); } FloatMatrix FloatMatrix::sum (int dim) const { - MX_REDUCTION_OP (FloatMatrix, +=, 0.0, 0.0); + return do_mx_red_op<FloatMatrix> (*this, dim, mx_inline_sum); } FloatMatrix FloatMatrix::sumsq (int dim) const { -#define ROW_EXPR \ - float d = elem (i, j); \ - retval.elem (i, 0) += d * d - -#define COL_EXPR \ - float d = elem (i, j); \ - retval.elem (0, j) += d * d - - MX_BASE_REDUCTION_OP (FloatMatrix, ROW_EXPR, COL_EXPR, 0.0, 0.0); - -#undef ROW_EXPR -#undef COL_EXPR + return do_mx_red_op<FloatMatrix> (*this, dim, mx_inline_sumsq); } FloatMatrix
--- a/liboctave/fNDArray.cc +++ b/liboctave/fNDArray.cc @@ -658,31 +658,31 @@ FloatNDArray FloatNDArray::cumprod (int dim) const { - MX_ND_CUMULATIVE_OP (FloatNDArray, float, 1, *); + return do_mx_cum_op<FloatNDArray> (*this, dim, mx_inline_cumprod); } FloatNDArray FloatNDArray::cumsum (int dim) const { - MX_ND_CUMULATIVE_OP (FloatNDArray, float, 0, +); + return do_mx_cum_op<FloatNDArray> (*this, dim, mx_inline_cumsum); } FloatNDArray FloatNDArray::prod (int dim) const { - MX_ND_REDUCTION (retval(result_idx) *= elem (iter_idx), 1, FloatNDArray); + return do_mx_red_op<FloatNDArray> (*this, dim, mx_inline_prod); +} + +FloatNDArray +FloatNDArray::sum (int dim) const +{ + return do_mx_red_op<FloatNDArray> (*this, dim, mx_inline_sum); } FloatNDArray FloatNDArray::sumsq (int dim) const { - MX_ND_REDUCTION (retval(result_idx) += std::pow (elem (iter_idx), 2), 0, FloatNDArray); -} - -FloatNDArray -FloatNDArray::sum (int dim) const -{ - MX_ND_REDUCTION (retval(result_idx) += elem (iter_idx), 0, FloatNDArray); + return do_mx_red_op<FloatNDArray> (*this, dim, mx_inline_sumsq); } FloatNDArray
--- a/liboctave/mx-inlines.cc +++ b/liboctave/mx-inlines.cc @@ -282,6 +282,220 @@ OP_DUP_FCN (imag, mx_inline_imag_dup, float, FloatComplex) OP_DUP_FCN (conj, mx_inline_conj_dup, FloatComplex, FloatComplex) +// NOTE: std::norm is NOT equivalent +template <class T> +T cabsq (const std::complex<T>& c) +{ return c.real () * c.real () + c.imag () * c.imag (); } + +#define OP_RED_SUM(ac, el) ac += el +#define OP_RED_PROD(ac, el) ac *= el +#define OP_RED_SUMSQ(ac, el) ac += el*el +#define OP_RED_SUMSQC(ac, el) ac += cabsq (el) + +#define OP_RED_FCN(F, TSRC, OP, ZERO) \ +template <class T> \ +inline T \ +F (const TSRC* v, octave_idx_type n) \ +{ \ + T ac = ZERO; \ + for (octave_idx_type i = 0; i < n; i++) \ + OP(ac, v[i]); \ + return ac; \ +} + +OP_RED_FCN (mx_inline_sum, T, OP_RED_SUM, 0) +OP_RED_FCN (mx_inline_prod, T, OP_RED_PROD, 1) +OP_RED_FCN (mx_inline_sumsq, T, OP_RED_SUMSQ, 0) +OP_RED_FCN (mx_inline_sumsq, std::complex<T>, OP_RED_SUMSQC, 0) + +#define OP_RED_FCN2(F, TSRC, OP, ZERO) \ +template <class T> \ +inline void \ +F (const TSRC* v, T *r, octave_idx_type m, octave_idx_type n) \ +{ \ + for (octave_idx_type i = 0; i < m; i++) \ + r[i] = ZERO; \ + for (octave_idx_type j = 0; j < n; j++) \ + { \ + for (octave_idx_type i = 0; i < m; i++) \ + OP(r[i], v[i]); \ + v += m; \ + } \ +} + +OP_RED_FCN2 (mx_inline_sum, T, OP_RED_SUM, 0) +OP_RED_FCN2 (mx_inline_prod, T, OP_RED_PROD, 1) +OP_RED_FCN2 (mx_inline_sumsq, T, OP_RED_SUMSQ, 0) +OP_RED_FCN2 (mx_inline_sumsq, std::complex<T>, OP_RED_SUMSQC, 0) + +#define OP_RED_FCNN(F, TSRC) \ +template <class T> \ +inline void \ +F (const TSRC *v, T *r, octave_idx_type l, \ + octave_idx_type n, octave_idx_type u) \ +{ \ + if (l == 1) \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + r[i] = F (v, n); \ + v += n; \ + } \ + } \ + else \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, l, n); \ + v += l*n; \ + r += l; \ + } \ + } \ +} + +OP_RED_FCNN (mx_inline_sum, T) +OP_RED_FCNN (mx_inline_prod, T) +OP_RED_FCNN (mx_inline_sumsq, T) +OP_RED_FCNN (mx_inline_sumsq, std::complex<T>) + +#define OP_CUM_FCN(F, OP) \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type n) \ +{ \ + if (n) \ + { \ + T t = r[0] = v[0]; \ + for (octave_idx_type i = 1; i < n; i++) \ + r[i] = t = t OP v[i]; \ + } \ +} + +OP_CUM_FCN (mx_inline_cumsum, +) +OP_CUM_FCN (mx_inline_cumprod, *) + +#define OP_CUM_FCN2(F, OP) \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \ +{ \ + if (n) \ + { \ + for (octave_idx_type i = 0; i < m; i++) \ + r[i] = v[i]; \ + const T *r0 = r; \ + for (octave_idx_type j = 1; j < n; j++) \ + { \ + r += m; v += m; \ + for (octave_idx_type i = 0; i < m; i++) \ + r[i] = v[i] OP r0[i]; \ + r0 += m; \ + } \ + } \ +} + +OP_CUM_FCN2 (mx_inline_cumsum, +) +OP_CUM_FCN2 (mx_inline_cumprod, *) + +#define OP_CUM_FCNN(F) \ +template <class T> \ +inline void \ +F (const T *v, T *r, octave_idx_type l, \ + octave_idx_type n, octave_idx_type u) \ +{ \ + if (l == 1) \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, n); \ + v += n; r += n; \ + } \ + } \ + else \ + { \ + for (octave_idx_type i = 0; i < u; i++) \ + { \ + F (v, r, l, n); \ + v += l*n; \ + r += l*n; \ + } \ + } \ +} + +OP_CUM_FCNN (mx_inline_cumsum) +OP_CUM_FCNN (mx_inline_cumprod) + +// Assistant function + +inline void +get_extent_triplet (const dim_vector& dims, int& dim, + octave_idx_type& l, octave_idx_type& n, + octave_idx_type& u) +{ + octave_idx_type ndims = dims.length (); + if (dim >= ndims) + { + l = dims.numel (); + n = 1; + u = 1; + } + else + { + if (dim < 0) + { + // find first non-singleton dim + for (dim = 0; dims(dim) == 1 && dim < ndims - 1; dim++) ; + } + // calculate extent triplet. + l = 1, n = dims(dim), u = 1; + for (octave_idx_type i = 0; i < dim; i++) + l *= dims (i); + for (octave_idx_type i = dim + 1; i < ndims; i++) + u *= dims (i); + } +} + +// Appliers. +// FIXME: is this the best design? C++ gives a lot of options here... +// maybe it can be done without an explicit parameter? + +template <class ArrayType, class T> +inline ArrayType +do_mx_red_op (const Array<T>& src, int dim, + void (*mx_red_op) (const T *, typename ArrayType::element_type *, + octave_idx_type, octave_idx_type, octave_idx_type)) +{ + octave_idx_type l, n, u; + dim_vector dims = src.dims (); + get_extent_triplet (dims, dim, l, n, u); + + // Reduction operation reduces the array size. + if (dim < dims.length ()) dims(dim) = 1; + dims.chop_trailing_singletons (); + + ArrayType ret (dims); + mx_red_op (src.data (), ret.fortran_vec (), l, n, u); + + return ret; +} + +template <class ArrayType, class T> +inline ArrayType +do_mx_cum_op (const Array<T>& src, int dim, + void (*mx_cum_op) (const T *, typename ArrayType::element_type *, + octave_idx_type, octave_idx_type, octave_idx_type)) +{ + octave_idx_type l, n, u; + dim_vector dims = src.dims (); + get_extent_triplet (dims, dim, l, n, u); + + // Cumulative operation doesn't reduce the array size. + ArrayType ret (dims); + mx_cum_op (src.data (), ret.fortran_vec (), l, n, u); + + return ret; +} + // Avoid some code duplication. Maybe we should use templates. #define MX_CUMULATIVE_OP(RET_TYPE, ELT_TYPE, OP) \