Mercurial > octave-nkf

changeset 8736:53b4fdeacc2e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/CMatrix.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/CNDArray.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/ChangeLog	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/dMatrix.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/dNDArray.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/fCMatrix.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/fCNDArray.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/fMatrix.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/fNDArray.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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	Fri Feb 13 13:03:04 2009 -0500
+++ b/liboctave/mx-inlines.cc	Fri Feb 13 21:04:50 2009 +0100
@@ -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) \