--- a/liboctave/mx-inlines.cc
+++ b/liboctave/mx-inlines.cc
@@ -2,6 +2,8 @@
 
 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001, 2002,
               2003, 2004, 2005, 2006, 2007 John W. Eaton
+Copyright (C) 2009 Jaroslav Hajek
+Copyright (C) 2009 VZLU Prague
 
 This file is part of Octave.
 
@@ -656,6 +658,186 @@
 OP_MINMAX_FCNN (mx_inline_min)
 OP_MINMAX_FCNN (mx_inline_max)
 
+#define OP_CUMMINMAX_FCN(F, OP) \
+template <class T> \
+void F (const T *v, T *r, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  T tmp = v[0]; \
+  octave_idx_type i = 1, j = 0; \
+  if (xisnan (tmp)) \
+    { \
+      for (; i < n && xisnan (v[i]); i++) ; \
+      for (; j < i; j++) r[j] = tmp; \
+      if (i < n) tmp = v[i]; \
+    } \
+  for (; i < n; i++) \
+    if (v[i] OP tmp) \
+      { \
+        for (; j < i; j++) r[j] = tmp; \
+        tmp = v[i]; \
+      } \
+  for (; j < i; j++) r[j] = tmp; \
+} \
+template <class T> \
+void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  T tmp = v[0]; octave_idx_type tmpi = 0; \
+  octave_idx_type i = 1, j = 0; \
+  if (xisnan (tmp)) \
+    { \
+      for (; i < n && xisnan (v[i]); i++) ; \
+      for (; j < i; j++) { r[j] = tmp; ri[j] = tmpi; } \
+      if (i < n) { tmp = v[i]; tmpi = i; } \
+    } \
+  for (; i < n; i++) \
+    if (v[i] OP tmp) \
+      { \
+        for (; j < i; j++) { r[j] = tmp; ri[j] = tmpi; } \
+        tmp = v[i]; tmpi = i; \
+      } \
+  for (; j < i; j++) { r[j] = tmp; ri[j] = tmpi; } \
+}
+
+OP_CUMMINMAX_FCN (mx_inline_cummin, <)
+OP_CUMMINMAX_FCN (mx_inline_cummax, >)
+
+// Row reductions will be slightly complicated.  We will proceed with checks
+// for NaNs until we detect that no row will yield a NaN, in which case we
+// proceed to a faster code.
+
+#define OP_CUMMINMAX_FCN2(F, OP) \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type m, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  bool nan = false; \
+  const T *r0; \
+  octave_idx_type j = 0; \
+  for (octave_idx_type i = 0; i < m; i++) \
+    {  \
+      r[i] = v[i]; \
+      if (xisnan (v[i])) nan = true;  \
+    } \
+  j++; v += m; r0 = r; r += m; \
+  while (nan && j < n) \
+    { \
+      nan = false; \
+      for (octave_idx_type i = 0; i < m; i++) \
+        {  \
+          if (xisnan (v[i])) \
+            { r[i] = r0[i]; nan = true; } \
+          else if (xisnan (r[i]) || v[i] OP r[i]) \
+            r[i] = v[i]; \
+        } \
+      j++; v += m; r0 = r; r += m; \
+    } \
+  while (j < n) \
+    { \
+      for (octave_idx_type i = 0; i < m; i++) \
+        if (v[i] OP r[i]) \
+          r[i] = v[i]; \
+        else \
+          r[i] = r0[i]; \
+      j++; v += m; r0 = r; r += m; \
+    } \
+} \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type *ri, \
+   octave_idx_type m, octave_idx_type n) \
+{ \
+  if (! n) return; \
+  bool nan = false; \
+  const T *r0; const octave_idx_type *r0i; \
+  octave_idx_type j = 0; \
+  for (octave_idx_type i = 0; i < m; i++) \
+    {  \
+      r[i] = v[i]; ri[i] = 0; \
+      if (xisnan (v[i])) nan = true;  \
+    } \
+  j++; v += m; r0 = r; r += m; r0i = ri; ri += m;  \
+  while (nan && j < n) \
+    { \
+      nan = false; \
+      for (octave_idx_type i = 0; i < m; i++) \
+        {  \
+          if (xisnan (v[i])) \
+            { r[i] = r0[i]; ri[i] = r0i[i]; nan = true; } \
+          else if (xisnan (r[i]) || v[i] OP r[i]) \
+            { r[i] = v[i]; ri[i] = j; }\
+        } \
+      j++; v += m; r0 = r; r += m; r0i = ri; ri += m;  \
+    } \
+  while (j < n) \
+    { \
+      for (octave_idx_type i = 0; i < m; i++) \
+        if (v[i] OP r[i]) \
+          { r[i] = v[i]; ri[i] = j; } \
+        else \
+          { r[i] = r0[i]; ri[i] = r0i[i]; } \
+      j++; v += m; r0 = r; r += m; r0i = ri; ri += m;  \
+    } \
+}
+
+OP_CUMMINMAX_FCN2 (mx_inline_cummin, <)
+OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
+
+#define OP_CUMMINMAX_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 (! n) return; \
+  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; \
+        } \
+    } \
+} \
+template <class T> \
+inline void \
+F (const T *v, T *r, octave_idx_type *ri, \
+   octave_idx_type l, octave_idx_type n, octave_idx_type u) \
+{ \
+  if (! n) return; \
+  if (l == 1) \
+    { \
+      for (octave_idx_type i = 0; i < u; i++) \
+        { \
+          F (v, r, ri, n); \
+          v += n; r += n; ri += n; \
+        } \
+    } \
+  else \
+    { \
+      for (octave_idx_type i = 0; i < u; i++) \
+        { \
+          F (v, r, ri, l, n); \
+          v += l*n; \
+          r += l*n; ri += l*n; \
+        } \
+    } \
+}
+
+OP_CUMMINMAX_FCNN (mx_inline_cummin)
+OP_CUMMINMAX_FCNN (mx_inline_cummax)
+
 // Assistant function
 
 inline void
@@ -777,6 +959,44 @@
   return ret;
 }
 
+template <class ArrayType>
+inline ArrayType
+do_mx_cumminmax_op (const ArrayType& src, int dim,
+                    void (*mx_cumminmax_op) (const typename ArrayType::element_type *, 
+                                             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);
+
+  ArrayType ret (dims);
+  mx_cumminmax_op (src.data (), ret.fortran_vec (), l, n, u);
+
+  return ret;
+}
+
+template <class ArrayType>
+inline ArrayType
+do_mx_cumminmax_op (const ArrayType& src, Array<octave_idx_type>& idx, int dim,
+                    void (*mx_cumminmax_op) (const typename ArrayType::element_type *, 
+                                             typename ArrayType::element_type *,
+                                             octave_idx_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);
+
+  ArrayType ret (dims);
+  if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
+
+  mx_cumminmax_op (src.data (), ret.fortran_vec (), idx.fortran_vec (),
+                   l, n, u);
+
+  return ret;
+}
+
 #endif
 
 /*