# HG changeset patch
# User Stefan Mirea <mireastefangabriel@gmail.com>
# Date 1488972095 -7200
#      Wed Mar 08 13:21:35 2017 +0200
# Node ID 09c845ae95651be784c9953f9e1e6dcbe050adf2
# Parent  4890b1c4a6bd53f409f994d327054ddd61a5e548
Add nanflag argument to min, max, cummin, cummax (bug #50007)

* max.cc (update_minmax_bin_op_inan_full,
update_minmax_bin_op_inan_sparse): New static helper functions.
(do_minmax_body, do_cumminmax_body): Process "omitnan" and
"includenan" options and handle includenan case.
(do_cumminmax_body): New calls to do_cumminmax_red_op<SparseMatrix>
and do_cumminmax_red_op<SparseComplexMatrix>.
(Fmin, Fmax, Fcummin, Fcummax): Update docstrings with nanflag usage.
New tests for the nanflag argument.
Fix two incorrect tests. Change xtests to tests as complex values with
the same magnitude are now compared by phase angle.

* op-s-m.cc, op-s-s.cc (install_s_m_ops, install_s_s_ops): Install
assignment conversion functions.

* CMatrix.cc, fCMatrix.cc (min, max): Handle matrix - NaN scalar
comparisons correctly. Remove columns_are_real_only to keep complex
ordering.

* dMatrix.cc, fMatrix.cc (min, max): Share code between matrix-scalar
and scalar-matrix overloads as with ComplexMatrix and
FloatComplexMatrix whereas checking if d is NaN ensures safety.

* CSparse.cc (SparseComplexMatrix::max, SparseComplexMatrix::min):
Compare numbers using operators instead of std::abs to take phase
angles into account.
(SparseComplexMatrix::cummax, SparseComplexMatrix::cummin): New
methods and new test for them.
(min, max): Handle NaN scalar - matrix, null scalar - matrix and
matrix - NaN scalar comparisons correctly.

* dSparse.cc (SparseMatrix::max, SparseMatrix::min): Omit NaNs which
appear on the first position on a row with all elements nonzero.
(SparseMatrix::cummax, SparseMatrix::cummin): New methods and new test
for them.
(min, max): Handle matrix - NaN scalar comparisons correctly.

* CSparse.h, dSparse.h (SparseComplexMatrix::cummax,
SparseComplexMatrix::cummin, SparseMatrix::cummax,
SparseMatrix::cummin): New methods.

* lo-mappers.h (min, max): Omit NaN values for complex types. Use
operators instead of abs to take phase angles into account.

* mx-inlines.cc (mx_inline_xmin<double>, mx_inline_xmin<float>,
mx_inline_xmax<double>, mx_inline_xmax<float>): Remove scalar-vector
specializations to preserve NaN scalar value.
(mx_inline_min, mx_inline_max, mx_inline_cummin, mx_inline_cummax):
Proceed to the faster code in OP_MINMAX_FCN2 and OP_CUMMINMAX_FCN2 as
soon as r contains no NaNs.

diff --git a/libinterp/corefcn/max.cc b/libinterp/corefcn/max.cc
--- a/libinterp/corefcn/max.cc
+++ b/libinterp/corefcn/max.cc
@@ -40,6 +40,7 @@ along with Octave; see the file COPYING.
 #include "ov-cx-mat.h"
 #include "ov-re-sparse.h"
 #include "ov-cx-sparse.h"
+#include "parse.h"
 
 template <typename ArrayType>
 static octave_value_list
@@ -232,19 +233,128 @@ do_minmax_bin_op<charNDArray> (const oct
   return retval;
 }
 
+template <typename ArrayType>
+static void
+update_minmax_bin_op_inan_full (octave_value &result,
+                                const octave_value& argx,
+                                const octave_value& argy)
+{
+  typedef typename ArrayType::element_type ScalarType;
+
+  ArrayType x = octave_value_extract<ArrayType> (argx);
+  ArrayType y = octave_value_extract<ArrayType> (argy);
+
+  auto op_vv = [](size_t n, ScalarType *r, const ScalarType *x,
+                  const ScalarType *y)
+    {
+      for (size_t i = 0; i < n; i++)
+        r[i] = octave::math::isnan (x[i]) ? x[i] : y[i];
+    };
+
+  auto op_sv = [](size_t n, ScalarType *r, ScalarType x, const ScalarType *y)
+    {
+      if (octave::math::isnan (x))
+        std::fill_n (r, n, x);
+      else
+        std::memcpy (r, y, n * sizeof (ScalarType));
+    };
+
+  auto op_vs = [](size_t n, ScalarType *r, const ScalarType *x, ScalarType y)
+    {
+      if (octave::math::isnan (y))
+        std::fill_n (r, n, y);
+      else
+        std::memcpy (r, x, n * sizeof (ScalarType));
+    };
+
+  ArrayType r = do_bsxfun_op<ScalarType, ScalarType, ScalarType> (x, y,
+                                                                  op_vv, op_sv,
+                                                                  op_vs);
+
+  octave_value_list args;
+  args(0) = argx.isnan ();
+  args(1) = argy.isnan ();
+  octave_value_list res = octave::feval ("xor", args, 1);
+  idx_vector update_idx = res(0).bool_array_value ();
+
+  std::list<octave_value_list> update_idx_l = {ovl (update_idx)};
+  result = result.subsasgn ("(", update_idx_l, r.index (update_idx));
+}
+
+template <typename SparseType>
+static void
+update_minmax_bin_op_inan_sparse (octave_value &result,
+                                  const octave_value& argx,
+                                  const octave_value& argy)
+{
+  typedef typename SparseType::element_type ScalarType;
+
+  SparseType x = octave_value_extract<SparseType> (argx);
+  SparseType y = octave_value_extract<SparseType> (argy);
+
+  if (argx.is_scalar_type ())
+    {
+      ScalarType s = octave_value_extract<ScalarType> (argx);
+      if (octave::math::isnan (s))
+        {
+          result = SparseType (y.rows (), y.cols (), s);
+          return;
+        }
+    }
+
+  idx_vector x_nans_idx = argx.isnan ().sparse_bool_matrix_value ();
+  idx_vector y_nans_idx = argy.isnan ().sparse_bool_matrix_value ();
+
+  octave_value nan_vals;
+  std::list<octave_value_list> r_nans_idx;
+
+  if (argy.is_scalar_type ())
+    {
+      ScalarType s = octave_value_extract<ScalarType> (argy);
+      if (octave::math::isnan (s))
+        result = SparseType (x.rows (), x.cols (), s);
+    }
+  else
+    {
+      nan_vals = y.index (y_nans_idx);
+      r_nans_idx = {ovl (y_nans_idx)};
+      result = result.subsasgn ("(", r_nans_idx, nan_vals);
+    }
+
+  nan_vals = x.index (x_nans_idx);
+  r_nans_idx = {ovl (x_nans_idx)};
+  result = result.subsasgn ("(", r_nans_idx, nan_vals);
+}
+
 static octave_value_list
 do_minmax_body (const octave_value_list& args,
                 int nargout, bool ismin)
 {
   int nargin = args.length ();
 
+  bool includenan = false;
+
+  const char *func = ismin ? "min" : "max";
+
+  if (nargin >= 3 && args(nargin - 1).is_string ())
+    {
+      std::string str = args(nargin - 1).string_value ();
+
+      if (str == "includenan")
+        includenan = true;
+      else if (str != "omitnan")
+        error ("%s: unrecognized NaN flag '%s'", func, str.c_str ());
+
+      nargin--;
+      if (nargin == 2 && args(1).is_empty ())
+        nargin--;
+    }
+
   if (nargin < 1 || nargin > 3)
     print_usage ();
 
   octave_value_list retval (nargout > 1 ? 2 : 1);
 
-  const char *func = ismin ? "min" : "max";
-
   if (nargin == 3 || nargin == 1)
     {
       octave_value arg = args(0);
@@ -344,6 +454,44 @@ do_minmax_body (const octave_value_list&
         default:
           err_wrong_type_arg (func, arg);
         }
+
+      if (includenan && btyp_isfloat (arg.builtin_type ())
+          && dim < arg.ndims ())
+        {
+          Array<octave_idx_type> idx;
+          idx_vector nans_out_idx;
+
+          if (! arg.is_sparse_type ())
+            {
+              int8NDArray nans_in = arg.isnan ().int8_array_value ();
+              int8NDArray nans_out = nans_in.max (idx, dim);
+              nans_out_idx = octave_value (nans_out).bool_array_value ();
+            }
+          else
+            {
+              SparseMatrix nans_in = arg.isnan ().sparse_matrix_value ();
+              SparseMatrix nans_out = nans_in.max (idx, dim);
+              nans_out_idx =
+                octave_value (nans_out).sparse_bool_matrix_value ();
+            }
+
+          idx_vector idx_nan_only = idx.index (nans_out_idx);
+
+          dim_vector dv = idx.dims ().redim (arg.ndims ());
+          Array<idx_vector> idxa = ind2sub (dv, nans_out_idx);
+
+          if (dim < 0)
+            dim = arg.dims ().first_non_singleton ();
+          idxa(dim) = idx_nan_only;
+
+          idx_vector first_nan_idx = sub2ind (arg.dims (), idxa);
+          octave_value first_nans = arg.do_index_op (ovl (first_nan_idx));
+
+          std::list<octave_value_list> nans_out_idx_l = {ovl (nans_out_idx)};
+          retval(0) = retval(0).subsasgn ("(", nans_out_idx_l, first_nans);
+          if (nargout > 1)
+            retval(1) = retval(1).subsasgn ("(", nans_out_idx_l, idx_nan_only);
+        }
     }
   else
     {
@@ -427,6 +575,51 @@ do_minmax_body (const octave_value_list&
       if (xtyp == btyp_bool && ytyp == btyp_bool)
         retval(0) = retval(0).bool_array_value ();
 
+      if (includenan)
+        {
+          switch (rtyp)
+            {
+            case btyp_double:
+              {
+                if ((argx.is_sparse_type ()
+                     && (argy.is_sparse_type () || argy.is_scalar_type ()))
+                    || (argy.is_sparse_type () && argx.is_scalar_type ()))
+                  update_minmax_bin_op_inan_sparse<SparseMatrix> (retval(0),
+                                                                  argx, argy);
+                else
+                  update_minmax_bin_op_inan_full<NDArray> (retval(0),
+                                                           argx, argy);
+              }
+              break;
+
+            case btyp_complex:
+              {
+                if ((argx.is_sparse_type ()
+                     && (argy.is_sparse_type () || argy.is_scalar_type ()))
+                    || (argy.is_sparse_type () && argx.is_scalar_type ()))
+                  update_minmax_bin_op_inan_sparse<SparseComplexMatrix> (
+                    retval(0), argx, argy);
+                else
+                  update_minmax_bin_op_inan_full<ComplexNDArray> (retval(0),
+                                                                  argx, argy);
+              }
+              break;
+
+            case btyp_float:
+              update_minmax_bin_op_inan_full<FloatNDArray> (retval(0),
+                                                            argx, argy);
+              break;
+
+            case btyp_float_complex:
+              update_minmax_bin_op_inan_full<FloatComplexNDArray> (retval(0),
+                                                                   argx, argy);
+              break;
+
+            default:
+              // NaN values not supported.
+              break;
+            }
+        }
     }
 
   return retval;
@@ -438,6 +631,8 @@ DEFUN (min, args, nargout,
 @deftypefnx {} {} min (@var{x}, [], @var{dim})
 @deftypefnx {} {[@var{w}, @var{iw}] =} min (@var{x})
 @deftypefnx {} {} min (@var{x}, @var{y})
+@deftypefnx {} {} min (@dots{}, "omitnan")
+@deftypefnx {} {} min (@dots{}, "includenan")
 Find minimum values in the array @var{x}.
 
 For a vector argument, return the minimum value.  For a matrix argument,
@@ -451,6 +646,9 @@ set to the empty matrix.
 
 For two matrices (or a matrix and a scalar), return the pairwise minimum.
 
+The optional @qcode{"nanflag"} argument specifies if NaN values found in the
+input array(s) are to be ignored or propagated to the result.
+
 Thus,
 
 @example
@@ -538,9 +736,7 @@ the first index of the minimum value(s).
 ## Sparse double values
 %!assert (min (sparse ([1, 4, 2, 3])), sparse (1))
 %!assert (min (sparse ([1; -10; 5; -2])), sparse(-10))
-## FIXME: sparse doesn't order complex values by phase angle
-%!xtest
-%! assert (min (sparse ([4, 2i 4.999; -2, 2, 3+4i])), sparse ([-2, 2, 4.999]));
+%!assert (min (sparse ([4, 2i 4.999; -2, 2, 3+4i])), sparse ([-2, 2, 4.999]));
 
 ## Test dimension argument
 %!test
@@ -553,6 +749,26 @@ the first index of the minimum value(s).
 %! assert (ndims (i), 2);
 %! assert (i, [1, 1; 1, 1]);
 
+## Test nanflag argument
+%!assert (min ([3, NaN, 1, 2], [], "omitnan"), 1)
+%!test
+%! [y, i] = min (single ([3, NA, 1, NaN, 2]), [], "includenan");
+%! assert (y, single (NA));
+%! assert (i, 2);
+%!test
+%! x = sparse ([NA, 0, NaN; NaN, -2, 0]);
+%! assert (min (x, [], "includenan"), sparse ([NA, -2, NaN]));
+%!test
+%! x = [1, NA+4i; NaN+6i, 3];
+%! y = [2; NaN+5i];
+%! assert (min (x, y, "omitnan"), [1, 2; NaN+6i, 3]);
+%! assert (min (x, y, "includenan"), [1, NA+4i; NaN+6i, NaN+5i]);
+%!test
+%! x = reshape ([1, 2, 3, NaN, NaN, NA, 4, NaN], [2,2,2]);
+%! [y, i] = min (x, [], 2, "includenan");
+%! assert (y, reshape ([1, NaN, NaN, NA], [2,1,2]));
+%! assert (i, reshape ([1, 2, 1, 1], [2,1,2]));
+
 ## Test 2-output forms for various arg types
 ## Special routines for char arrays
 %!test
@@ -584,12 +800,9 @@ the first index of the minimum value(s).
 %! assert (min (x, 3), [1 2 3 3]);
 %! assert (min (2, x), [1 2 2 2]);
 %! assert (min (x, 2.1i), [1 2 2.1i 2.1i]);
-## FIXME: Ordering of complex results with equal magnitude is not by phase
-##        angle in the 2-input form.  Instead, it is in the order in which it
-##        appears in the argument list.
-%!xtest
+%!test
 %! x = [1, 2, 3, 4];  y = fliplr (x);
-%! assert (min (x, 2i), [2i 2i 3 4]);
+%! assert (min (x, 2i), [1 2 2i 2i]);
 ## Special routines for char arrays
 %!assert (min ("abc", "b"), [97 98 98])
 %!assert (min ("b", "cba"), [98 98 97])
@@ -648,6 +861,7 @@ the first index of the minimum value(s).
 %!error min ()
 %!error min (1, 2, 3, 4)
 %!error <DIM must be a valid dimension> min ([1 2; 3 4], [], -3)
+%!error <unrecognized NaN flag 'foobar'> min (1, [], "foobar")
 %!warning <second argument is ignored> min ([1 2 3 4], 2, 2);
 %!error <wrong type argument 'cell'> min ({1 2 3 4})
 %!error <cannot compute min \(cell, scalar\)> min ({1, 2, 3}, 2)
@@ -659,6 +873,8 @@ DEFUN (max, args, nargout,
 @deftypefnx {} {} max (@var{x}, [], @var{dim})
 @deftypefnx {} {[@var{w}, @var{iw}] =} max (@var{x})
 @deftypefnx {} {} max (@var{x}, @var{y})
+@deftypefnx {} {} max (@dots{}, "omitnan")
+@deftypefnx {} {} max (@dots{}, "includenan")
 Find maximum values in the array @var{x}.
 
 For a vector argument, return the maximum value.  For a matrix argument,
@@ -672,6 +888,9 @@ set to the empty matrix.
 
 For two matrices (or a matrix and a scalar), return the pairwise maximum.
 
+The optional @qcode{"nanflag"} argument specifies if NaN values found in the
+input array(s) are to be ignored or propagated to the result.
+
 Thus,
 
 @example
@@ -764,13 +983,33 @@ the first index of the maximum value(s).
 ## Test dimension argument
 %!test
 %! x = reshape (1:8, [2,2,2]);
-%! assert (min (x, [], 1), reshape ([1, 3, 5, 7], [1,2,2]));
-%! assert (min (x, [], 2), reshape ([1, 2, 5, 6], [2,1,2]));
-%! [y, i] = min (x, [], 3);
+%! assert (max (x, [], 1), reshape ([2, 4, 6, 8], [1,2,2]));
+%! assert (max (x, [], 2), reshape ([3, 4, 7, 8], [2,1,2]));
+%! [y, i] = max (x, [], 3);
 %! assert (ndims (y), 2);
-%! assert (y, [1, 3; 2, 4]);
+%! assert (y, [5, 7; 6, 8]);
 %! assert (ndims (i), 2);
-%! assert (i, [1, 1; 1, 1]);
+%! assert (i, [2, 2; 2, 2]);
+
+## Test nanflag argument
+%!assert (max ([3, NaN, 1, 2], [], "omitnan"), 3)
+%!test
+%! [y, i] = max (single ([3, NA, 1, NaN, 2]), [], "includenan");
+%! assert (y, single (NA));
+%! assert (i, 2);
+%!test
+%! x = sparse ([NA, -2, 0; NaN, 0, NaN]);
+%! assert (max (x, [], "includenan"), sparse ([NA, 0, NaN]));
+%!test
+%! x = [1, NA+4i; NaN+6i, 3];
+%! y = [2; NaN+5i];
+%! assert (max (x, y, "omitnan"), [2, 2; NaN+6i, 3]);
+%! assert (max (x, y, "includenan"), [2, NA+4i; NaN+6i, NaN+5i]);
+%!test
+%! x = reshape ([1, 2, 3, NaN, NaN, NA, 4, NaN], [2,2,2]);
+%! [y, i] = max (x, [], 2, "includenan");
+%! assert (y, reshape ([3, NaN, NaN, NA], [2,1,2]));
+%! assert (i, reshape ([2, 2, 1, 1], [2,1,2]));
 
 ## Test 2-output forms for various arg types
 ## Special routines for char arrays
@@ -803,10 +1042,7 @@ the first index of the maximum value(s).
 %! assert (max (x, 3), [3 3 3 4]);
 %! assert (max (2, x), [2 2 3 4]);
 %! assert (max (x, 2.1i), [2.1i 2.1i 3 4]);
-## FIXME: Ordering of complex results with equal magnitude is not by phase
-##        angle in the 2-input form.  Instead, it is in the order in which it
-##        appears in the argument list.
-%!xtest
+%!test
 %! x = [1, 2, 3, 4];  y = fliplr (x);
 %! assert (max (x, 2i), [2i 2i 3 4]);
 ## Special routines for char arrays
@@ -872,6 +1108,7 @@ the first index of the maximum value(s).
 %!error max ()
 %!error max (1, 2, 3, 4)
 %!error <DIM must be a valid dimension> max ([1 2; 3 4], [], -3)
+%!error <unrecognized NaN flag 'foobar'> max (1, [], "foobar")
 %!warning <second argument is ignored> max ([1 2 3 4], 2, 2);
 %!error <wrong type argument 'cell'> max ({1 2 3 4})
 %!error <cannot compute max \(cell, scalar\)> max ({1, 2, 3}, 2)
@@ -914,11 +1151,25 @@ do_cumminmax_body (const octave_value_li
 {
   int nargin = args.length ();
 
+  bool includenan = false;
+
+  const char *func = ismin ? "cummin" : "cummax";
+
+  if (nargin >= 2 && args(nargin - 1).is_string ())
+    {
+      std::string str = args(nargin - 1).string_value ();
+
+      if (str == "includenan")
+        includenan = true;
+      else if (str != "omitnan")
+        error ("%s: unrecognized NaN flag '%s'", func, str.c_str ());
+
+      nargin--;
+    }
+
   if (nargin < 1 || nargin > 2)
     print_usage ();
 
-  const char *func = ismin ? "cummin" : "cummax";
-
   octave_value arg = args(0);
   int dim = -1;
   if (nargin == 2)
@@ -934,12 +1185,19 @@ do_cumminmax_body (const octave_value_li
   switch (arg.builtin_type ())
     {
     case btyp_double:
-      retval = do_cumminmax_red_op<NDArray> (arg, nargout, dim, ismin);
+      if (arg.is_sparse_type ())
+        retval = do_cumminmax_red_op<SparseMatrix> (arg, nargout, dim, ismin);
+      else
+        retval = do_cumminmax_red_op<NDArray> (arg, nargout, dim, ismin);
       break;
 
     case btyp_complex:
-      retval = do_cumminmax_red_op<ComplexNDArray> (arg, nargout, dim,
-                                                    ismin);
+      if (arg.is_sparse_type ())
+        retval = do_cumminmax_red_op<SparseComplexMatrix> (arg, nargout, dim,
+                                                           ismin);
+      else
+        retval = do_cumminmax_red_op<ComplexNDArray> (arg, nargout, dim,
+                                                      ismin);
       break;
 
     case btyp_float:
@@ -980,6 +1238,141 @@ do_cumminmax_body (const octave_value_li
       err_wrong_type_arg (func, arg);
     }
 
+  if (dim < 0)
+    dim = arg.dims ().first_non_singleton ();
+
+  if (includenan && btyp_isfloat (arg.builtin_type ()) && dim < arg.ndims ()
+      && arg.dims ()(dim) > 0)
+    {
+      Array<octave_idx_type> cummax_idx, max_idx;
+      idx_vector nans_out_idx, nans_dim_idx;
+      int8NDArray nans_out_full, nans_dim_full;
+      SparseMatrix nans_out_sparse, nans_dim_sparse;
+
+      bool extract_first_nans = false;
+      bool full = ! arg.is_sparse_type ();
+      int nd = arg.ndims ();
+
+      if (full)
+        {
+          int8NDArray nans_in = arg.isnan ().int8_array_value ();
+          nans_out_full = nans_in.cummax (cummax_idx, dim);
+          nans_out_idx = octave_value (nans_out_full).bool_array_value ();
+          extract_first_nans = true;
+        }
+      else
+        {
+          SparseMatrix nans_in = arg.isnan ().sparse_matrix_value ();
+          if (nargout > 1)
+            {
+              nans_out_sparse = nans_in.cummax (cummax_idx, dim);
+              extract_first_nans = true;
+            }
+          else
+            {
+              nans_out_sparse = nans_in.cummax (dim);
+              nans_dim_sparse = nans_in.max (max_idx, dim);
+              nans_dim_idx =
+                octave_value (nans_dim_sparse).sparse_bool_matrix_value ();
+            }
+          nans_out_idx =
+            octave_value (nans_out_sparse).sparse_bool_matrix_value ();
+        }
+
+      if (extract_first_nans)
+        {
+          Array<idx_vector> last (dim_vector (nd, 1));
+          for (octave_idx_type i = 0; i < nd; i++)
+            {
+              if (i == dim)
+                last(i) = idx_vector (arg.dims ()(i) - 1);
+              else
+                last(i) = idx_vector::colon;
+            }
+
+          if (full)
+            {
+              nans_dim_full = nans_out_full.index (last);
+              nans_dim_idx = octave_value (nans_dim_full).bool_array_value ();
+              max_idx = cummax_idx.index (last);
+            }
+          else
+            {
+              nans_dim_sparse = nans_out_sparse.index (last(0), last(1));
+              nans_dim_idx =
+                octave_value (nans_dim_sparse).sparse_bool_matrix_value ();
+              max_idx = cummax_idx.index (last(0), last(1));
+            }
+        }
+
+      Array<idx_vector> idxa = ind2sub (max_idx.dims ().redim (nd),
+                                        nans_dim_idx);
+      idxa(dim) = max_idx.index (nans_dim_idx);
+      idx_vector first_nan_idx = sub2ind (arg.dims (), idxa);
+      octave_value first_nans = arg.do_index_op (ovl (first_nan_idx));
+
+      octave_value nan_values;
+      if (full)
+        {
+          nan_values = NDArray (max_idx.dims ());
+          std::list<octave_value_list> nans_dim_idx_l = {ovl (nans_dim_idx)};
+          nan_values = nan_values.subsasgn ("(", nans_dim_idx_l, first_nans);
+
+          dim_vector dv (1, 1);
+          dv.resize (dim + 1, 1);
+          dv(dim) = arg.dims ()(dim);
+          dv.chop_trailing_singletons ();
+          nan_values = do_binary_op (octave_value::op_add, nan_values,
+                                     NDArray (dv, 0.));
+        }
+      else
+        {
+          nans_dim_sparse.resize (arg.dims ());
+          idx_vector nans_dim_ext =
+            octave_value (nans_dim_sparse).sparse_bool_matrix_value ();
+          SparseMatrix nan_stripes = nans_dim_sparse.cummax (dim);
+          idx_vector stripes_idx =
+            octave_value (nan_stripes).sparse_bool_matrix_value ();
+
+          octave_value generic_nan;
+          if (arg.builtin_type () == btyp_double
+              || arg.builtin_type () == btyp_complex)
+            generic_nan = octave::numeric_limits<double>::NaN ();
+          else
+            generic_nan = octave::numeric_limits<float>::NaN ();
+
+          std::list<octave_value_list> stripes_idx_l = {ovl (stripes_idx)};
+          std::list<octave_value_list> nans_dim_ext_l = {ovl (nans_dim_ext)};
+
+          nan_values = SparseMatrix (arg.dims ());
+          nan_values = nan_values.subsasgn ("(", stripes_idx_l, generic_nan);
+          nan_values = nan_values.subsasgn ("(", nans_dim_ext_l, first_nans);
+
+          if (nan_values.is_complex_matrix ())
+            {
+              SparseComplexMatrix nsp =
+                octave_value_extract<SparseComplexMatrix> (nan_values);
+              nan_values = nsp.cummax (dim);
+            }
+          else
+            {
+              SparseMatrix nsp =
+                octave_value_extract<SparseMatrix> (nan_values);
+              nan_values = nsp.cummax (dim);
+            }
+        }
+
+      octave_value r_rhs = nan_values.do_index_op (ovl (nans_out_idx));
+      std::list<octave_value_list> nans_out_idx_l = {ovl (nans_out_idx)};
+      retval(0) = retval(0).subsasgn ("(", nans_out_idx_l, r_rhs);
+
+      if (nargout > 1)
+        {
+          idx_vector i_rhs = cummax_idx.index (nans_out_idx);
+          retval(1) = retval(1).subsasgn ("(", nans_out_idx_l, i_rhs);
+        }
+    }
+
   return retval;
 }
 
@@ -988,6 +1381,8 @@ DEFUN (cummin, args, nargout,
 @deftypefn  {} {} cummin (@var{x})
 @deftypefnx {} {} cummin (@var{x}, @var{dim})
 @deftypefnx {} {[@var{w}, @var{iw}] =} cummin (@var{x})
+@deftypefnx {} {} cummin (@dots{}, "omitnan")
+@deftypefnx {} {} cummin (@dots{}, "includenan")
 Return the cumulative minimum values along dimension @var{dim}.
 
 If @var{dim} is unspecified it defaults to column-wise operation.  For
@@ -1012,6 +1407,9 @@ iw = 1  2  2  4  4  6
 @end group
 @end example
 
+The optional @qcode{"nanflag"} argument specifies if NaN values found in the
+input array are to be ignored or propagated in the result.
+
 @seealso{cummax, min, max}
 @end deftypefn */)
 {
@@ -1034,6 +1432,21 @@ iw = 1  2  2  4  4  6
 %! assert (ndims (iw), 3);
 %! assert (iw, ones (2,2,2));
 
+## Test nanflag argument
+%!assert (cummin ([NA, 3, 1, 2], "omitnan"), [NA, 3, 1, 1])
+%!test
+%! [w, iw] = cummin (single ([3, 1, NA, 2, NaN]), "includenan");
+%! assert (w, single ([3, 1, NA, NA, NA]));
+%! assert (iw, [1, 2, 3, 3, 3]);
+%!test
+%! x = sparse ([NaN+i, 2, 0; NaN+2i, 0, NA]);
+%! assert (cummin (x, "includenan"), sparse ([NaN+i, 2, 0; NaN+i, 0, NA]));
+%!test
+%! x = reshape ([1, 2, 3, NaN, NaN, NA, 4, NaN], [2,2,2]);
+%! [w, iw] = cummin (x, 2, "includenan");
+%! assert (w, reshape ([1, 2, 1, NaN, NaN, NA, NaN, NA], [2,2,2]));
+%! assert (iw, reshape ([1, 1, 1, 2, 1, 1, 1, 1], [2,2,2]));
+
 %!error cummin ()
 %!error cummin (1, 2, 3)
 */
@@ -1043,6 +1456,8 @@ DEFUN (cummax, args, nargout,
 @deftypefn  {} {} cummax (@var{x})
 @deftypefnx {} {} cummax (@var{x}, @var{dim})
 @deftypefnx {} {[@var{w}, @var{iw}] =} cummax (@dots{})
+@deftypefnx {} {} cummax (@dots{}, "omitnan")
+@deftypefnx {} {} cummax (@dots{}, "includenan")
 Return the cumulative maximum values along dimension @var{dim}.
 
 If @var{dim} is unspecified it defaults to column-wise operation.  For
@@ -1067,6 +1482,9 @@ iw = 1  2  2  4  4  4
 @end group
 @end example
 
+The optional @qcode{"nanflag"} argument specifies if NaN values found in the
+input array are to be ignored or propagated in the result.
+
 @seealso{cummin, max, min}
 @end deftypefn */)
 {
@@ -1089,6 +1507,21 @@ iw = 1  2  2  4  4  4
 %! assert (ndims (iw), 3);
 %! assert (iw, ones (2,2,2));
 
+## Test nanflag argument
+%!assert (cummax ([NA, 3, 1, 2], "omitnan"), [NA, 3, 3, 3])
+%!test
+%! [w, iw] = cummax (single ([3, 1, NA, 2, NaN]), "includenan");
+%! assert (w, single ([3, 3, NA, NA, NA]));
+%! assert (iw, [1, 1, 3, 3, 3]);
+%!test
+%! x = sparse ([NaN+i, 2, 0; NaN+2i, 0, NA]);
+%! assert (cummax (x, "includenan"), sparse ([NaN+i, 2, 0; NaN+i, 2, NA]));
+%!test
+%! x = reshape ([1, 2, 3, NaN, NaN, NA, 4, NaN], [2,2,2]);
+%! [w, iw] = cummax (x, 2, "includenan");
+%! assert (w, reshape ([1, 2, 3, NaN, NaN, NA, NaN, NA], [2,2,2]));
+%! assert (iw, reshape ([1, 1, 2, 2, 1, 1, 1, 1], [2,2,2]));
+
 %!error cummax ()
 %!error cummax (1, 2, 3)
 */
diff --git a/libinterp/operators/op-s-m.cc b/libinterp/operators/op-s-m.cc
--- a/libinterp/operators/op-s-m.cc
+++ b/libinterp/operators/op-s-m.cc
@@ -29,6 +29,7 @@ along with Octave; see the file COPYING.
 #include "ov.h"
 #include "ov-scalar.h"
 #include "ov-float.h"
+#include "ov-flt-complex.h"
 #include "ov-re-mat.h"
 #include "ov-flt-re-mat.h"
 #include "ov-typeinfo.h"
@@ -134,6 +135,8 @@ install_s_m_ops (void)
 
   INSTALL_ASSIGNCONV (octave_scalar, octave_matrix, octave_matrix);
   INSTALL_ASSIGNCONV (octave_float_scalar, octave_matrix, octave_float_matrix);
+  INSTALL_ASSIGNCONV (octave_float_complex_scalar, octave_matrix,
+                      octave_float_complex_matrix);
 
   INSTALL_WIDENOP (octave_scalar, octave_matrix, matrix_conv);
 }
diff --git a/libinterp/operators/op-s-s.cc b/libinterp/operators/op-s-s.cc
--- a/libinterp/operators/op-s-s.cc
+++ b/libinterp/operators/op-s-s.cc
@@ -29,8 +29,10 @@ along with Octave; see the file COPYING.
 #include "errwarn.h"
 #include "ovl.h"
 #include "ov.h"
+#include "ov-complex.h"
 #include "ov-scalar.h"
 #include "ov-float.h"
+#include "ov-flt-complex.h"
 #include "ov-re-mat.h"
 #include "ov-flt-re-mat.h"
 #include "ov-typeinfo.h"
@@ -174,4 +176,9 @@ install_s_s_ops (void)
   INSTALL_ASSIGNCONV (octave_scalar, octave_null_matrix, octave_matrix);
   INSTALL_ASSIGNCONV (octave_scalar, octave_null_str, octave_matrix);
   INSTALL_ASSIGNCONV (octave_scalar, octave_null_sq_str, octave_matrix);
+
+  INSTALL_ASSIGNCONV (octave_float_complex_scalar, octave_scalar,
+                      octave_float_complex_matrix);
+  INSTALL_ASSIGNCONV (octave_float_complex_scalar, octave_complex_scalar,
+                      octave_float_complex_matrix);
 }
diff --git a/liboctave/array/CMatrix.cc b/liboctave/array/CMatrix.cc
--- a/liboctave/array/CMatrix.cc
+++ b/liboctave/array/CMatrix.cc
@@ -3641,6 +3641,9 @@ min (const Complex& c, const ComplexMatr
 ComplexMatrix
 min (const ComplexMatrix& m, const Complex& c)
 {
+  if (octave::math::isnan (c))
+    return m;
+
   return min (c, m);
 }
 
@@ -3659,33 +3662,11 @@ min (const ComplexMatrix& a, const Compl
   ComplexMatrix result (nr, nc);
 
   for (octave_idx_type j = 0; j < nc; j++)
-    {
-      bool columns_are_real_only = true;
-      for (octave_idx_type i = 0; i < nr; i++)
-        {
-          octave_quit ();
-          if (octave::math::imag (a(i, j)) != 0.0 || octave::math::imag (b(i, j)) != 0.0)
-            {
-              columns_are_real_only = false;
-              break;
-            }
-        }
-
-      if (columns_are_real_only)
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            result(i, j) = octave::math::min (octave::math::real (a(i, j)),
-                                              octave::math::real (b(i, j)));
-        }
-      else
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            {
-              octave_quit ();
-              result(i, j) = octave::math::min (a(i, j), b(i, j));
-            }
-        }
-    }
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+        octave_quit ();
+        result(i, j) = octave::math::min (a(i, j), b(i, j));
+      }
 
   return result;
 }
@@ -3713,6 +3694,9 @@ max (const Complex& c, const ComplexMatr
 ComplexMatrix
 max (const ComplexMatrix& m, const Complex& c)
 {
+  if (octave::math::isnan (c))
+    return m;
+
   return max (c, m);
 }
 
@@ -3731,37 +3715,11 @@ max (const ComplexMatrix& a, const Compl
   ComplexMatrix result (nr, nc);
 
   for (octave_idx_type j = 0; j < nc; j++)
-    {
-      bool columns_are_real_only = true;
-      for (octave_idx_type i = 0; i < nr; i++)
-        {
-          octave_quit ();
-          if (octave::math::imag (a(i, j)) != 0.0 || octave::math::imag (b(i, j)) != 0.0)
-            {
-              columns_are_real_only = false;
-              break;
-            }
-        }
-
-      // FIXME: is it so much faster?
-      if (columns_are_real_only)
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            {
-              octave_quit ();
-              result(i, j) = octave::math::max (octave::math::real (a(i, j)),
-                                                octave::math::real (b(i, j)));
-            }
-        }
-      else
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            {
-              octave_quit ();
-              result(i, j) = octave::math::max (a(i, j), b(i, j));
-            }
-        }
-    }
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+        octave_quit ();
+        result(i, j) = octave::math::max (a(i, j), b(i, j));
+      }
 
   return result;
 }
diff --git a/liboctave/array/CSparse.cc b/liboctave/array/CSparse.cc
--- a/liboctave/array/CSparse.cc
+++ b/liboctave/array/CSparse.cc
@@ -217,8 +217,7 @@ SparseComplexMatrix::max (Array<octave_i
       octave_idx_type nel = 0;
       for (octave_idx_type j = 0; j < nc; j++)
         {
-          Complex tmp_max;
-          double abs_max = octave::numeric_limits<double>::NaN ();
+          Complex tmp_max = Complex_NaN_result;
           octave_idx_type idx_j = 0;
           for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
             {
@@ -229,10 +228,7 @@ SparseComplexMatrix::max (Array<octave_i
             }
 
           if (idx_j != nr)
-            {
-              tmp_max = 0.;
-              abs_max = 0.;
-            }
+            tmp_max = 0.;
 
           for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
             {
@@ -241,18 +237,15 @@ SparseComplexMatrix::max (Array<octave_i
               if (octave::math::isnan (tmp))
                 continue;
 
-              double abs_tmp = std::abs (tmp);
-
-              if (octave::math::isnan (abs_max) || abs_tmp > abs_max)
+              if (octave::math::isnan (tmp_max) || tmp > tmp_max)
                 {
                   idx_j = ridx (i);
                   tmp_max = tmp;
-                  abs_max = abs_tmp;
                 }
             }
 
           idx_arg.elem (j) = octave::math::isnan (tmp_max) ? 0 : idx_j;
-          if (abs_max != 0.)
+          if (tmp_max != 0.)
             nel++;
         }
 
@@ -299,10 +292,11 @@ SparseComplexMatrix::max (Array<octave_i
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               Complex tmp = data (i);
+              Complex best = elem (ir, ix);
 
               if (octave::math::isnan (tmp))
                 continue;
-              else if (ix == -1 || std::abs (tmp) > std::abs (elem (ir, ix)))
+              else if (octave::math::isnan (best) || tmp > best)
                 idx_arg.elem (ir) = j;
             }
         }
@@ -374,8 +368,7 @@ SparseComplexMatrix::min (Array<octave_i
       octave_idx_type nel = 0;
       for (octave_idx_type j = 0; j < nc; j++)
         {
-          Complex tmp_min;
-          double abs_min = octave::numeric_limits<double>::NaN ();
+          Complex tmp_min = Complex_NaN_result;
           octave_idx_type idx_j = 0;
           for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
             {
@@ -386,10 +379,7 @@ SparseComplexMatrix::min (Array<octave_i
             }
 
           if (idx_j != nr)
-            {
-              tmp_min = 0.;
-              abs_min = 0.;
-            }
+            tmp_min = 0.;
 
           for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
             {
@@ -398,18 +388,15 @@ SparseComplexMatrix::min (Array<octave_i
               if (octave::math::isnan (tmp))
                 continue;
 
-              double abs_tmp = std::abs (tmp);
-
-              if (octave::math::isnan (abs_min) || abs_tmp < abs_min)
+              if (octave::math::isnan (tmp_min) || tmp < tmp_min)
                 {
                   idx_j = ridx (i);
                   tmp_min = tmp;
-                  abs_min = abs_tmp;
                 }
             }
 
           idx_arg.elem (j) = octave::math::isnan (tmp_min) ? 0 : idx_j;
-          if (abs_min != 0.)
+          if (tmp_min != 0.)
             nel++;
         }
 
@@ -456,10 +443,11 @@ SparseComplexMatrix::min (Array<octave_i
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               Complex tmp = data (i);
+              Complex best = elem (ir, ix);
 
               if (octave::math::isnan (tmp))
                 continue;
-              else if (ix == -1 || std::abs (tmp) < std::abs (elem (ir, ix)))
+              else if (octave::math::isnan (best) || tmp < best)
                 idx_arg.elem (ir) = j;
             }
         }
@@ -512,6 +500,731 @@ SparseComplexMatrix::min (Array<octave_i
 
 */
 
+SparseComplexMatrix
+SparseComplexMatrix::cummax (int dim) const
+{
+  SparseComplexMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    return *this;
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseComplexMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type expected_r = 0;
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (! octave::math::isnan (data (i)))
+                    break;
+                  expected_r++;
+                }
+
+              if (i < cidx (j+1) && ridx (i) == expected_r)
+                {
+                  nel += nr;
+                  continue;
+                }
+              nel += expected_r;
+            }
+          for (; i < cidx (j+1); i++)
+            if (! octave::math::isnan (data (i)))
+              {
+                nel += nr - ridx (i);
+                break;
+              }
+        }
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1))
+            {
+              octave_idx_type last = 0;
+              bool zero_found = ridx (i) > 0;
+              Complex tmp_max = zero_found ? 0. : data (i++);
+
+              if (! zero_found
+                  && ((i < cidx (j+1) && ridx (i) > 1)
+                      || (i == cidx (j+1) && nr >= 2)))
+                {
+                  zero_found = true;
+                  if (octave::math::isnan (tmp_max))
+                    {
+                      result.xdata (ii) = tmp_max;
+                      result.xridx (ii++) = 0;
+                      last = 1;
+                      tmp_max = 0.;
+                    }
+                }
+
+              for (; i < cidx (j+1); i++)
+                {
+                  if (data (i) > tmp_max
+                      || (! octave::math::isnan (data (i))
+                          && octave::math::isnan (tmp_max)))
+                    {
+                      if (tmp_max != 0.)
+                        for (; last < ridx (i); last++)
+                          {
+                            result.xdata (ii) = tmp_max;
+                            result.xridx (ii++) = last;
+                          }
+                      last = ridx (i);
+                      tmp_max = data (i);
+                    }
+                  if (! zero_found
+                      && ((i + 1 < cidx (j+1) && ridx (i+1) > ridx (i) + 1)
+                          || (i + 1 == cidx (j+1) && ridx (i) < nr - 1)))
+                    {
+                      zero_found = true;
+                      if (octave::math::isnan (tmp_max))
+                        {
+                          for (; last < ridx (i) + 1; last++)
+                            {
+                              result.xdata (ii) = tmp_max;
+                              result.xridx (ii++) = last;
+                            }
+                          last = ridx (i) + 1;
+                          tmp_max = 0.;
+                        }
+                    }
+                }
+              if (tmp_max != 0.)
+                for (; last < nr; last++)
+                  {
+                    result.xdata (ii) = tmp_max;
+                    result.xridx (ii++) = last;
+                  }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (int8_t, max_sgn, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          max_sgn[i] = -1;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (max_sgn[ridx (i)] <= 0)
+              {
+                bool end_of_zeros = ! octave::math::isnan (data (i));
+
+                if (max_sgn[ridx (i)] < 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || end_of_zeros)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          max_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (end_of_zeros)
+                  {
+                    max_sgn[ridx (i)] = 1;
+                    nel += nc - j;
+                  }
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (max_sgn[i] < 0)
+          nel += expected_visit[i];
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type i = cidx (0); i < cidx (1); i++)
+        {
+          result.xdata (ii) = data (i);
+          result.xridx (ii++) = ridx (i);
+        }
+      result.xcidx (1) = ii;
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          octave_idx_type i_max = cidx (j+1);
+          bool i_lt_max = i < i_max;
+
+          octave_idx_type r = result.cidx (j-1);
+          octave_idx_type r_max = result.cidx (j);
+          bool r_lt_max = r < r_max;
+
+          while (i_lt_max || r_lt_max)
+            {
+              octave_quit ();
+              Complex i_val, r_val;
+              octave_idx_type row;
+
+              if (! r_lt_max || (i_lt_max && ridx (i) < result.ridx (r)))
+                {
+                  i_val = data (i);
+                  r_val = 0.;
+                  row = ridx (i++);
+                  i_lt_max = i < i_max;
+                }
+              else if (! i_lt_max || (r_lt_max && result.ridx (r) < ridx (i)))
+                {
+                  i_val = 0.;
+                  r_val = result.data (r);
+                  row = result.ridx (r++);
+                  r_lt_max = r < r_max;
+                }
+              else
+                {
+                  i_val = data (i);
+                  r_val = result.data (r);
+                  row = ridx (i++);
+                  i_lt_max = i < i_max;
+                  r++;
+                  r_lt_max = r < r_max;
+                }
+
+              Complex tmp = octave::math::max (r_val, i_val);
+              if (tmp != 0.)
+                {
+                  result.data (ii) = tmp;
+                  result.ridx (ii++) = row;
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+{
+  SparseComplexMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    {
+      idx_arg.resize (dim_vector (nr, nc), 0);
+      return *this;
+    }
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  idx_arg.resize (dim_vector (nr, nc));
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseComplexMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type expected_r = 0;
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (! octave::math::isnan (data (i)))
+                    break;
+                  expected_r++;
+                }
+
+              if (i < cidx (j+1) && ridx (i) == expected_r)
+                {
+                  nel += nr;
+                  continue;
+                }
+              nel += expected_r;
+            }
+          for (; i < cidx (j+1); i++)
+            if (! octave::math::isnan (data (i)))
+              {
+                nel += nr - ridx (i);
+                break;
+              }
+        }
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      octave_idx_type i = cidx (0);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp_max = i < cidx (j+1) && ridx (i) == 0 ? data (i++) : 0.;
+          idx_arg.elem (0, j) = 0;
+          if (tmp_max != 0.)
+            {
+              result.xdata (ii) = tmp_max;
+              result.xridx (ii++) = 0;
+            }
+
+          for (octave_idx_type ir = 1; ir < nr; ir++)
+            {
+              Complex tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+
+              if (tmp > tmp_max
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_max)))
+                {
+                  tmp_max = tmp;
+                  idx_arg.elem (ir, j) = ir;
+                }
+              else
+                idx_arg.elem (ir, j) = idx_arg.elem (ir - 1, j);
+              if (tmp_max != 0.)
+                {
+                  result.xdata (ii) = tmp_max;
+                  result.xridx (ii++) = ir;
+                }
+            }
+
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (int8_t, max_sgn, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          max_sgn[i] = -1;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (max_sgn[ridx (i)] <= 0)
+              {
+                bool end_of_zeros = ! octave::math::isnan (data (i));
+
+                if (max_sgn[ridx (i)] < 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || end_of_zeros)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          max_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (end_of_zeros)
+                  {
+                    max_sgn[ridx (i)] = 1;
+                    nel += nc - j;
+                  }
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (max_sgn[i] < 0)
+          nel += expected_visit[i];
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      octave_idx_type i = cidx (0);
+      for (octave_idx_type ir = 0; ir < nr; ir++)
+        {
+          idx_arg.elem (ir, 0) = 0;
+          if (i < cidx (1) && ridx (i) == ir)
+            {
+              result.xdata (ii) = data (i++);
+              result.xridx (ii++) = ir;
+            }
+        }
+      result.xcidx (1) = ii;
+      octave_idx_type r = result.cidx (0);
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          for (octave_idx_type ir = 0; ir < nr; ir++)
+            {
+              Complex tmp_max = (r < result.cidx (j) && result.ridx (r) == ir
+                                 ? result.data (r++) : 0.);
+              Complex tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+              if (tmp > tmp_max
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_max)))
+                {
+                  idx_arg.elem (ir, j) = j;
+                  if (tmp != 0.)
+                    {
+                      result.xdata (ii) = tmp;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+              else
+                {
+                  idx_arg.elem (ir, j) = idx_arg.elem (ir, j - 1);
+                  if (tmp_max != 0.)
+                    {
+                      result.xdata (ii) = tmp_max;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::cummin (int dim) const
+{
+  SparseComplexMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    return *this;
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseComplexMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type expected_r = 0;
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              if (ridx (i) != expected_r)
+                break;
+              expected_r++;
+            }
+          nel += expected_r;
+        }
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type last = 0;
+              Complex tmp_min = data (i++);
+              octave_idx_type expected_r = 1;
+
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (data (i) < tmp_min
+                      || (! octave::math::isnan (data (i))
+                          && octave::math::isnan (tmp_min)))
+                    {
+                      for (; last < ridx (i); last++)
+                        {
+                          result.xdata (ii) = tmp_min;
+                          result.xridx (ii++) = last;
+                        }
+                      last = ridx (i);
+                      tmp_min = data (i);
+                    }
+                  expected_r++;
+                }
+              for (; last < expected_r; last++)
+                {
+                  result.xdata (ii) = tmp_min;
+                  result.xridx (ii++) = last;
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (bool, leading_nonzeros, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          leading_nonzeros[i] = true;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (leading_nonzeros[ridx (i)])
+              {
+                if (expected_visit[ridx (i)] < j)
+                  {
+                    nel += expected_visit[ridx (i)];
+                    leading_nonzeros[ridx (i)] = false;
+                  }
+                else
+                  expected_visit[ridx (i)]++;
+              }
+          }
+
+      octave_idx_type max_leading_nonzeros = 0;
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          if (leading_nonzeros[i])
+            nel += expected_visit[i];
+          if (expected_visit[i] > max_leading_nonzeros)
+            max_leading_nonzeros = expected_visit[i];
+        }
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type i = cidx (0); i < cidx (1); i++)
+        {
+          result.xdata (ii) = data (i);
+          result.xridx (ii++) = ridx (i);
+          expected_visit[ridx (i)]--;
+        }
+      result.xcidx (1) = ii;
+
+      for (octave_idx_type j = 1; j < max_leading_nonzeros; j++)
+        {
+          octave_idx_type r = result.cidx (j-1);
+
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            if (expected_visit[ridx (i)] > 0)
+              {
+                expected_visit[ridx (i)]--;
+
+                for (; result.ridx (r) < ridx (i); r++);
+                result.data (ii) = octave::math::min (result.data (r),
+                                                      data (i));
+                result.ridx (ii++) = ridx (i);
+              }
+
+          result.xcidx (j+1) = ii;
+        }
+
+      for (octave_idx_type j = max_leading_nonzeros + 1; j <= nc; j++)
+        result.xcidx (j) = ii;
+    }
+
+  return result;
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+{
+  SparseComplexMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    {
+      idx_arg.resize (dim_vector (nr, nc), 0);
+      return *this;
+    }
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  idx_arg.resize (dim_vector (nr, nc));
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseComplexMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type expected_r = 0;
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              if (ridx (i) != expected_r)
+                break;
+              expected_r++;
+            }
+          nel += expected_r;
+        }
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      octave_idx_type i = cidx (0);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp_min = i < cidx (j+1) && ridx (i) == 0 ? data (i++) : 0.;
+          idx_arg.elem (0, j) = 0;
+          if (tmp_min != 0.)
+            {
+              result.xdata (ii) = tmp_min;
+              result.xridx (ii++) = 0;
+            }
+
+          for (octave_idx_type ir = 1; ir < nr; ir++)
+            {
+              Complex tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+
+              if (tmp < tmp_min
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_min)))
+                {
+                  tmp_min = tmp;
+                  idx_arg.elem (ir, j) = ir;
+                }
+              else
+                idx_arg.elem (ir, j) = idx_arg.elem (ir - 1, j);
+              if (tmp_min != 0.)
+                {
+                  result.xdata (ii) = tmp_min;
+                  result.xridx (ii++) = ir;
+                }
+            }
+
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (bool, leading_nonzeros, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          leading_nonzeros[i] = true;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (leading_nonzeros[ridx (i)])
+              {
+                if (expected_visit[ridx (i)] < j)
+                  {
+                    nel += expected_visit[ridx (i)];
+                    leading_nonzeros[ridx (i)] = false;
+                  }
+                else
+                  expected_visit[ridx (i)]++;
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (leading_nonzeros[i])
+          nel += expected_visit[i];
+
+      result = SparseComplexMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      octave_idx_type i = cidx (0);
+      for (octave_idx_type ir = 0; ir < nr; ir++)
+        {
+          idx_arg.elem (ir, 0) = 0;
+          if (i < cidx (1) && ridx (i) == ir)
+            {
+              result.xdata (ii) = data (i++);
+              result.xridx (ii++) = ir;
+            }
+        }
+      result.xcidx (1) = ii;
+      octave_idx_type r = result.cidx (0);
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          for (octave_idx_type ir = 0; ir < nr; ir++)
+            {
+              Complex tmp_min = (r < result.cidx (j) && result.ridx (r) == ir
+                                 ? result.data (r++) : 0.);
+              Complex tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+              if (tmp < tmp_min
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_min)))
+                {
+                  idx_arg.elem (ir, j) = j;
+                  if (tmp != 0.)
+                    {
+                      result.xdata (ii) = tmp;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+              else
+                {
+                  idx_arg.elem (ir, j) = idx_arg.elem (ir, j - 1);
+                  if (tmp_min != 0.)
+                    {
+                      result.xdata (ii) = tmp_min;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+/*
+
+%!test
+%! x = kron (eye (4), sparse (1, 1, 1, 2^14, 2^14) * i);
+%! assert (sum (sum (cummax (x, 1))), sparse ((2^17 + 2^15) * i));
+%! assert (sum (sum (cummax (x, 2))), sparse ((2^17 + 2^15) * i));
+%! assert (sum (sum (cummin (x, 1))), sparse (i));
+%! assert (sum (sum (cummin (x, 2))), sparse (i));
+
+*/
+
 ComplexRowVector
 SparseComplexMatrix::row (octave_idx_type i) const
 {
@@ -7694,6 +8407,9 @@ min (const Complex& c, const SparseCompl
 SparseComplexMatrix
 min (const SparseComplexMatrix& m, const Complex& c)
 {
+  if (octave::math::isnan (c))
+    return m;
+
   return min (c, m);
 }
 
@@ -7795,7 +8511,15 @@ max (const Complex& c, const SparseCompl
   EMPTY_RETURN_CHECK (SparseComplexMatrix);
 
   // Count the number of nonzero elements
-  if (octave::math::max (c, 0.) != 0.)
+  if (octave::math::isnan (c))
+    {
+      result = SparseComplexMatrix (m);
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+          if (octave::math::isnan (m.data (i)))
+            result.data (i) = c;
+    }
+  else if (octave::math::max (c, 0.) != 0.)
     {
       result = SparseComplexMatrix (nr, nc, c);
       for (octave_idx_type j = 0; j < nc; j++)
@@ -7803,7 +8527,24 @@ max (const Complex& c, const SparseCompl
           result.xdata (m.ridx (i) + j * nr) = octave::math::max (c, m.data (i));
     }
   else
-    result = SparseComplexMatrix (m);
+    {
+      result = SparseComplexMatrix (nr, nc, m.nnz ());
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+            if (! octave::math::isnan (m.data (i)))
+              {
+                result.xdata (ii) = m.data (i);
+                result.xridx (ii++) = m.ridx (i);
+              }
+          result.xcidx (j+1) = ii;
+        }
+
+      result.maybe_compress ();
+    }
 
   return result;
 }
@@ -7811,6 +8552,9 @@ max (const Complex& c, const SparseCompl
 SparseComplexMatrix
 max (const SparseComplexMatrix& m, const Complex& c)
 {
+  if (octave::math::isnan (c))
+    return m;
+
   return max (c, m);
 }
 
diff --git a/liboctave/array/CSparse.h b/liboctave/array/CSparse.h
--- a/liboctave/array/CSparse.h
+++ b/liboctave/array/CSparse.h
@@ -121,6 +121,13 @@ public:
   SparseComplexMatrix min (int dim = -1) const;
   SparseComplexMatrix min (Array<octave_idx_type>& index, int dim = -1) const;
 
+  SparseComplexMatrix cummax (int dim = -1) const;
+  SparseComplexMatrix cummax (Array<octave_idx_type>& index,
+                              int dim = -1) const;
+  SparseComplexMatrix cummin (int dim = -1) const;
+  SparseComplexMatrix cummin (Array<octave_idx_type>& index,
+                              int dim = -1) const;
+
   SparseComplexMatrix& insert (const SparseComplexMatrix& a,
                                octave_idx_type r, octave_idx_type c);
   SparseComplexMatrix& insert (const SparseMatrix& a,
diff --git a/liboctave/array/dMatrix.cc b/liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc
+++ b/liboctave/array/dMatrix.cc
@@ -3048,21 +3048,10 @@ min (double d, const Matrix& m)
 Matrix
 min (const Matrix& m, double d)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (Matrix);
-
-  Matrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::min (m(i, j), d);
-      }
-
-  return result;
+  if (octave::math::isnan (d))
+    return m;
+
+  return min (d, m);
 }
 
 Matrix
@@ -3112,21 +3101,10 @@ max (double d, const Matrix& m)
 Matrix
 max (const Matrix& m, double d)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (Matrix);
-
-  Matrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::max (m(i, j), d);
-      }
-
-  return result;
+  if (octave::math::isnan (d))
+    return m;
+
+  return max (d, m);
 }
 
 Matrix
diff --git a/liboctave/array/dSparse.cc b/liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc
+++ b/liboctave/array/dSparse.cc
@@ -295,10 +295,11 @@ SparseMatrix::max (Array<octave_idx_type
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               double tmp = data (i);
+              double best = elem (ir, ix);
 
               if (octave::math::isnan (tmp))
                 continue;
-              else if (ix == -1 || tmp > elem (ir, ix))
+              else if (octave::math::isnan (best) || tmp > best)
                 idx_arg.elem (ir) = j;
             }
         }
@@ -446,10 +447,11 @@ SparseMatrix::min (Array<octave_idx_type
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               double tmp = data (i);
+              double best = elem (ir, ix);
 
               if (octave::math::isnan (tmp))
                 continue;
-              else if (ix == -1 || tmp < elem (ir, ix))
+              else if (octave::math::isnan (best) || tmp < best)
                 idx_arg.elem (ir) = j;
             }
         }
@@ -502,6 +504,831 @@ SparseMatrix::min (Array<octave_idx_type
 
 */
 
+SparseMatrix
+SparseMatrix::cummax (int dim) const
+{
+  SparseMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    return *this;
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type expected_r = 0;
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (data (i) > 0.)
+                    break;
+                  expected_r++;
+                }
+
+              if (i < cidx (j+1) && ridx (i) == expected_r)
+                {
+                  nel += nr;
+                  continue;
+                }
+              nel += expected_r;
+            }
+          for (; i < cidx (j+1); i++)
+            if (data (i) > 0.)
+              {
+                nel += nr - ridx (i);
+                break;
+              }
+        }
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1))
+            {
+              octave_idx_type last = 0;
+              bool zero_found = ridx (i) > 0;
+              double tmp_max = zero_found ? 0. : data (i++);
+
+              if (! zero_found
+                  && ((i < cidx (j+1) && ridx (i) > 1)
+                      || (i == cidx (j+1) && nr >= 2)))
+                {
+                  zero_found = true;
+                  if (0. > tmp_max || octave::math::isnan (tmp_max))
+                    {
+                      result.xdata (ii) = tmp_max;
+                      result.xridx (ii++) = 0;
+                      last = 1;
+                      tmp_max = 0.;
+                    }
+                }
+
+              for (; i < cidx (j+1); i++)
+                {
+                  if (data (i) > tmp_max
+                      || (! octave::math::isnan (data (i))
+                          && octave::math::isnan (tmp_max)))
+                    {
+                      if (tmp_max != 0.)
+                        for (; last < ridx (i); last++)
+                          {
+                            result.xdata (ii) = tmp_max;
+                            result.xridx (ii++) = last;
+                          }
+                      last = ridx (i);
+                      tmp_max = data (i);
+                    }
+                  if (! zero_found
+                      && ((i + 1 < cidx (j+1) && ridx (i+1) > ridx (i) + 1)
+                          || (i + 1 == cidx (j+1) && ridx (i) < nr - 1)))
+                    {
+                      zero_found = true;
+                      if (0. > tmp_max || octave::math::isnan (tmp_max))
+                        {
+                          for (; last < ridx (i) + 1; last++)
+                            {
+                              result.xdata (ii) = tmp_max;
+                              result.xridx (ii++) = last;
+                            }
+                          last = ridx (i) + 1;
+                          tmp_max = 0.;
+                        }
+                    }
+                }
+              if (tmp_max != 0.)
+                for (; last < nr; last++)
+                  {
+                    result.xdata (ii) = tmp_max;
+                    result.xridx (ii++) = last;
+                  }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (int8_t, max_sgn, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          max_sgn[i] = -1;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (max_sgn[ridx (i)] <= 0)
+              {
+                if (max_sgn[ridx (i)] < 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || data (i) > 0.)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          max_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (data (i) > 0.)
+                  {
+                    max_sgn[ridx (i)] = 1;
+                    nel += nc - j;
+                  }
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (max_sgn[i] < 0)
+          nel += expected_visit[i];
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type i = cidx (0); i < cidx (1); i++)
+        {
+          result.xdata (ii) = data (i);
+          result.xridx (ii++) = ridx (i);
+        }
+      result.xcidx (1) = ii;
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          octave_idx_type i_max = cidx (j+1);
+          bool i_lt_max = i < i_max;
+
+          octave_idx_type r = result.cidx (j-1);
+          octave_idx_type r_max = result.cidx (j);
+          bool r_lt_max = r < r_max;
+
+          while (i_lt_max || r_lt_max)
+            {
+              octave_quit ();
+              double i_val, r_val;
+              octave_idx_type row;
+
+              if (! r_lt_max || (i_lt_max && ridx (i) < result.ridx (r)))
+                {
+                  i_val = data (i);
+                  r_val = 0.;
+                  row = ridx (i++);
+                  i_lt_max = i < i_max;
+                }
+              else if (! i_lt_max || (r_lt_max && result.ridx (r) < ridx (i)))
+                {
+                  i_val = 0.;
+                  r_val = result.data (r);
+                  row = result.ridx (r++);
+                  r_lt_max = r < r_max;
+                }
+              else
+                {
+                  i_val = data (i);
+                  r_val = result.data (r);
+                  row = ridx (i++);
+                  i_lt_max = i < i_max;
+                  r++;
+                  r_lt_max = r < r_max;
+                }
+
+              double tmp = octave::math::max (r_val, i_val);
+              if (tmp != 0.)
+                {
+                  result.data (ii) = tmp;
+                  result.ridx (ii++) = row;
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseMatrix
+SparseMatrix::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+{
+  SparseMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    {
+      idx_arg.resize (dim_vector (nr, nc), 0);
+      return *this;
+    }
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  idx_arg.resize (dim_vector (nr, nc));
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type expected_r = 0;
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (data (i) > 0.)
+                    break;
+                  expected_r++;
+                }
+
+              if (i < cidx (j+1) && ridx (i) == expected_r)
+                {
+                  nel += nr;
+                  continue;
+                }
+              nel += expected_r;
+            }
+          for (; i < cidx (j+1); i++)
+            if (data (i) > 0.)
+              {
+                nel += nr - ridx (i);
+                break;
+              }
+        }
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      octave_idx_type i = cidx (0);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp_max = i < cidx (j+1) && ridx (i) == 0 ? data (i++) : 0.;
+          idx_arg.elem (0, j) = 0;
+          if (tmp_max != 0.)
+            {
+              result.xdata (ii) = tmp_max;
+              result.xridx (ii++) = 0;
+            }
+
+          for (octave_idx_type ir = 1; ir < nr; ir++)
+            {
+              double tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+
+              if (tmp > tmp_max
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_max)))
+                {
+                  tmp_max = tmp;
+                  idx_arg.elem (ir, j) = ir;
+                }
+              else
+                idx_arg.elem (ir, j) = idx_arg.elem (ir - 1, j);
+              if (tmp_max != 0.)
+                {
+                  result.xdata (ii) = tmp_max;
+                  result.xridx (ii++) = ir;
+                }
+            }
+
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (int8_t, max_sgn, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          max_sgn[i] = -1;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (max_sgn[ridx (i)] <= 0)
+              {
+                if (max_sgn[ridx (i)] < 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || data (i) > 0.)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          max_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (data (i) > 0.)
+                  {
+                    max_sgn[ridx (i)] = 1;
+                    nel += nc - j;
+                  }
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (max_sgn[i] < 0)
+          nel += expected_visit[i];
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      octave_idx_type i = cidx (0);
+      for (octave_idx_type ir = 0; ir < nr; ir++)
+        {
+          idx_arg.elem (ir, 0) = 0;
+          if (i < cidx (1) && ridx (i) == ir)
+            {
+              result.xdata (ii) = data (i++);
+              result.xridx (ii++) = ir;
+            }
+        }
+      result.xcidx (1) = ii;
+      octave_idx_type r = result.cidx (0);
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          for (octave_idx_type ir = 0; ir < nr; ir++)
+            {
+              double tmp_max = (r < result.cidx (j) && result.ridx (r) == ir
+                                ? result.data (r++) : 0.);
+              double tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+              if (tmp > tmp_max
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_max)))
+                {
+                  idx_arg.elem (ir, j) = j;
+                  if (tmp != 0.)
+                    {
+                      result.xdata (ii) = tmp;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+              else
+                {
+                  idx_arg.elem (ir, j) = idx_arg.elem (ir, j - 1);
+                  if (tmp_max != 0.)
+                    {
+                      result.xdata (ii) = tmp_max;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseMatrix
+SparseMatrix::cummin (int dim) const
+{
+  SparseMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    return *this;
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type expected_r = 0;
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (data (i) < 0.)
+                    break;
+                  expected_r++;
+                }
+
+              if (i < cidx (j+1) && ridx (i) == expected_r)
+                {
+                  nel += nr;
+                  continue;
+                }
+              nel += expected_r;
+            }
+          for (; i < cidx (j+1); i++)
+            if (data (i) < 0.)
+              {
+                nel += nr - ridx (i);
+                break;
+              }
+        }
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1))
+            {
+              octave_idx_type last = 0;
+              bool zero_found = ridx (i) > 0;
+              double tmp_min = zero_found ? 0. : data (i++);
+
+              if (! zero_found
+                  && ((i < cidx (j+1) && ridx (i) > 1)
+                      || (i == cidx (j+1) && nr >= 2)))
+                {
+                  zero_found = true;
+                  if (0. < tmp_min || octave::math::isnan (tmp_min))
+                    {
+                      result.xdata (ii) = tmp_min;
+                      result.xridx (ii++) = 0;
+                      last = 1;
+                      tmp_min = 0.;
+                    }
+                }
+
+              for (; i < cidx (j+1); i++)
+                {
+                  if (data (i) < tmp_min
+                      || (! octave::math::isnan (data (i))
+                          && octave::math::isnan (tmp_min)))
+                    {
+                      if (tmp_min != 0.)
+                        for (; last < ridx (i); last++)
+                          {
+                            result.xdata (ii) = tmp_min;
+                            result.xridx (ii++) = last;
+                          }
+                      last = ridx (i);
+                      tmp_min = data (i);
+                    }
+                  if (! zero_found
+                      && ((i + 1 < cidx (j+1) && ridx (i+1) > ridx (i) + 1)
+                          || (i + 1 == cidx (j+1) && ridx (i) < nr - 1)))
+                    {
+                      zero_found = true;
+                      if (0. < tmp_min || octave::math::isnan (tmp_min))
+                        {
+                          for (; last < ridx (i) + 1; last++)
+                            {
+                              result.xdata (ii) = tmp_min;
+                              result.xridx (ii++) = last;
+                            }
+                          last = ridx (i) + 1;
+                          tmp_min = 0.;
+                        }
+                    }
+                }
+              if (tmp_min != 0.)
+                for (; last < nr; last++)
+                  {
+                    result.xdata (ii) = tmp_min;
+                    result.xridx (ii++) = last;
+                  }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (int8_t, min_sgn, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          min_sgn[i] = 1;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (min_sgn[ridx (i)] >= 0)
+              {
+                if (min_sgn[ridx (i)] > 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || data (i) < 0.)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          min_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (data (i) < 0.)
+                  {
+                    min_sgn[ridx (i)] = 1;
+                    nel += nc - j;
+                  }
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (min_sgn[i] > 0)
+          nel += expected_visit[i];
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      for (octave_idx_type i = cidx (0); i < cidx (1); i++)
+        {
+          result.xdata (ii) = data (i);
+          result.xridx (ii++) = ridx (i);
+        }
+      result.xcidx (1) = ii;
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          octave_idx_type i_max = cidx (j+1);
+          bool i_lt_max = i < i_max;
+
+          octave_idx_type r = result.cidx (j-1);
+          octave_idx_type r_max = result.cidx (j);
+          bool r_lt_max = r < r_max;
+
+          while (i_lt_max || r_lt_max)
+            {
+              octave_quit ();
+              double i_val, r_val;
+              octave_idx_type row;
+
+              if (! r_lt_max || (i_lt_max && ridx (i) < result.ridx (r)))
+                {
+                  i_val = data (i);
+                  r_val = 0.;
+                  row = ridx (i++);
+                  i_lt_max = i < i_max;
+                }
+              else if (! i_lt_max || (r_lt_max && result.ridx (r) < ridx (i)))
+                {
+                  i_val = 0.;
+                  r_val = result.data (r);
+                  row = result.ridx (r++);
+                  r_lt_max = r < r_max;
+                }
+              else
+                {
+                  i_val = data (i);
+                  r_val = result.data (r);
+                  row = ridx (i++);
+                  i_lt_max = i < i_max;
+                  r++;
+                  r_lt_max = r < r_max;
+                }
+
+              double tmp = octave::math::min (r_val, i_val);
+              if (tmp != 0.)
+                {
+                  result.data (ii) = tmp;
+                  result.ridx (ii++) = row;
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseMatrix
+SparseMatrix::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+{
+  SparseMatrix result;
+  dim_vector dv = dims ();
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim >= dv.ndims ())
+    {
+      idx_arg.resize (dim_vector (nr, nc), 0);
+      return *this;
+    }
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  idx_arg.resize (dim_vector (nr, nc));
+
+  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
+    return SparseMatrix (nr, nc);
+
+  if (dim == 0)
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_idx_type i = cidx (j);
+          if (i < cidx (j+1) && ridx (i) == 0)
+            {
+              octave_idx_type expected_r = 0;
+              for (; i < cidx (j+1) && ridx (i) == expected_r; i++)
+                {
+                  if (data (i) < 0.)
+                    break;
+                  expected_r++;
+                }
+
+              if (i < cidx (j+1) && ridx (i) == expected_r)
+                {
+                  nel += nr;
+                  continue;
+                }
+              nel += expected_r;
+            }
+          for (; i < cidx (j+1); i++)
+            if (data (i) < 0.)
+              {
+                nel += nr - ridx (i);
+                break;
+              }
+        }
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      octave_idx_type i = cidx (0);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp_min = i < cidx (j+1) && ridx (i) == 0 ? data (i++) : 0.;
+          idx_arg.elem (0, j) = 0;
+          if (tmp_min != 0.)
+            {
+              result.xdata (ii) = tmp_min;
+              result.xridx (ii++) = 0;
+            }
+
+          for (octave_idx_type ir = 1; ir < nr; ir++)
+            {
+              double tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+
+              if (tmp < tmp_min
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_min)))
+                {
+                  tmp_min = tmp;
+                  idx_arg.elem (ir, j) = ir;
+                }
+              else
+                idx_arg.elem (ir, j) = idx_arg.elem (ir - 1, j);
+              if (tmp_min != 0.)
+                {
+                  result.xdata (ii) = tmp_min;
+                  result.xridx (ii++) = ir;
+                }
+            }
+
+          result.xcidx (j+1) = ii;
+        }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, expected_visit, nr);
+      OCTAVE_LOCAL_BUFFER (int8_t, min_sgn, nr);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          expected_visit[i] = 0;
+          min_sgn[i] = 1;
+        }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+          {
+            if (min_sgn[ridx (i)] >= 0)
+              {
+                if (min_sgn[ridx (i)] > 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || data (i) < 0.)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          min_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (data (i) < 0.)
+                  {
+                    min_sgn[ridx (i)] = 1;
+                    nel += nc - j;
+                  }
+              }
+          }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        if (min_sgn[i] > 0)
+          nel += expected_visit[i];
+
+      result = SparseMatrix (nr, nc, nel);
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+
+      octave_idx_type i = cidx (0);
+      for (octave_idx_type ir = 0; ir < nr; ir++)
+        {
+          idx_arg.elem (ir, 0) = 0;
+          if (i < cidx (1) && ridx (i) == ir)
+            {
+              result.xdata (ii) = data (i++);
+              result.xridx (ii++) = ir;
+            }
+        }
+      result.xcidx (1) = ii;
+      octave_idx_type r = result.cidx (0);
+
+      for (octave_idx_type j = 1; j < nc; j++)
+        {
+          for (octave_idx_type ir = 0; ir < nr; ir++)
+            {
+              double tmp_min = (r < result.cidx (j) && result.ridx (r) == ir
+                                ? result.data (r++) : 0.);
+              double tmp = i < cidx (j+1) && ridx (i) == ir ? data (i++) : 0.;
+              if (tmp < tmp_min
+                  || (! octave::math::isnan (tmp)
+                      && octave::math::isnan (tmp_min)))
+                {
+                  idx_arg.elem (ir, j) = j;
+                  if (tmp != 0.)
+                    {
+                      result.xdata (ii) = tmp;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+              else
+                {
+                  idx_arg.elem (ir, j) = idx_arg.elem (ir, j - 1);
+                  if (tmp_min != 0.)
+                    {
+                      result.xdata (ii) = tmp_min;
+                      result.xridx (ii++) = ir;
+                    }
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+/*
+
+%!test
+%! x = kron (eye (4), sparse (1, 1, 1, 2^14, 2^14));
+%! assert (sum (sum (cummax (x, 1))), sparse (2^17 + 2^15));
+%! assert (sum (sum (cummax (x, 2))), sparse (2^17 + 2^15));
+%! assert (sum (sum (cummin (x, 1))), sparse (1));
+%! assert (sum (sum (cummin (x, 2))), sparse (1));
+
+*/
+
 RowVector
 SparseMatrix::row (octave_idx_type i) const
 {
@@ -7793,6 +8620,9 @@ min (double d, const SparseMatrix& m)
 SparseMatrix
 min (const SparseMatrix& m, double d)
 {
+  if (octave::math::isnan (d))
+    return m;
+
   return min (d, m);
 }
 
@@ -7943,6 +8773,9 @@ max (double d, const SparseMatrix& m)
 SparseMatrix
 max (const SparseMatrix& m, double d)
 {
+  if (octave::math::isnan (d))
+    return m;
+
   return max (d, m);
 }
 
diff --git a/liboctave/array/dSparse.h b/liboctave/array/dSparse.h
--- a/liboctave/array/dSparse.h
+++ b/liboctave/array/dSparse.h
@@ -110,6 +110,11 @@ public:
   SparseMatrix min (int dim = -1) const;
   SparseMatrix min (Array<octave_idx_type>& index, int dim = -1) const;
 
+  SparseMatrix cummax (int dim = -1) const;
+  SparseMatrix cummax (Array<octave_idx_type>& index, int dim = -1) const;
+  SparseMatrix cummin (int dim = -1) const;
+  SparseMatrix cummin (Array<octave_idx_type>& index, int dim = -1) const;
+
   // destructive insert/delete/reorder operations
 
   SparseMatrix& insert (const SparseMatrix& a, octave_idx_type r,
diff --git a/liboctave/array/fCMatrix.cc b/liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc
+++ b/liboctave/array/fCMatrix.cc
@@ -3663,6 +3663,9 @@ min (const FloatComplex& c, const FloatC
 FloatComplexMatrix
 min (const FloatComplexMatrix& m, const FloatComplex& c)
 {
+  if (octave::math::isnan (c))
+    return m;
+
   return min (c, m);
 }
 
@@ -3681,33 +3684,11 @@ min (const FloatComplexMatrix& a, const 
   FloatComplexMatrix result (nr, nc);
 
   for (octave_idx_type j = 0; j < nc; j++)
-    {
-      bool columns_are_real_only = true;
-      for (octave_idx_type i = 0; i < nr; i++)
-        {
-          octave_quit ();
-          if (octave::math::imag (a(i, j)) != 0.0 || octave::math::imag (b(i, j)) != 0.0)
-            {
-              columns_are_real_only = false;
-              break;
-            }
-        }
-
-      if (columns_are_real_only)
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            result(i, j) = octave::math::min (octave::math::real (a(i, j)),
-                                              octave::math::real (b(i, j)));
-        }
-      else
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            {
-              octave_quit ();
-              result(i, j) = octave::math::min (a(i, j), b(i, j));
-            }
-        }
-    }
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+        octave_quit ();
+        result(i, j) = octave::math::min (a(i, j), b(i, j));
+      }
 
   return result;
 }
@@ -3735,6 +3716,9 @@ max (const FloatComplex& c, const FloatC
 FloatComplexMatrix
 max (const FloatComplexMatrix& m, const FloatComplex& c)
 {
+  if (octave::math::isnan (c))
+    return m;
+
   return max (c, m);
 }
 
@@ -3753,36 +3737,11 @@ max (const FloatComplexMatrix& a, const 
   FloatComplexMatrix result (nr, nc);
 
   for (octave_idx_type j = 0; j < nc; j++)
-    {
-      bool columns_are_real_only = true;
-      for (octave_idx_type i = 0; i < nr; i++)
-        {
-          octave_quit ();
-          if (octave::math::imag (a(i, j)) != 0.0 || octave::math::imag (b(i, j)) != 0.0)
-            {
-              columns_are_real_only = false;
-              break;
-            }
-        }
-
-      if (columns_are_real_only)
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            {
-              octave_quit ();
-              result(i, j) = octave::math::max (octave::math::real (a(i, j)),
-                                                octave::math::real (b(i, j)));
-            }
-        }
-      else
-        {
-          for (octave_idx_type i = 0; i < nr; i++)
-            {
-              octave_quit ();
-              result(i, j) = octave::math::max (a(i, j), b(i, j));
-            }
-        }
-    }
+    for (octave_idx_type i = 0; i < nr; i++)
+      {
+        octave_quit ();
+        result(i, j) = octave::math::max (a(i, j), b(i, j));
+      }
 
   return result;
 }
diff --git a/liboctave/array/fMatrix.cc b/liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc
+++ b/liboctave/array/fMatrix.cc
@@ -3052,21 +3052,10 @@ min (float d, const FloatMatrix& m)
 FloatMatrix
 min (const FloatMatrix& m, float d)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (FloatMatrix);
-
-  FloatMatrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::min (m(i, j), d);
-      }
-
-  return result;
+  if (octave::math::isnan (d))
+    return m;
+
+  return min (d, m);
 }
 
 FloatMatrix
@@ -3116,21 +3105,10 @@ max (float d, const FloatMatrix& m)
 FloatMatrix
 max (const FloatMatrix& m, float d)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (FloatMatrix);
-
-  FloatMatrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::max (m(i, j), d);
-      }
-
-  return result;
+  if (octave::math::isnan (d))
+    return m;
+
+  return max (d, m);
 }
 
 FloatMatrix
diff --git a/liboctave/numeric/lo-mappers.h b/liboctave/numeric/lo-mappers.h
--- a/liboctave/numeric/lo-mappers.h
+++ b/liboctave/numeric/lo-mappers.h
@@ -424,25 +424,25 @@ namespace octave
     inline std::complex<double>
     min (const std::complex<double>& x, const std::complex<double>& y)
     {
-      return abs (x) <= abs (y) ? x : (isnan (x) ? x : y);
+      return isnan (y) ? x : (x <= y ? x : y);
     }
 
     inline std::complex<float>
     min (const std::complex<float>& x, const std::complex<float>& y)
     {
-      return abs (x) <= abs (y) ? x : (isnan (x) ? x : y);
+      return isnan (y) ? x : (x <= y ? x : y);
     }
 
     inline std::complex<double>
     max (const std::complex<double>& x, const std::complex<double>& y)
     {
-      return abs (x) >= abs (y) ? x : (isnan (x) ? x : y);
+      return isnan (y) ? x : (x >= y ? x : y);
     }
 
     inline std::complex<float>
     max (const std::complex<float>& x, const std::complex<float>& y)
     {
-      return abs (x) >= abs (y) ? x : (isnan (x) ? x : y);
+      return isnan (y) ? x : (x >= y ? x : y);
     }
 
     template <typename T>
diff --git a/liboctave/operators/mx-inlines.cc b/liboctave/operators/mx-inlines.cc
--- a/liboctave/operators/mx-inlines.cc
+++ b/liboctave/operators/mx-inlines.cc
@@ -392,15 +392,6 @@ mx_inline_xmax (size_t n, T *r, T x, con
     else                                                        \
       for (size_t i = 0; i < n; i++)                            \
         r[i] = (x[i] OP y) ? x[i] : y;                          \
-  }                                                             \
-  template <>                                                   \
-  inline void F<T> (size_t n, T *r, T x, const T *y) throw ()   \
-  {                                                             \
-    if (octave::math::isnan (x))                                \
-      std::memcpy (r, y, n * sizeof (T));                       \
-    else                                                        \
-      for (size_t i = 0; i < n; i++)                            \
-        r[i] = (y[i] OP x) ? y[i] : x;                          \
   }
 
 DEFMINMAXSPEC (double, mx_inline_xmin, <=)
@@ -1008,7 +999,10 @@ OP_MINMAX_FCN (mx_inline_max, >)
         for (octave_idx_type i = 0; i < m; i++)                         \
           {                                                             \
             if (octave::math::isnan (v[i]))                             \
-              nan = true;                                               \
+              {                                                         \
+                if (octave::math::isnan (r[i]))                         \
+                  nan = true;                                           \
+              }                                                         \
             else if (octave::math::isnan (r[i]) || v[i] OP r[i])        \
               r[i] = v[i];                                              \
           }                                                             \
@@ -1048,7 +1042,10 @@ OP_MINMAX_FCN (mx_inline_max, >)
         for (octave_idx_type i = 0; i < m; i++)                         \
           {                                                             \
             if (octave::math::isnan (v[i]))                             \
-              nan = true;                                               \
+              {                                                         \
+                if (octave::math::isnan (r[i]))                         \
+                  nan = true;                                           \
+              }                                                         \
             else if (octave::math::isnan (r[i]) || v[i] OP r[i])        \
               {                                                         \
                 r[i] = v[i];                                            \
@@ -1235,7 +1232,8 @@ OP_CUMMINMAX_FCN (mx_inline_cummax, >)
             if (octave::math::isnan (v[i]))                             \
               {                                                         \
                 r[i] = r0[i];                                           \
-                nan = true;                                             \
+                if (octave::math::isnan (r[i]))                         \
+                  nan = true;                                           \
               }                                                         \
             else if (octave::math::isnan (r0[i]) || v[i] OP r0[i])      \
               r[i] = v[i];                                              \
@@ -1292,7 +1290,8 @@ OP_CUMMINMAX_FCN (mx_inline_cummax, >)
               {                                                         \
                 r[i] = r0[i];                                           \
                 ri[i] = r0i[i];                                         \
-                nan = true;                                             \
+                if (octave::math::isnan (r[i]))                         \
+                  nan = true;                                           \
               }                                                         \
             else if (octave::math::isnan (r0[i]) || v[i] OP r0[i])      \
               {                                                         \
