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

* max.cc (do_minmax_red_op, do_minmax_red_op<charNDArray>,
do_minmax_red_op<boolNDArray>, do_minmax_bin_op,
do_minmax_bin_op<charNDArray>, do_cumminmax_red_op): New includenan
argument.
(do_minmax_body, do_cumminmax_body): Process "omitnan" and
"includenan" options and pass the flag to do_minmax_red_op,
do_minmax_bin_op or do_cumminmax_red_op.
(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.h, CMatrix.cc, fCMatrix.h, fCMatrix.cc (min, max): New
includenan argument. Handle matrix - NaN scalar comparisons correctly.
Remove columns_are_real_only to keep complex ordering.
* dMatrix.h, dMatrix.cc, fMatrix.h, fMatrix.cc (min, max): New
includenan argument.

* CNDArray.cc, CNDArray.h, dNDArray.cc, dNDArray.h, fCNDArray.cc,
fCNDArray.h, fNDArray.cc, fNDArray.h, intNDArray.cc, intNDArray.h
(ComplexNDArray::max, ComplexNDArray::min, ComplexNDArray::cummax,
ComplexNDArray::cummin, NDArray::max, NDArray::min, NDArray::cummax,
NDArray::cummin, FloatComplexNDArray::max, FloatComplexNDArray::min,
FloatComplexNDArray::cummax, FloatComplexNDArray::cummin,
FloatNDArray::max, FloatNDArray::min, FloatNDArray::cummax,
FloatNDArray::cummin, intNDArray<T>::max, intNDArray<T>::min,
intNDArray<T>::cummax, intNDArray<T>::cummin): New includenan
argument.

* CSparse.cc, CSparse.h (SparseComplexMatrix::max,
SparseComplexMatrix::min): New includenan argument. 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): New includenan argument. Handle matrix - NaN scalar
comparisons correctly.

* dSparse.cc, dSparse.h (SparseMatrix::max, SparseMatrix::min): New
includenan argument.
(SparseMatrix::cummax, SparseMatrix::cummin): New
methods and new test for them.
(min, max): New includenan argument. Handle matrix - NaN scalar
comparisons correctly.

* 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_xmin_inan, mx_inline_xmax_inan): New functions.
(mx_inline_min, mx_inline_max, mx_inline_cummin, mx_inline_cummax):
New includenan argument. Proceed to the faster code in OP_MINMAX_FCN2
and OP_CUMMINMAX_FCN2 as soon as r contains no NaNs.
(do_mx_minmax_op, do_mx_cumminmax_op): New includenan argument. Update
callback function type.

* mx-op-decl.h, mx-op-defs.h: New includenan argument.

diff --git a/libinterp/corefcn/max.cc b/libinterp/corefcn/max.cc
--- a/libinterp/corefcn/max.cc
+++ b/libinterp/corefcn/max.cc
@@ -44,7 +44,7 @@ along with Octave; see the file COPYING.
 template <typename ArrayType>
 static octave_value_list
 do_minmax_red_op (const octave_value& arg,
-                  int nargout, int dim, bool ismin)
+                  int nargout, int dim, bool includenan, bool ismin)
 {
   octave_value_list retval (nargout > 1 ? 2 : 1);
   ArrayType array = octave_value_extract<ArrayType> (arg);
@@ -52,17 +52,17 @@ do_minmax_red_op (const octave_value& ar
   if (nargout <= 1)
     {
       if (ismin)
-        retval(0) = array.min (dim);
+        retval(0) = array.min (dim, includenan);
       else
-        retval(0) = array.max (dim);
+        retval(0) = array.max (dim, includenan);
     }
   else
     {
       Array<octave_idx_type> idx;
       if (ismin)
-        retval(0) = array.min (idx, dim);
+        retval(0) = array.min (idx, dim, includenan);
       else
-        retval(0) = array.max (idx, dim);
+        retval(0) = array.max (idx, dim, includenan);
 
       retval(1) = octave_value (idx, true, true);
     }
@@ -80,9 +80,11 @@ do_minmax_red_op (const octave_value& ar
 
 template <>
 octave_value_list
-do_minmax_red_op<charNDArray> (const octave_value& arg,
-                               int nargout, int dim, bool ismin)
+do_minmax_red_op<charNDArray> (const octave_value& arg, int nargout, int dim,
+                               bool includenan, bool ismin)
 {
+  octave_unused_parameter (includenan);
+
   octave_value_list retval (nargout > 1 ? 2 : 1);
   charNDArray array = octave_value_extract<charNDArray> (arg);
 
@@ -110,8 +112,8 @@ do_minmax_red_op<charNDArray> (const oct
 // Specialization for bool arrays (dense or sparse).
 template <>
 octave_value_list
-do_minmax_red_op<boolNDArray> (const octave_value& arg,
-                               int nargout, int dim, bool ismin)
+do_minmax_red_op<boolNDArray> (const octave_value& arg, int nargout, int dim,
+                               bool includenan, bool ismin)
 {
   octave_value_list retval;
 
@@ -132,7 +134,8 @@ do_minmax_red_op<boolNDArray> (const oct
       else
         {
           // any/all don't have indexed versions, so do it via a conversion.
-          retval = do_minmax_red_op<int8NDArray> (arg, nargout, dim, ismin);
+          retval = do_minmax_red_op<int8NDArray> (arg, nargout, dim,
+                                                  includenan, ismin);
 
           retval(0) = retval(0).bool_array_value ();
         }
@@ -141,7 +144,8 @@ do_minmax_red_op<boolNDArray> (const oct
     {
       // Sparse: Don't use any/all trick, as full matrix could exceed memory.
       // Instead, convert to double.
-      retval = do_minmax_red_op<SparseMatrix> (arg, nargout, dim, ismin);
+      retval = do_minmax_red_op<SparseMatrix> (arg, nargout, dim, includenan,
+                                               ismin);
 
       retval(0) = retval(0).sparse_bool_matrix_value ();
     }
@@ -152,7 +156,7 @@ do_minmax_red_op<boolNDArray> (const oct
 template <typename ArrayType>
 static octave_value
 do_minmax_bin_op (const octave_value& argx, const octave_value& argy,
-                  bool ismin)
+                  bool includenan, bool ismin)
 {
   typedef typename ArrayType::element_type ScalarType;
 
@@ -164,9 +168,9 @@ do_minmax_bin_op (const octave_value& ar
       ArrayType y = octave_value_extract<ArrayType> (argy);
 
       if (ismin)
-        retval = min (x, y);
+        retval = min (x, y, includenan);
       else
-        retval = max (x, y);
+        retval = max (x, y, includenan);
     }
   else if (argy.is_scalar_type ())
     {
@@ -174,9 +178,9 @@ do_minmax_bin_op (const octave_value& ar
       ScalarType y = octave_value_extract<ScalarType> (argy);
 
       if (ismin)
-        retval = min (x, y);
+        retval = min (x, y, includenan);
       else
-        retval = max (x, y);
+        retval = max (x, y, includenan);
     }
   else
     {
@@ -184,9 +188,9 @@ do_minmax_bin_op (const octave_value& ar
       ArrayType y = octave_value_extract<ArrayType> (argy);
 
       if (ismin)
-        retval = min (x, y);
+        retval = min (x, y, includenan);
       else
-        retval = max (x, y);
+        retval = max (x, y, includenan);
     }
 
   return retval;
@@ -203,8 +207,11 @@ do_minmax_bin_op (const octave_value& ar
 template <>
 octave_value
 do_minmax_bin_op<charNDArray> (const octave_value& argx,
-                               const octave_value& argy, bool ismin)
+                               const octave_value& argy, bool includenan,
+                               bool ismin)
 {
+  octave_unused_parameter (includenan);
+
   octave_value retval;
 
   charNDArray x = octave_value_extract<charNDArray> (argx);
@@ -238,13 +245,29 @@ do_minmax_body (const octave_value_list&
 {
   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);
@@ -290,9 +313,10 @@ do_minmax_body (const octave_value_list&
               }
             else if (arg.is_sparse_type ())
               retval = do_minmax_red_op<SparseMatrix> (arg, nargout, dim,
-                                                       ismin);
+                                                       includenan, ismin);
             else
-              retval = do_minmax_red_op<NDArray> (arg, nargout, dim, ismin);
+              retval = do_minmax_red_op<NDArray> (arg, nargout, dim,
+                                                  includenan, ismin);
 
           }
           break;
@@ -301,29 +325,34 @@ do_minmax_body (const octave_value_list&
           {
             if (arg.is_sparse_type ())
               retval = do_minmax_red_op<SparseComplexMatrix> (arg, nargout, dim,
+                                                              includenan,
                                                               ismin);
             else
               retval = do_minmax_red_op<ComplexNDArray> (arg, nargout, dim,
+                                                         includenan,
                                                          ismin);
           }
           break;
 
         case btyp_float:
-          retval = do_minmax_red_op<FloatNDArray> (arg, nargout, dim, ismin);
+          retval = do_minmax_red_op<FloatNDArray> (arg, nargout, dim,
+                                                   includenan, ismin);
           break;
 
         case btyp_float_complex:
           retval = do_minmax_red_op<FloatComplexNDArray> (arg, nargout, dim,
-                                                          ismin);
+                                                          includenan, ismin);
           break;
 
         case btyp_char:
-          retval = do_minmax_red_op<charNDArray> (arg, nargout, dim, ismin);
+          retval = do_minmax_red_op<charNDArray> (arg, nargout, dim,
+                                                  includenan, ismin);
           break;
 
 #define MAKE_INT_BRANCH(X)                                                    \
         case btyp_ ## X:                                                      \
-          retval = do_minmax_red_op<X ## NDArray> (arg, nargout, dim, ismin); \
+          retval = do_minmax_red_op<X ## NDArray> (arg, nargout, dim,         \
+                                                   includenan, ismin);        \
           break;
 
         MAKE_INT_BRANCH (int8);
@@ -338,7 +367,8 @@ do_minmax_body (const octave_value_list&
 #undef MAKE_INT_BRANCH
 
         case btyp_bool:
-          retval = do_minmax_red_op<boolNDArray> (arg, nargout, dim, ismin);
+          retval = do_minmax_red_op<boolNDArray> (arg, nargout, dim,
+                                                  includenan, ismin);
           break;
 
         default:
@@ -367,9 +397,11 @@ do_minmax_body (const octave_value_list&
             if ((argx.is_sparse_type ()
                  && (argy.is_sparse_type () || argy.is_scalar_type ()))
                 || (argy.is_sparse_type () && argx.is_scalar_type ()))
-              retval = do_minmax_bin_op<SparseMatrix> (argx, argy, ismin);
+              retval = do_minmax_bin_op<SparseMatrix> (argx, argy, includenan,
+                                                       ismin);
             else
-              retval = do_minmax_bin_op<NDArray> (argx, argy, ismin);
+              retval = do_minmax_bin_op<NDArray> (argx, argy, includenan,
+                                                  ismin);
           }
           break;
 
@@ -379,27 +411,33 @@ do_minmax_body (const octave_value_list&
                  && (argy.is_sparse_type () || argy.is_scalar_type ()))
                 || (argy.is_sparse_type () && argx.is_scalar_type ()))
               retval = do_minmax_bin_op<SparseComplexMatrix> (argx, argy,
+                                                              includenan,
                                                               ismin);
             else
-              retval = do_minmax_bin_op<ComplexNDArray> (argx, argy, ismin);
+              retval = do_minmax_bin_op<ComplexNDArray> (argx, argy,
+                                                         includenan, ismin);
           }
           break;
 
         case btyp_float:
-          retval = do_minmax_bin_op<FloatNDArray> (argx, argy, ismin);
+          retval = do_minmax_bin_op<FloatNDArray> (argx, argy, includenan,
+                                                   ismin);
           break;
 
         case btyp_float_complex:
-          retval = do_minmax_bin_op<FloatComplexNDArray> (argx, argy, ismin);
+          retval = do_minmax_bin_op<FloatComplexNDArray> (argx, argy,
+                                                          includenan, ismin);
           break;
 
         case btyp_char:
-          retval = do_minmax_bin_op<charNDArray> (argx, argy, ismin);
+          retval = do_minmax_bin_op<charNDArray> (argx, argy, includenan,
+                                                  ismin);
           break;
 
-#define MAKE_INT_BRANCH(X)                                             \
-        case btyp_ ## X:                                               \
-          retval = do_minmax_bin_op<X ## NDArray> (argx, argy, ismin); \
+#define MAKE_INT_BRANCH(X)                                                 \
+        case btyp_ ## X:                                                   \
+          retval = do_minmax_bin_op<X ## NDArray> (argx, argy, includenan, \
+                                                   ismin);                 \
           break;
 
         MAKE_INT_BRANCH (int8);
@@ -415,7 +453,8 @@ do_minmax_body (const octave_value_list&
 
         // FIXME: This is what should happen when boolNDArray has max()
         // case btyp_bool:
-        //   retval = do_minmax_bin_op<boolNDArray> (argx, argy, ismin);
+        //   retval = do_minmax_bin_op<boolNDArray> (argx, argy, includenan,
+        //                                           ismin);
         //   break;
 
         default:
@@ -438,6 +477,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 +492,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 +582,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 +595,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 +646,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 +707,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 +719,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 +734,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 +829,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 +888,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 +954,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)
@@ -881,7 +964,7 @@ the first index of the maximum value(s).
 template <typename ArrayType>
 static octave_value_list
 do_cumminmax_red_op (const octave_value& arg,
-                     int nargout, int dim, bool ismin)
+                     int nargout, int dim, bool includenan, bool ismin)
 {
   octave_value_list retval (nargout > 1 ? 2 : 1);
   ArrayType array = octave_value_extract<ArrayType> (arg);
@@ -889,18 +972,18 @@ do_cumminmax_red_op (const octave_value&
   if (nargout <= 1)
     {
       if (ismin)
-        retval(0) = array.cummin (dim);
+        retval(0) = array.cummin (dim, includenan);
       else
-        retval(0) = array.cummax (dim);
+        retval(0) = array.cummax (dim, includenan);
     }
   else
     {
       retval.resize (2);
       Array<octave_idx_type> idx;
       if (ismin)
-        retval(0) = array.cummin (idx, dim);
+        retval(0) = array.cummin (idx, dim, includenan);
       else
-        retval(0) = array.cummax (idx, dim);
+        retval(0) = array.cummax (idx, dim, includenan);
 
       retval(1) = octave_value (idx, true, true);
     }
@@ -914,11 +997,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,26 +1031,37 @@ 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,
+                                                    includenan, ismin);
+      else
+        retval = do_cumminmax_red_op<NDArray> (arg, nargout, dim, includenan,
+                                               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,
+                                                           includenan, ismin);
+      else
+        retval = do_cumminmax_red_op<ComplexNDArray> (arg, nargout, dim,
+                                                      includenan, ismin);
       break;
 
     case btyp_float:
-      retval = do_cumminmax_red_op<FloatNDArray> (arg, nargout, dim, ismin);
+      retval = do_cumminmax_red_op<FloatNDArray> (arg, nargout, dim,
+                                                  includenan, ismin);
       break;
 
     case btyp_float_complex:
       retval = do_cumminmax_red_op<FloatComplexNDArray> (arg, nargout, dim,
-                                                         ismin);
+                                                         includenan, ismin);
       break;
 
 #define MAKE_INT_BRANCH(X)                                                   \
     case btyp_ ## X:                                                         \
-      retval = do_cumminmax_red_op<X ## NDArray> (arg, nargout, dim, ismin); \
+      retval = do_cumminmax_red_op<X ## NDArray> (arg, nargout, dim,         \
+                                                  includenan, ismin);        \
       break;
 
       MAKE_INT_BRANCH (int8);
@@ -970,7 +1078,7 @@ do_cumminmax_body (const octave_value_li
     case btyp_bool:
       {
         retval = do_cumminmax_red_op<int8NDArray> (arg, nargout, dim,
-                                                   ismin);
+                                                   includenan, ismin);
         if (retval.length () > 0)
           retval(0) = retval(0).bool_array_value ();
       }
@@ -988,6 +1096,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 +1122,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 +1147,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 +1171,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 +1197,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 +1222,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
@@ -3619,33 +3619,68 @@ operator * (const ComplexMatrix& a, cons
     return T (nr, nc);
 
 ComplexMatrix
-min (const Complex& c, const ComplexMatrix& m)
+min (const Complex& c, const ComplexMatrix& m, bool includenan)
 {
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
   EMPTY_RETURN_CHECK (ComplexMatrix);
 
+  if (octave::math::isnan (c) && includenan)
+    return ComplexMatrix (nr, nc, c);
+
+  ComplexMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = c <= m(i, j) ? c : m(i, j);
+          }
+    }
+  else
+    {
+      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 (c, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+ComplexMatrix
+min (const ComplexMatrix& m, const Complex& c, bool includenan)
+{
+  if (! octave::math::isnan (c))
+    return min (c, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (ComplexMatrix);
+
   ComplexMatrix 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 (c, m(i, j));
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : c;
       }
 
   return result;
 }
 
 ComplexMatrix
-min (const ComplexMatrix& m, const Complex& c)
-{
-  return min (c, m);
-}
-
-ComplexMatrix
-min (const ComplexMatrix& a, const ComplexMatrix& b)
+min (const ComplexMatrix& a, const ComplexMatrix& b, bool includenan)
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.columns ();
@@ -3658,66 +3693,93 @@ min (const ComplexMatrix& a, const Compl
 
   ComplexMatrix result (nr, nc);
 
-  for (octave_idx_type j = 0; j < nc; j++)
+  if (includenan)
     {
-      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 j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) <= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
     }
 
   return result;
 }
 
 ComplexMatrix
-max (const Complex& c, const ComplexMatrix& m)
+max (const Complex& c, const ComplexMatrix& m, bool includenan)
 {
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
   EMPTY_RETURN_CHECK (ComplexMatrix);
 
+  if (octave::math::isnan (c) && includenan)
+    return ComplexMatrix (nr, nc, c);
+
+  ComplexMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = c >= m(i, j) ? c : m(i, j);
+          }
+    }
+  else
+    {
+      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 (c, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+ComplexMatrix
+max (const ComplexMatrix& m, const Complex& c, bool includenan)
+{
+  if (! octave::math::isnan (c))
+    return max (c, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (ComplexMatrix);
+
   ComplexMatrix 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 (c, m(i, j));
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : c;
       }
 
   return result;
 }
 
 ComplexMatrix
-max (const ComplexMatrix& m, const Complex& c)
-{
-  return max (c, m);
-}
-
-ComplexMatrix
-max (const ComplexMatrix& a, const ComplexMatrix& b)
+max (const ComplexMatrix& a, const ComplexMatrix& b, bool includenan)
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.columns ();
@@ -3730,37 +3792,25 @@ max (const ComplexMatrix& a, const Compl
 
   ComplexMatrix result (nr, nc);
 
-  for (octave_idx_type j = 0; j < nc; j++)
+  if (includenan)
     {
-      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 j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) >= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
     }
 
   return result;
diff --git a/liboctave/array/CMatrix.h b/liboctave/array/CMatrix.h
--- a/liboctave/array/CMatrix.h
+++ b/liboctave/array/CMatrix.h
@@ -443,15 +443,21 @@ extern OCTAVE_API ComplexMatrix operator
 extern OCTAVE_API ComplexMatrix operator * (const ComplexMatrix&,
                                             const ComplexMatrix&);
 
-extern OCTAVE_API ComplexMatrix min (const Complex& c, const ComplexMatrix& m);
-extern OCTAVE_API ComplexMatrix min (const ComplexMatrix& m, const Complex& c);
+extern OCTAVE_API ComplexMatrix min (const Complex& c, const ComplexMatrix& m,
+                                     bool includenan = false);
+extern OCTAVE_API ComplexMatrix min (const ComplexMatrix& m, const Complex& c,
+                                     bool includenan = false);
 extern OCTAVE_API ComplexMatrix min (const ComplexMatrix& a,
-                                     const ComplexMatrix& b);
+                                     const ComplexMatrix& b,
+                                     bool includenan = false);
 
-extern OCTAVE_API ComplexMatrix max (const Complex& c, const ComplexMatrix& m);
-extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& m, const Complex& c);
+extern OCTAVE_API ComplexMatrix max (const Complex& c, const ComplexMatrix& m,
+                                     bool includenan = false);
+extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& m, const Complex& c,
+                                     bool includenan = false);
 extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& a,
-                                     const ComplexMatrix& b);
+                                     const ComplexMatrix& b,
+                                     bool includenan = false);
 
 extern OCTAVE_API ComplexMatrix linspace (const ComplexColumnVector& x1,
                                           const ComplexColumnVector& x2,
diff --git a/liboctave/array/CNDArray.cc b/liboctave/array/CNDArray.cc
--- a/liboctave/array/CNDArray.cc
+++ b/liboctave/array/CNDArray.cc
@@ -664,51 +664,61 @@ static const Complex Complex_NaN_result 
                                          octave::numeric_limits<double>::NaN ());
 
 ComplexNDArray
-ComplexNDArray::max (int dim) const
+ComplexNDArray::max (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<Complex> (*this, dim, mx_inline_max);
+  return do_mx_minmax_op<Complex> (*this, dim, mx_inline_max, includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::max (Array<octave_idx_type>& idx_arg, int dim) const
+ComplexNDArray::max (Array<octave_idx_type>& idx_arg,
+                     int dim, bool includenan) const
 {
-  return do_mx_minmax_op<Complex> (*this, idx_arg, dim, mx_inline_max);
+  return do_mx_minmax_op<Complex> (*this, idx_arg, dim, mx_inline_max,
+                                   includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::min (int dim) const
+ComplexNDArray::min (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<Complex> (*this, dim, mx_inline_min);
+  return do_mx_minmax_op<Complex> (*this, dim, mx_inline_min, includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::min (Array<octave_idx_type>& idx_arg, int dim) const
+ComplexNDArray::min (Array<octave_idx_type>& idx_arg,
+                     int dim, bool includenan) const
 {
-  return do_mx_minmax_op<Complex> (*this, idx_arg, dim, mx_inline_min);
+  return do_mx_minmax_op<Complex> (*this, idx_arg, dim, mx_inline_min,
+                                   includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::cummax (int dim) const
+ComplexNDArray::cummax (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<Complex> (*this, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<Complex> (*this, dim, mx_inline_cummax,
+                                      includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+ComplexNDArray::cummax (Array<octave_idx_type>& idx_arg,
+                        int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<Complex> (*this, idx_arg, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<Complex> (*this, idx_arg, dim, mx_inline_cummax,
+                                      includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::cummin (int dim) const
+ComplexNDArray::cummin (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<Complex> (*this, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<Complex> (*this, dim, mx_inline_cummin,
+                                      includenan);
 }
 
 ComplexNDArray
-ComplexNDArray::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+ComplexNDArray::cummin (Array<octave_idx_type>& idx_arg,
+                        int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<Complex> (*this, idx_arg, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<Complex> (*this, idx_arg, dim, mx_inline_cummin,
+                                      includenan);
 }
 
 NDArray
diff --git a/liboctave/array/CNDArray.h b/liboctave/array/CNDArray.h
--- a/liboctave/array/CNDArray.h
+++ b/liboctave/array/CNDArray.h
@@ -86,15 +86,19 @@ public:
   ComplexNDArray concat (const NDArray& rb,
                          const Array<octave_idx_type>& ra_idx);
 
-  ComplexNDArray max (int dim = -1) const;
-  ComplexNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
-  ComplexNDArray min (int dim = -1) const;
-  ComplexNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
+  ComplexNDArray max (int dim = -1, bool includenan = false) const;
+  ComplexNDArray max (Array<octave_idx_type>& index,
+                      int dim = -1, bool includenan = false) const;
+  ComplexNDArray min (int dim = -1, bool includenan = false) const;
+  ComplexNDArray min (Array<octave_idx_type>& index,
+                      int dim = -1, bool includenan = false) const;
 
-  ComplexNDArray cummax (int dim = -1) const;
-  ComplexNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
-  ComplexNDArray cummin (int dim = -1) const;
-  ComplexNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
+  ComplexNDArray cummax (int dim = -1, bool includenan = false) const;
+  ComplexNDArray cummax (Array<octave_idx_type>& index,
+                         int dim = -1, bool includenan = false) const;
+  ComplexNDArray cummin (int dim = -1, bool includenan = false) const;
+  ComplexNDArray cummin (Array<octave_idx_type>& index,
+                         int dim = -1, bool includenan = false) const;
 
   ComplexNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
diff --git a/liboctave/array/CSparse.cc b/liboctave/array/CSparse.cc
--- a/liboctave/array/CSparse.cc
+++ b/liboctave/array/CSparse.cc
@@ -184,14 +184,15 @@ static const Complex
                       octave::numeric_limits<double>::NaN ());
 
 SparseComplexMatrix
-SparseComplexMatrix::max (int dim) const
+SparseComplexMatrix::max (int dim, bool includenan) const
 {
   Array<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return max (dummy_idx, dim, includenan);
 }
 
 SparseComplexMatrix
-SparseComplexMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
+SparseComplexMatrix::max (Array<octave_idx_type>& idx_arg,
+                          int dim, bool includenan) const
 {
   SparseComplexMatrix result;
   dim_vector dv = dims ();
@@ -217,8 +218,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,30 +229,36 @@ 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++)
             {
               Complex tmp = data (i);
 
               if (octave::math::isnan (tmp))
-                continue;
-
-              double abs_tmp = std::abs (tmp);
-
-              if (octave::math::isnan (abs_max) || abs_tmp > abs_max)
+                {
+                  if (includenan)
+                    {
+                      idx_j = ridx (i);
+                      tmp_max = tmp;
+                      break;
+                    }
+                  else
+                    continue;
+                }
+
+              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 (! includenan && octave::math::isnan (tmp_max))
+            idx_arg.elem (j) = 0;
+          else
+            idx_arg.elem (j) = idx_j;
+          if (tmp_max != 0.)
             nel++;
         }
 
@@ -299,10 +305,17 @@ SparseComplexMatrix::max (Array<octave_i
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               Complex tmp = data (i);
-
-              if (octave::math::isnan (tmp))
+              Complex best = elem (ir, ix);
+
+              if (includenan)
+                {
+                  if (! octave::math::isnan (best)
+                      && (octave::math::isnan (tmp) || tmp > best))
+                    idx_arg.elem (ir) = j;
+                }
+              else 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;
             }
         }
@@ -341,14 +354,15 @@ SparseComplexMatrix::max (Array<octave_i
 }
 
 SparseComplexMatrix
-SparseComplexMatrix::min (int dim) const
+SparseComplexMatrix::min (int dim, bool includenan) const
 {
   Array<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return min (dummy_idx, dim, includenan);
 }
 
 SparseComplexMatrix
-SparseComplexMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
+SparseComplexMatrix::min (Array<octave_idx_type>& idx_arg,
+                          int dim, bool includenan) const
 {
   SparseComplexMatrix result;
   dim_vector dv = dims ();
@@ -374,8 +388,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,30 +399,36 @@ 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++)
             {
               Complex tmp = data (i);
 
               if (octave::math::isnan (tmp))
-                continue;
-
-              double abs_tmp = std::abs (tmp);
-
-              if (octave::math::isnan (abs_min) || abs_tmp < abs_min)
+                {
+                  if (includenan)
+                    {
+                      idx_j = ridx (i);
+                      tmp_min = tmp;
+                      break;
+                    }
+                  else
+                    continue;
+                }
+
+              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 (! includenan && octave::math::isnan (tmp_min))
+            idx_arg.elem (j) = 0;
+          else
+            idx_arg.elem (j) = idx_j;
+          if (tmp_min != 0.)
             nel++;
         }
 
@@ -456,10 +475,17 @@ SparseComplexMatrix::min (Array<octave_i
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               Complex tmp = data (i);
-
-              if (octave::math::isnan (tmp))
+              Complex best = elem (ir, ix);
+
+              if (includenan)
+                {
+                  if (! octave::math::isnan (best)
+                      && (octave::math::isnan (tmp) || tmp < best))
+                    idx_arg.elem (ir) = j;
+                }
+              else 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 +538,871 @@ SparseComplexMatrix::min (Array<octave_i
 
 */
 
+SparseComplexMatrix
+SparseComplexMatrix::cummax (int dim, bool includenan) 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 (includenan)
+            {
+              if (i < cidx (j+1))
+                nel += nr - ridx (i);
+            }
+          else
+            {
+              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) && ! includenan)
+                    {
+                      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
+                      || (includenan == octave::math::isnan (data (i))
+                          && includenan != 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) && ! includenan)
+                        {
+                          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))
+                                     || includenan);
+
+                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;
+              if (includenan)
+                tmp = r_val >= i_val ? r_val : (octave::math::isnan (r_val)
+                                                ? r_val : i_val);
+              else
+                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, bool includenan) 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 (includenan)
+            {
+              if (i < cidx (j+1))
+                nel += nr - ridx (i);
+            }
+          else
+            {
+              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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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))
+                                     || includenan);
+
+                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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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, bool includenan) 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)) && includenan)
+                    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)) && includenan)
+              {
+                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_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 (! octave::math::isnan (tmp_min) || ! includenan)
+                    {
+                      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
+                      || (includenan == octave::math::isnan (data (i))
+                          && includenan != 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 (! octave::math::isnan (tmp_min) || ! includenan)
+                        {
+                          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)
+              {
+                bool end_of_zeros = (octave::math::isnan (data (i))
+                                     && includenan);
+
+                if (min_sgn[ridx (i)] > 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || end_of_zeros)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          min_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (end_of_zeros)
+                  {
+                    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 = 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;
+              if (includenan)
+                tmp = r_val <= i_val ? r_val : (octave::math::isnan (r_val)
+                                                ? r_val : i_val);
+              else
+                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;
+}
+
+SparseComplexMatrix
+SparseComplexMatrix::cummin (Array<octave_idx_type>& idx_arg,
+                             int dim, bool includenan) 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)) && includenan)
+                    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)) && includenan)
+              {
+                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_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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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)
+              {
+                bool end_of_zeros = (octave::math::isnan (data (i))
+                                     && includenan);
+
+                if (min_sgn[ridx (i)] > 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || end_of_zeros)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          min_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (end_of_zeros)
+                  {
+                    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 = 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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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
 {
@@ -7668,7 +8559,7 @@ operator * (const SparseComplexMatrix& a
     return T (nr, nc);
 
 SparseComplexMatrix
-min (const Complex& c, const SparseComplexMatrix& m)
+min (const Complex& c, const SparseComplexMatrix& m, bool includenan)
 {
   SparseComplexMatrix result;
 
@@ -7678,7 +8569,43 @@ min (const Complex& c, const SparseCompl
   EMPTY_RETURN_CHECK (SparseComplexMatrix);
 
   if (abs (c) == 0.)
-    return SparseComplexMatrix (nr, nc);
+    {
+      if (includenan)
+        {
+          octave_idx_type nel = 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)))
+                nel++;
+          result = SparseComplexMatrix (nr, nc, nel);
+
+          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;
+            }
+        }
+      else
+        return SparseComplexMatrix (nr, nc);
+    }
+  else if (includenan)
+    {
+      if (octave::math::isnan (c))
+        return SparseComplexMatrix (nr, nc, 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 (c < m.data (i))
+            result.data (i) = c;
+    }
   else
     {
       result = SparseComplexMatrix (m);
@@ -7692,13 +8619,31 @@ min (const Complex& c, const SparseCompl
 }
 
 SparseComplexMatrix
-min (const SparseComplexMatrix& m, const Complex& c)
-{
-  return min (c, m);
+min (const SparseComplexMatrix& m, const Complex& c, bool includenan)
+{
+  if (! octave::math::isnan (c))
+    return min (c, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (SparseComplexMatrix);
+
+  SparseComplexMatrix result (nr, nc, c);
+  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 (m.ridx (i) + j * nr) = m.data (i);
+
+  return result;
 }
 
 SparseComplexMatrix
-min (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
+min (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+     bool includenan)
 {
   SparseComplexMatrix r;
 
@@ -7728,11 +8673,10 @@ min (const SparseComplexMatrix& a, const
               octave_quit ();
               if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
                 {
-                  Complex tmp = octave::math::min (a.data (ja), 0.);
-                  if (tmp != 0.)
+                  if (octave::math::isnan (a.data (ja)) && includenan)
                     {
                       r.ridx (jx) = a.ridx (ja);
-                      r.data (jx) = tmp;
+                      r.data (jx) = a.data (ja);
                       jx++;
                     }
                   ja++;
@@ -7741,11 +8685,10 @@ min (const SparseComplexMatrix& a, const
               else if ((! ja_lt_max)
                        || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
                 {
-                  Complex tmp = octave::math::min (0., b.data (jb));
-                  if (tmp != 0.)
+                  if (octave::math::isnan (b.data (jb)) && includenan)
                     {
                       r.ridx (jx) = b.ridx (jb);
-                      r.data (jx) = tmp;
+                      r.data (jx) = b.data (jb);
                       jx++;
                     }
                   jb++;
@@ -7753,7 +8696,13 @@ min (const SparseComplexMatrix& a, const
                 }
               else
                 {
-                  Complex tmp = octave::math::min (a.data (ja), b.data (jb));
+                  Complex tmp;
+                  if (includenan)
+                    tmp = a.data (ja) <= b.data (jb) ? a.data (ja)
+                          : (octave::math::isnan (a.data (ja)) ? a.data (ja)
+                             : b.data (jb));
+                  else
+                    tmp = octave::math::min (a.data (ja), b.data (jb));
                   if (tmp != 0.)
                     {
                       r.data (jx) = tmp;
@@ -7785,7 +8734,7 @@ min (const SparseComplexMatrix& a, const
 }
 
 SparseComplexMatrix
-max (const Complex& c, const SparseComplexMatrix& m)
+max (const Complex& c, const SparseComplexMatrix& m, bool includenan)
 {
   SparseComplexMatrix result;
 
@@ -7794,28 +8743,89 @@ max (const Complex& c, const SparseCompl
 
   EMPTY_RETURN_CHECK (SparseComplexMatrix);
 
+  if (octave::math::isnan (c))
+    {
+      if (includenan)
+        result = SparseComplexMatrix (nr, nc, c);
+      else
+        {
+          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;
+        }
+    }
   // Count the number of nonzero elements
-  if (octave::math::max (c, 0.) != 0.)
+  else if (octave::math::max (c, 0.) != 0.)
     {
       result = SparseComplexMatrix (nr, nc, c);
+      if (includenan)
+        {
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+              if (m.data (i) > c || octave::math::isnan (m.data (i)))
+                result.xdata (m.ridx (i) + j * nr) = m.data (i);
+        }
+      else
+        {
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+              result.xdata (m.ridx (i) + j * nr) =
+                octave::math::max (c, m.data (i));
+        }
+    }
+  else if (includenan)
+    result = SparseComplexMatrix (m);
+  else
+    {
+      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++)
-          result.xdata (m.ridx (i) + j * nr) = octave::math::max (c, m.data (i));
+        {
+          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 ();
     }
-  else
-    result = SparseComplexMatrix (m);
 
   return result;
 }
 
 SparseComplexMatrix
-max (const SparseComplexMatrix& m, const Complex& c)
-{
-  return max (c, m);
+max (const SparseComplexMatrix& m, const Complex& c, bool includenan)
+{
+  if (! octave::math::isnan (c))
+    return max (c, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (SparseComplexMatrix);
+
+  SparseComplexMatrix result (nr, nc, c);
+  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 (m.ridx (i) + j * nr) = m.data (i);
+
+  return result;
 }
 
 SparseComplexMatrix
-max (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
+max (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+     bool includenan)
 {
   SparseComplexMatrix r;
 
@@ -7845,7 +8855,11 @@ max (const SparseComplexMatrix& a, const
               octave_quit ();
               if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
                 {
-                  Complex tmp = octave::math::max (a.data (ja), 0.);
+                  Complex tmp;
+                  if (octave::math::isnan (a.data (ja)) && ! includenan)
+                    tmp = 0.;
+                  else
+                    tmp = a.data (ja);
                   if (tmp != 0.)
                     {
                       r.ridx (jx) = a.ridx (ja);
@@ -7858,7 +8872,11 @@ max (const SparseComplexMatrix& a, const
               else if ((! ja_lt_max)
                        || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
                 {
-                  Complex tmp = octave::math::max (0., b.data (jb));
+                  Complex tmp;
+                  if (octave::math::isnan (b.data (jb)) && ! includenan)
+                    tmp = 0.;
+                  else
+                    tmp = b.data (jb);
                   if (tmp != 0.)
                     {
                       r.ridx (jx) = b.ridx (jb);
@@ -7870,7 +8888,13 @@ max (const SparseComplexMatrix& a, const
                 }
               else
                 {
-                  Complex tmp = octave::math::max (a.data (ja), b.data (jb));
+                  Complex tmp;
+                  if (includenan)
+                    tmp = a.data (ja) >= b.data (jb) ? a.data (ja)
+                          : (octave::math::isnan (a.data (ja)) ? a.data (ja)
+                             : b.data (jb));
+                  else
+                    tmp = octave::math::max (a.data (ja), b.data (jb));
                   if (tmp != 0.)
                     {
                       r.data (jx) = tmp;
diff --git a/liboctave/array/CSparse.h b/liboctave/array/CSparse.h
--- a/liboctave/array/CSparse.h
+++ b/liboctave/array/CSparse.h
@@ -116,10 +116,19 @@ public:
 
   bool is_hermitian (void) const;
 
-  SparseComplexMatrix max (int dim = -1) const;
-  SparseComplexMatrix max (Array<octave_idx_type>& index, int dim = -1) const;
-  SparseComplexMatrix min (int dim = -1) const;
-  SparseComplexMatrix min (Array<octave_idx_type>& index, int dim = -1) const;
+  SparseComplexMatrix max (int dim = -1, bool includenan = false) const;
+  SparseComplexMatrix max (Array<octave_idx_type>& index,
+                           int dim = -1, bool includenan = false) const;
+  SparseComplexMatrix min (int dim = -1, bool includenan = false) const;
+  SparseComplexMatrix min (Array<octave_idx_type>& index,
+                           int dim = -1, bool includenan = false) const;
+
+  SparseComplexMatrix cummax (int dim = -1, bool includenan = false) const;
+  SparseComplexMatrix cummax (Array<octave_idx_type>& index,
+                              int dim = -1, bool includenan = false) const;
+  SparseComplexMatrix cummin (int dim = -1, bool includenan = false) const;
+  SparseComplexMatrix cummin (Array<octave_idx_type>& index,
+                              int dim = -1, bool includenan = false) const;
 
   SparseComplexMatrix& insert (const SparseComplexMatrix& a,
                                octave_idx_type r, octave_idx_type c);
@@ -527,18 +536,24 @@ extern OCTAVE_API SparseComplexMatrix op
                                                   const PermMatrix&);
 
 extern OCTAVE_API SparseComplexMatrix min (const Complex& c,
-                                           const SparseComplexMatrix& m);
+                                           const SparseComplexMatrix& m,
+                                           bool includenan = false);
 extern OCTAVE_API SparseComplexMatrix min (const SparseComplexMatrix& m,
-                                           const Complex& c);
+                                           const Complex& c,
+                                           bool includenan = false);
 extern OCTAVE_API SparseComplexMatrix min (const SparseComplexMatrix& a,
-                                           const SparseComplexMatrix& b);
+                                           const SparseComplexMatrix& b,
+                                           bool includenan = false);
 
 extern OCTAVE_API SparseComplexMatrix max (const Complex& c,
-                                           const SparseComplexMatrix& m);
+                                           const SparseComplexMatrix& m,
+                                           bool includenan = false);
 extern OCTAVE_API SparseComplexMatrix max (const SparseComplexMatrix& m,
-                                           const Complex& c);
+                                           const Complex& c,
+                                           bool includenan = false);
 extern OCTAVE_API SparseComplexMatrix max (const SparseComplexMatrix& a,
-                                           const SparseComplexMatrix& b);
+                                           const SparseComplexMatrix& b,
+                                           bool includenan = false);
 
 SPARSE_SMS_CMP_OP_DECLS (SparseComplexMatrix, Complex, OCTAVE_API)
 SPARSE_SMS_BOOL_OP_DECLS (SparseComplexMatrix, Complex, OCTAVE_API)
diff --git a/liboctave/array/dMatrix.cc b/liboctave/array/dMatrix.cc
--- a/liboctave/array/dMatrix.cc
+++ b/liboctave/array/dMatrix.cc
@@ -3026,28 +3026,148 @@ operator * (const Matrix& a, const Matri
     return T (nr, nc);
 
 Matrix
-min (double d, const Matrix& m)
+min (double d, const Matrix& m, bool includenan)
 {
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
   EMPTY_RETURN_CHECK (Matrix);
 
+  if (octave::math::isnan (d) && includenan)
+    return Matrix (nr, nc, d);
+
+  Matrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = d <= m(i, j) ? d : m(i, j);
+          }
+    }
+  else
+    {
+      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 (d, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+Matrix
+min (const Matrix& m, double d, bool includenan)
+{
+  if (! octave::math::isnan (d))
+    return min (d, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  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 (d, m(i, j));
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : d;
       }
 
   return result;
 }
 
 Matrix
-min (const Matrix& m, double d)
+min (const Matrix& a, const Matrix& b, bool includenan)
 {
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.columns ();
+
+  if (nr != b.rows () || nc != b.columns ())
+    (*current_liboctave_error_handler)
+      ("two-arg min requires same size arguments");
+
+  EMPTY_RETURN_CHECK (Matrix);
+
+  Matrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) <= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
+    }
+
+  return result;
+}
+
+Matrix
+max (double d, const Matrix& m, bool includenan)
+{
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (Matrix);
+
+  if (octave::math::isnan (d) && includenan)
+    return Matrix (nr, nc, d);
+
+  Matrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = d >= m(i, j) ? d : m(i, j);
+          }
+    }
+  else
+    {
+      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 (d, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+Matrix
+max (const Matrix& m, double d, bool includenan)
+{
+  if (! octave::math::isnan (d))
+    return max (d, m, includenan);
+
+  if (! includenan)
+    return m;
+
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
@@ -3059,96 +3179,46 @@ min (const Matrix& m, double d)
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = octave::math::min (m(i, j), d);
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : d;
       }
 
   return result;
 }
 
 Matrix
-min (const Matrix& a, const Matrix& b)
+max (const Matrix& a, const Matrix& b, bool includenan)
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.columns ();
 
   if (nr != b.rows () || nc != b.columns ())
     (*current_liboctave_error_handler)
-      ("two-arg min requires same size arguments");
+      ("two-arg max requires same size arguments");
 
   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 (a(i, j), b(i, j));
-      }
-
-  return result;
-}
-
-Matrix
-max (double d, const Matrix& m)
-{
-  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 (d, m(i, j));
-      }
-
-  return result;
-}
-
-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;
-}
-
-Matrix
-max (const Matrix& a, const Matrix& b)
-{
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.columns ();
-
-  if (nr != b.rows () || nc != b.columns ())
-    (*current_liboctave_error_handler)
-      ("two-arg max requires same size arguments");
-
-  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 (a(i, j), b(i, j));
-      }
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) >= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
+    }
 
   return result;
 }
diff --git a/liboctave/array/dMatrix.h b/liboctave/array/dMatrix.h
--- a/liboctave/array/dMatrix.h
+++ b/liboctave/array/dMatrix.h
@@ -375,13 +375,19 @@ extern OCTAVE_API Matrix xgemm (const Ma
 
 extern OCTAVE_API Matrix operator * (const Matrix& a, const Matrix& b);
 
-extern OCTAVE_API Matrix min (double d, const Matrix& m);
-extern OCTAVE_API Matrix min (const Matrix& m, double d);
-extern OCTAVE_API Matrix min (const Matrix& a, const Matrix& b);
+extern OCTAVE_API Matrix min (double d, const Matrix& m,
+                              bool includenan = false);
+extern OCTAVE_API Matrix min (const Matrix& m, double d,
+                              bool includenan = false);
+extern OCTAVE_API Matrix min (const Matrix& a, const Matrix& b,
+                              bool includenan = false);
 
-extern OCTAVE_API Matrix max (double d, const Matrix& m);
-extern OCTAVE_API Matrix max (const Matrix& m, double d);
-extern OCTAVE_API Matrix max (const Matrix& a, const Matrix& b);
+extern OCTAVE_API Matrix max (double d, const Matrix& m,
+                              bool includenan = false);
+extern OCTAVE_API Matrix max (const Matrix& m, double d,
+                              bool includenan = false);
+extern OCTAVE_API Matrix max (const Matrix& a, const Matrix& b,
+                              bool includenan = false);
 
 extern OCTAVE_API Matrix linspace (const ColumnVector& x1,
                                    const ColumnVector& x2,
diff --git a/liboctave/array/dNDArray.cc b/liboctave/array/dNDArray.cc
--- a/liboctave/array/dNDArray.cc
+++ b/liboctave/array/dNDArray.cc
@@ -677,51 +677,57 @@ NDArray::sumsq (int dim) const
 }
 
 NDArray
-NDArray::max (int dim) const
+NDArray::max (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<double> (*this, dim, mx_inline_max);
+  return do_mx_minmax_op<double> (*this, dim, mx_inline_max, includenan);
 }
 
 NDArray
-NDArray::max (Array<octave_idx_type>& idx_arg, int dim) const
+NDArray::max (Array<octave_idx_type>& idx_arg, int dim, bool includenan) const
 {
-  return do_mx_minmax_op<double> (*this, idx_arg, dim, mx_inline_max);
+  return do_mx_minmax_op<double> (*this, idx_arg, dim, mx_inline_max,
+                                  includenan);
 }
 
 NDArray
-NDArray::min (int dim) const
+NDArray::min (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<double> (*this, dim, mx_inline_min);
+  return do_mx_minmax_op<double> (*this, dim, mx_inline_min, includenan);
 }
 
 NDArray
-NDArray::min (Array<octave_idx_type>& idx_arg, int dim) const
+NDArray::min (Array<octave_idx_type>& idx_arg, int dim, bool includenan) const
 {
-  return do_mx_minmax_op<double> (*this, idx_arg, dim, mx_inline_min);
+  return do_mx_minmax_op<double> (*this, idx_arg, dim, mx_inline_min,
+                                  includenan);
 }
 
 NDArray
-NDArray::cummax (int dim) const
+NDArray::cummax (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<double> (*this, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<double> (*this, dim, mx_inline_cummax, includenan);
 }
 
 NDArray
-NDArray::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+NDArray::cummax (Array<octave_idx_type>& idx_arg,
+                 int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<double> (*this, idx_arg, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<double> (*this, idx_arg, dim, mx_inline_cummax,
+                                     includenan);
 }
 
 NDArray
-NDArray::cummin (int dim) const
+NDArray::cummin (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<double> (*this, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<double> (*this, dim, mx_inline_cummin, includenan);
 }
 
 NDArray
-NDArray::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+NDArray::cummin (Array<octave_idx_type>& idx_arg,
+                 int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<double> (*this, idx_arg, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<double> (*this, idx_arg, dim, mx_inline_cummin,
+                                     includenan);
 }
 
 NDArray
diff --git a/liboctave/array/dNDArray.h b/liboctave/array/dNDArray.h
--- a/liboctave/array/dNDArray.h
+++ b/liboctave/array/dNDArray.h
@@ -105,15 +105,19 @@ public:
   charNDArray concat (const charNDArray& rb,
                       const Array<octave_idx_type>& ra_idx);
 
-  NDArray max (int dim = -1) const;
-  NDArray max (Array<octave_idx_type>& index, int dim = -1) const;
-  NDArray min (int dim = -1) const;
-  NDArray min (Array<octave_idx_type>& index, int dim = -1) const;
+  NDArray max (int dim = -1, bool includenan = false) const;
+  NDArray max (Array<octave_idx_type>& index,
+               int dim = -1, bool includenan = false) const;
+  NDArray min (int dim = -1, bool includenan = false) const;
+  NDArray min (Array<octave_idx_type>& index,
+               int dim = -1, bool includenan = false) const;
 
-  NDArray cummax (int dim = -1) const;
-  NDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
-  NDArray cummin (int dim = -1) const;
-  NDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
+  NDArray cummax (int dim = -1, bool includenan = false) const;
+  NDArray cummax (Array<octave_idx_type>& index,
+                  int dim = -1, bool includenan = false) const;
+  NDArray cummin (int dim = -1, bool includenan = false) const;
+  NDArray cummin (Array<octave_idx_type>& index,
+                  int dim = -1, bool includenan = false) const;
 
   NDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
diff --git a/liboctave/array/dSparse.cc b/liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc
+++ b/liboctave/array/dSparse.cc
@@ -186,14 +186,15 @@ SparseMatrix::insert (const SparseMatrix
 }
 
 SparseMatrix
-SparseMatrix::max (int dim) const
+SparseMatrix::max (int dim, bool includenan) const
 {
   Array<octave_idx_type> dummy_idx;
-  return max (dummy_idx, dim);
+  return max (dummy_idx, dim, includenan);
 }
 
 SparseMatrix
-SparseMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
+SparseMatrix::max (Array<octave_idx_type>& idx_arg,
+                   int dim, bool includenan) const
 {
   SparseMatrix result;
   dim_vector dv = dims ();
@@ -237,7 +238,16 @@ SparseMatrix::max (Array<octave_idx_type
               double tmp = data (i);
 
               if (octave::math::isnan (tmp))
-                continue;
+                {
+                  if (includenan)
+                    {
+                      idx_j = ridx (i);
+                      tmp_max = tmp;
+                      break;
+                    }
+                  else
+                    continue;
+                }
               else if (octave::math::isnan (tmp_max) || tmp > tmp_max)
                 {
                   idx_j = ridx (i);
@@ -246,7 +256,10 @@ SparseMatrix::max (Array<octave_idx_type
 
             }
 
-          idx_arg.elem (j) = octave::math::isnan (tmp_max) ? 0 : idx_j;
+          if (! includenan && octave::math::isnan (tmp_max))
+            idx_arg.elem (j) = 0;
+          else
+            idx_arg.elem (j) = idx_j;
           if (tmp_max != 0.)
             nel++;
         }
@@ -295,10 +308,17 @@ SparseMatrix::max (Array<octave_idx_type
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               double tmp = data (i);
-
-              if (octave::math::isnan (tmp))
+              double best = elem (ir, ix);
+
+              if (includenan)
+                {
+                  if (! octave::math::isnan (best)
+                      && (octave::math::isnan (tmp) || tmp > best))
+                    idx_arg.elem (ir) = j;
+                }
+              else 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;
             }
         }
@@ -337,14 +357,15 @@ SparseMatrix::max (Array<octave_idx_type
 }
 
 SparseMatrix
-SparseMatrix::min (int dim) const
+SparseMatrix::min (int dim, bool includenan) const
 {
   Array<octave_idx_type> dummy_idx;
-  return min (dummy_idx, dim);
+  return min (dummy_idx, dim, includenan);
 }
 
 SparseMatrix
-SparseMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
+SparseMatrix::min (Array<octave_idx_type>& idx_arg,
+                   int dim, bool includenan) const
 {
   SparseMatrix result;
   dim_vector dv = dims ();
@@ -388,7 +409,16 @@ SparseMatrix::min (Array<octave_idx_type
               double tmp = data (i);
 
               if (octave::math::isnan (tmp))
-                continue;
+                {
+                  if (includenan)
+                    {
+                      idx_j = ridx (i);
+                      tmp_min = tmp;
+                      break;
+                    }
+                  else
+                    continue;
+                }
               else if (octave::math::isnan (tmp_min) || tmp < tmp_min)
                 {
                   idx_j = ridx (i);
@@ -397,7 +427,10 @@ SparseMatrix::min (Array<octave_idx_type
 
             }
 
-          idx_arg.elem (j) = octave::math::isnan (tmp_min) ? 0 : idx_j;
+          if (! includenan && octave::math::isnan (tmp_min))
+            idx_arg.elem (j) = 0;
+          else
+            idx_arg.elem (j) = idx_j;
           if (tmp_min != 0.)
             nel++;
         }
@@ -446,10 +479,17 @@ SparseMatrix::min (Array<octave_idx_type
               octave_idx_type ir = ridx (i);
               octave_idx_type ix = idx_arg.elem (ir);
               double tmp = data (i);
-
-              if (octave::math::isnan (tmp))
+              double best = elem (ir, ix);
+
+              if (includenan)
+                {
+                  if (! octave::math::isnan (best)
+                      && (octave::math::isnan (tmp) || tmp < best))
+                    idx_arg.elem (ir) = j;
+                }
+              else 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 +542,871 @@ SparseMatrix::min (Array<octave_idx_type
 
 */
 
+SparseMatrix
+SparseMatrix::cummax (int dim, bool includenan) 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. || (octave::math::isnan (data (i))
+                                        && includenan))
+                    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. || (octave::math::isnan (data (i))
+                                  && includenan))
+              {
+                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) && ! includenan))
+                    {
+                      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
+                      || (includenan == octave::math::isnan (data (i))
+                          && includenan != 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) && ! includenan))
+                        {
+                          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 = (data (i) > 0.
+                                     || (octave::math::isnan (data (i))
+                                         && includenan));
+
+                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 = 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;
+              if (includenan)
+                tmp = r_val >= i_val ? r_val : (octave::math::isnan (r_val)
+                                                ? r_val : i_val);
+              else
+                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, bool includenan) 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. || (octave::math::isnan (data (i))
+                                        && includenan))
+                    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. || (octave::math::isnan (data (i))
+                                  && includenan))
+              {
+                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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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 = (data (i) > 0.
+                                     || (octave::math::isnan (data (i))
+                                         && includenan));
+
+                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 = 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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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, bool includenan) 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. || (octave::math::isnan (data (i))
+                                        && includenan))
+                    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. || (octave::math::isnan (data (i))
+                                  && includenan))
+              {
+                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) && ! includenan))
+                    {
+                      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
+                      || (includenan == octave::math::isnan (data (i))
+                          && includenan != 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) && ! includenan))
+                        {
+                          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)
+              {
+                bool end_of_zeros = (data (i) < 0.
+                                     || (octave::math::isnan (data (i))
+                                         && includenan));
+
+                if (min_sgn[ridx (i)] > 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || end_of_zeros)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          min_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (end_of_zeros)
+                  {
+                    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;
+              if (includenan)
+                tmp = r_val <= i_val ? r_val : (octave::math::isnan (r_val)
+                                                ? r_val : i_val);
+              else
+                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, bool includenan) 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. || (octave::math::isnan (data (i))
+                                        && includenan))
+                    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. || (octave::math::isnan (data (i))
+                                  && includenan))
+              {
+                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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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)
+              {
+                bool end_of_zeros = (data (i) < 0.
+                                     || (octave::math::isnan (data (i))
+                                         && includenan));
+
+                if (min_sgn[ridx (i)] > 0)
+                  {
+                    if (expected_visit[ridx (i)] < j || end_of_zeros)
+                      {
+                        nel += expected_visit[ridx (i)];
+                        if (expected_visit[ridx (i)] < j)
+                          min_sgn[ridx (i)] = 0;
+                      }
+                    else
+                      expected_visit[ridx (i)]++;
+                  }
+                if (end_of_zeros)
+                  {
+                    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
+                  || (includenan == octave::math::isnan (tmp)
+                      && includenan != 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
 {
@@ -7734,7 +8639,7 @@ operator * (const SparseMatrix& a, const
     return T (nr, nc);
 
 SparseMatrix
-min (double d, const SparseMatrix& m)
+min (double d, const SparseMatrix& m, bool includenan)
 {
   SparseMatrix result;
 
@@ -7750,7 +8655,11 @@ min (double d, const SparseMatrix& m)
       for (octave_idx_type j = 0; j < nc; j++)
         for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
           {
-            double tmp = octave::math::min (d, m.data (i));
+            double tmp = m.data (i);
+            if (! includenan)
+              tmp = octave::math::min (d, tmp);
+            else if (d < tmp)
+              tmp = d;
             if (tmp != 0.)
               {
                 octave_idx_type idx = m.ridx (i) + j * nr;
@@ -7761,11 +8670,21 @@ min (double d, const SparseMatrix& m)
     }
   else
     {
-      octave_idx_type nel = 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::min (d, m.data (i)) != 0.)
-            nel++;
+      if (octave::math::isnan (d) && includenan)
+        return SparseMatrix (nr, nc, d);
+
+      octave_idx_type nel;
+      if (d > 0. || octave::math::isnan (d))
+        nel = m.nnz ();
+      else
+        {
+          nel = 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 (m.data (i) < 0.
+                  || (includenan && octave::math::isnan (m.data (i))))
+                nel++;
+        }
 
       result = SparseMatrix (nr, nc, nel);
 
@@ -7775,7 +8694,11 @@ min (double d, const SparseMatrix& m)
         {
           for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
             {
-              double tmp = octave::math::min (d, m.data (i));
+              double tmp = m.data (i);
+              if (! includenan)
+                tmp = octave::math::min (d, tmp);
+              else if (d < tmp)
+                tmp = d;
 
               if (tmp != 0.)
                 {
@@ -7791,13 +8714,30 @@ min (double d, const SparseMatrix& m)
 }
 
 SparseMatrix
-min (const SparseMatrix& m, double d)
-{
-  return min (d, m);
+min (const SparseMatrix& m, double d, bool includenan)
+{
+  if (! octave::math::isnan (d))
+    return min (d, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (SparseMatrix);
+
+  SparseMatrix result (nr, nc, d);
+  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 (m.ridx (i) + j * nr) = m.data (i);
+
+  return result;
 }
 
 SparseMatrix
-min (const SparseMatrix& a, const SparseMatrix& b)
+min (const SparseMatrix& a, const SparseMatrix& b, bool includenan)
 {
   SparseMatrix r;
 
@@ -7827,7 +8767,11 @@ min (const SparseMatrix& a, const Sparse
               octave_quit ();
               if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
                 {
-                  double tmp = octave::math::min (a.data (ja), 0.);
+                  double tmp = a.data (ja);
+                  if (! includenan)
+                    tmp = octave::math::min (tmp, 0.);
+                  else if (0. < tmp)
+                    tmp = 0.;
                   if (tmp != 0.)
                     {
                       r.ridx (jx) = a.ridx (ja);
@@ -7840,7 +8784,11 @@ min (const SparseMatrix& a, const Sparse
               else if ((! ja_lt_max)
                        || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
                 {
-                  double tmp = octave::math::min (0., b.data (jb));
+                  double tmp = b.data (jb);
+                  if (! includenan)
+                    tmp = octave::math::min (0., tmp);
+                  else if (0. < tmp)
+                    tmp = 0.;
                   if (tmp != 0.)
                     {
                       r.ridx (jx) = b.ridx (jb);
@@ -7852,7 +8800,13 @@ min (const SparseMatrix& a, const Sparse
                 }
               else
                 {
-                  double tmp = octave::math::min (a.data (ja), b.data (jb));
+                  double da = a.data (ja);
+                  double db = b.data (jb);
+                  double tmp;
+                  if (includenan)
+                    tmp = da <= db ? da : (octave::math::isnan (da) ? da : db);
+                  else
+                    tmp = octave::math::min (da, db);
                   if (tmp != 0.)
                     {
                       r.data (jx) = tmp;
@@ -7884,7 +8838,7 @@ min (const SparseMatrix& a, const Sparse
 }
 
 SparseMatrix
-max (double d, const SparseMatrix& m)
+max (double d, const SparseMatrix& m, bool includenan)
 {
   SparseMatrix result;
 
@@ -7900,7 +8854,11 @@ max (double d, const SparseMatrix& m)
       for (octave_idx_type j = 0; j < nc; j++)
         for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
           {
-            double tmp = octave::math::max (d, m.data (i));
+            double tmp = m.data (i);
+            if (! includenan)
+              tmp = octave::math::max (d, tmp);
+            else if (d > tmp)
+              tmp = d;
 
             if (tmp != 0.)
               {
@@ -7912,11 +8870,21 @@ max (double d, const SparseMatrix& m)
     }
   else
     {
-      octave_idx_type nel = 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::max (d, m.data (i)) != 0.)
-            nel++;
+      if (octave::math::isnan (d) && includenan)
+        return SparseMatrix (nr, nc, d);
+
+      octave_idx_type nel;
+      if (d < 0. || octave::math::isnan (d))
+        nel = m.nnz ();
+      else
+        {
+          nel = 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 (m.data (i) > 0.
+                  || (includenan && octave::math::isnan (m.data (i))))
+                nel++;
+        }
 
       result = SparseMatrix (nr, nc, nel);
 
@@ -7926,7 +8894,11 @@ max (double d, const SparseMatrix& m)
         {
           for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
             {
-              double tmp = octave::math::max (d, m.data (i));
+              double tmp = m.data (i);
+              if (! includenan)
+                tmp = octave::math::max (d, tmp);
+              else if (d > tmp)
+                tmp = d;
               if (tmp != 0.)
                 {
                   result.xdata (ii) = tmp;
@@ -7941,13 +8913,30 @@ max (double d, const SparseMatrix& m)
 }
 
 SparseMatrix
-max (const SparseMatrix& m, double d)
-{
-  return max (d, m);
+max (const SparseMatrix& m, double d, bool includenan)
+{
+  if (! octave::math::isnan (d))
+    return max (d, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (SparseMatrix);
+
+  SparseMatrix result (nr, nc, d);
+  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 (m.ridx (i) + j * nr) = m.data (i);
+
+  return result;
 }
 
 SparseMatrix
-max (const SparseMatrix& a, const SparseMatrix& b)
+max (const SparseMatrix& a, const SparseMatrix& b, bool includenan)
 {
   SparseMatrix r;
 
@@ -7977,7 +8966,11 @@ max (const SparseMatrix& a, const Sparse
               octave_quit ();
               if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
                 {
-                  double tmp = octave::math::max (a.data (ja), 0.);
+                  double tmp = a.data (ja);
+                  if (! includenan)
+                    tmp = octave::math::max (tmp, 0.);
+                  else if (0. > tmp)
+                    tmp = 0.;
                   if (tmp != 0.)
                     {
                       r.ridx (jx) = a.ridx (ja);
@@ -7990,7 +8983,11 @@ max (const SparseMatrix& a, const Sparse
               else if ((! ja_lt_max)
                        || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
                 {
-                  double tmp = octave::math::max (0., b.data (jb));
+                  double tmp = b.data (jb);
+                  if (! includenan)
+                    tmp = octave::math::max (0., tmp);
+                  else if (0. > tmp)
+                    tmp = 0.;
                   if (tmp != 0.)
                     {
                       r.ridx (jx) = b.ridx (jb);
@@ -8002,7 +8999,13 @@ max (const SparseMatrix& a, const Sparse
                 }
               else
                 {
-                  double tmp = octave::math::max (a.data (ja), b.data (jb));
+                  double da = a.data (ja);
+                  double db = b.data (jb);
+                  double tmp;
+                  if (includenan)
+                    tmp = da >= db ? da : (octave::math::isnan (da) ? da : db);
+                  else
+                    tmp = octave::math::max (da, db);
                   if (tmp != 0.)
                     {
                       r.data (jx) = tmp;
diff --git a/liboctave/array/dSparse.h b/liboctave/array/dSparse.h
--- a/liboctave/array/dSparse.h
+++ b/liboctave/array/dSparse.h
@@ -105,10 +105,19 @@ public:
 
   bool is_symmetric (void) const;
 
-  SparseMatrix max (int dim = -1) const;
-  SparseMatrix max (Array<octave_idx_type>& index, int dim = -1) const;
-  SparseMatrix min (int dim = -1) const;
-  SparseMatrix min (Array<octave_idx_type>& index, int dim = -1) const;
+  SparseMatrix max (int dim = -1, bool includenan = false) const;
+  SparseMatrix max (Array<octave_idx_type>& index,
+                    int dim = -1, bool includenan = false) const;
+  SparseMatrix min (int dim = -1, bool includenan = false) const;
+  SparseMatrix min (Array<octave_idx_type>& index,
+                    int dim = -1, bool includenan = false) const;
+
+  SparseMatrix cummax (int dim = -1, bool includenan = false) const;
+  SparseMatrix cummax (Array<octave_idx_type>& index,
+                       int dim = -1, bool includenan = false) const;
+  SparseMatrix cummin (int dim = -1, bool includenan = false) const;
+  SparseMatrix cummin (Array<octave_idx_type>& index,
+                       int dim = -1, bool includenan = false) const;
 
   // destructive insert/delete/reorder operations
 
@@ -488,15 +497,21 @@ extern OCTAVE_API SparseMatrix operator 
 extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&,
                                            const PermMatrix&);
 
-extern OCTAVE_API SparseMatrix min (double d, const SparseMatrix& m);
-extern OCTAVE_API SparseMatrix min (const SparseMatrix& m, double d);
+extern OCTAVE_API SparseMatrix min (double d, const SparseMatrix& m,
+                                    bool includenan = false);
+extern OCTAVE_API SparseMatrix min (const SparseMatrix& m, double d,
+                                    bool includenan = false);
 extern OCTAVE_API SparseMatrix min (const SparseMatrix& a,
-                                    const SparseMatrix& b);
+                                    const SparseMatrix& b,
+                                    bool includenan = false);
 
-extern OCTAVE_API SparseMatrix max (double d, const SparseMatrix& m);
-extern OCTAVE_API SparseMatrix max (const SparseMatrix& m, double d);
+extern OCTAVE_API SparseMatrix max (double d, const SparseMatrix& m,
+                                    bool includenan = false);
+extern OCTAVE_API SparseMatrix max (const SparseMatrix& m, double d,
+                                    bool includenan = false);
 extern OCTAVE_API SparseMatrix max (const SparseMatrix& a,
-                                    const SparseMatrix& b);
+                                    const SparseMatrix& b,
+                                    bool includenan = false);
 
 SPARSE_SMS_CMP_OP_DECLS (SparseMatrix, double, OCTAVE_API)
 SPARSE_SMS_BOOL_OP_DECLS (SparseMatrix, double, OCTAVE_API)
diff --git a/liboctave/array/fCMatrix.cc b/liboctave/array/fCMatrix.cc
--- a/liboctave/array/fCMatrix.cc
+++ b/liboctave/array/fCMatrix.cc
@@ -3641,33 +3641,68 @@ operator * (const FloatComplexMatrix& a,
     return T (nr, nc);
 
 FloatComplexMatrix
-min (const FloatComplex& c, const FloatComplexMatrix& m)
+min (const FloatComplex& c, const FloatComplexMatrix& m, bool includenan)
 {
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
   EMPTY_RETURN_CHECK (FloatComplexMatrix);
 
+  if (octave::math::isnan (c) && includenan)
+    return FloatComplexMatrix (nr, nc, c);
+
+  FloatComplexMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = c <= m(i, j) ? c : m(i, j);
+          }
+    }
+  else
+    {
+      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 (c, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+FloatComplexMatrix
+min (const FloatComplexMatrix& m, const FloatComplex& c, bool includenan)
+{
+  if (! octave::math::isnan (c))
+    return min (c, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (FloatComplexMatrix);
+
   FloatComplexMatrix 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 (c, m(i, j));
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : c;
       }
 
   return result;
 }
 
 FloatComplexMatrix
-min (const FloatComplexMatrix& m, const FloatComplex& c)
-{
-  return min (c, m);
-}
-
-FloatComplexMatrix
-min (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
+min (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool includenan)
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.columns ();
@@ -3680,66 +3715,93 @@ min (const FloatComplexMatrix& a, const 
 
   FloatComplexMatrix result (nr, nc);
 
-  for (octave_idx_type j = 0; j < nc; j++)
+  if (includenan)
     {
-      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 j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) <= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
     }
 
   return result;
 }
 
 FloatComplexMatrix
-max (const FloatComplex& c, const FloatComplexMatrix& m)
+max (const FloatComplex& c, const FloatComplexMatrix& m, bool includenan)
 {
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
   EMPTY_RETURN_CHECK (FloatComplexMatrix);
 
+  if (octave::math::isnan (c) && includenan)
+    return FloatComplexMatrix (nr, nc, c);
+
+  FloatComplexMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = c >= m(i, j) ? c : m(i, j);
+          }
+    }
+  else
+    {
+      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 (c, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+FloatComplexMatrix
+max (const FloatComplexMatrix& m, const FloatComplex& c, bool includenan)
+{
+  if (! octave::math::isnan (c))
+    return max (c, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (FloatComplexMatrix);
+
   FloatComplexMatrix 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 (c, m(i, j));
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : c;
       }
 
   return result;
 }
 
 FloatComplexMatrix
-max (const FloatComplexMatrix& m, const FloatComplex& c)
-{
-  return max (c, m);
-}
-
-FloatComplexMatrix
-max (const FloatComplexMatrix& a, const FloatComplexMatrix& b)
+max (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool includenan)
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.columns ();
@@ -3752,36 +3814,25 @@ max (const FloatComplexMatrix& a, const 
 
   FloatComplexMatrix result (nr, nc);
 
-  for (octave_idx_type j = 0; j < nc; j++)
+  if (includenan)
     {
-      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 j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) >= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
     }
 
   return result;
diff --git a/liboctave/array/fCMatrix.h b/liboctave/array/fCMatrix.h
--- a/liboctave/array/fCMatrix.h
+++ b/liboctave/array/fCMatrix.h
@@ -467,18 +467,24 @@ extern OCTAVE_API FloatComplexMatrix ope
                                                  const FloatComplexMatrix&);
 
 extern OCTAVE_API FloatComplexMatrix min (const FloatComplex& c,
-                                          const FloatComplexMatrix& m);
+                                          const FloatComplexMatrix& m,
+                                          bool includenan = false);
 extern OCTAVE_API FloatComplexMatrix min (const FloatComplexMatrix& m,
-                                          const FloatComplex& c);
+                                          const FloatComplex& c,
+                                          bool includenan = false);
 extern OCTAVE_API FloatComplexMatrix min (const FloatComplexMatrix& a,
-                                          const FloatComplexMatrix& b);
+                                          const FloatComplexMatrix& b,
+                                          bool includenan = false);
 
 extern OCTAVE_API FloatComplexMatrix max (const FloatComplex& c,
-                                          const FloatComplexMatrix& m);
+                                          const FloatComplexMatrix& m,
+                                          bool includenan = false);
 extern OCTAVE_API FloatComplexMatrix max (const FloatComplexMatrix& m,
-                                          const FloatComplex& c);
+                                          const FloatComplex& c,
+                                          bool includenan = false);
 extern OCTAVE_API FloatComplexMatrix max (const FloatComplexMatrix& a,
-                                          const FloatComplexMatrix& b);
+                                          const FloatComplexMatrix& b,
+                                          bool includenan = false);
 
 extern OCTAVE_API FloatComplexMatrix
 linspace (const FloatComplexColumnVector& x1,
diff --git a/liboctave/array/fCNDArray.cc b/liboctave/array/fCNDArray.cc
--- a/liboctave/array/fCNDArray.cc
+++ b/liboctave/array/fCNDArray.cc
@@ -668,53 +668,61 @@ static const FloatComplex FloatComplex_N
                                                    octave::numeric_limits<float>::NaN ());
 
 FloatComplexNDArray
-FloatComplexNDArray::max (int dim) const
+FloatComplexNDArray::max (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<FloatComplex> (*this, dim, mx_inline_max);
+  return do_mx_minmax_op<FloatComplex> (*this, dim, mx_inline_max, includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::max (Array<octave_idx_type>& idx_arg, int dim) const
+FloatComplexNDArray::max (Array<octave_idx_type>& idx_arg,
+                          int dim, bool includenan) const
 {
-  return do_mx_minmax_op<FloatComplex> (*this, idx_arg, dim, mx_inline_max);
+  return do_mx_minmax_op<FloatComplex> (*this, idx_arg, dim, mx_inline_max,
+                                        includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::min (int dim) const
+FloatComplexNDArray::min (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<FloatComplex> (*this, dim, mx_inline_min);
+  return do_mx_minmax_op<FloatComplex> (*this, dim, mx_inline_min, includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::min (Array<octave_idx_type>& idx_arg, int dim) const
+FloatComplexNDArray::min (Array<octave_idx_type>& idx_arg,
+                          int dim, bool includenan) const
 {
-  return do_mx_minmax_op<FloatComplex> (*this, idx_arg, dim, mx_inline_min);
+  return do_mx_minmax_op<FloatComplex> (*this, idx_arg, dim, mx_inline_min,
+                                        includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::cummax (int dim) const
+FloatComplexNDArray::cummax (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<FloatComplex> (*this, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<FloatComplex> (*this, dim, mx_inline_cummax,
+                                           includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+FloatComplexNDArray::cummax (Array<octave_idx_type>& idx_arg,
+                             int dim, bool includenan) const
 {
   return do_mx_cumminmax_op<FloatComplex> (*this, idx_arg, dim,
-                                           mx_inline_cummax);
+                                           mx_inline_cummax, includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::cummin (int dim) const
+FloatComplexNDArray::cummin (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<FloatComplex> (*this, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<FloatComplex> (*this, dim, mx_inline_cummin,
+                                           includenan);
 }
 
 FloatComplexNDArray
-FloatComplexNDArray::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+FloatComplexNDArray::cummin (Array<octave_idx_type>& idx_arg,
+                             int dim, bool includenan) const
 {
   return do_mx_cumminmax_op<FloatComplex> (*this, idx_arg, dim,
-                                           mx_inline_cummin);
+                                           mx_inline_cummin, includenan);
 }
 
 FloatNDArray
diff --git a/liboctave/array/fCNDArray.h b/liboctave/array/fCNDArray.h
--- a/liboctave/array/fCNDArray.h
+++ b/liboctave/array/fCNDArray.h
@@ -88,17 +88,19 @@ public:
   FloatComplexNDArray concat (const FloatNDArray& rb,
                               const Array<octave_idx_type>& ra_idx);
 
-  FloatComplexNDArray max (int dim = -1) const;
-  FloatComplexNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
-  FloatComplexNDArray min (int dim = -1) const;
-  FloatComplexNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatComplexNDArray max (int dim = -1, bool includenan = false) const;
+  FloatComplexNDArray max (Array<octave_idx_type>& index,
+                           int dim = -1, bool includenan = false) const;
+  FloatComplexNDArray min (int dim = -1, bool includenan = false) const;
+  FloatComplexNDArray min (Array<octave_idx_type>& index,
+                           int dim = -1, bool includenan = false) const;
 
-  FloatComplexNDArray cummax (int dim = -1) const;
+  FloatComplexNDArray cummax (int dim = -1, bool includenan = false) const;
   FloatComplexNDArray cummax (Array<octave_idx_type>& index,
-                              int dim = -1) const;
-  FloatComplexNDArray cummin (int dim = -1) const;
+                              int dim = -1, bool includenan = false) const;
+  FloatComplexNDArray cummin (int dim = -1, bool includenan = false) const;
   FloatComplexNDArray cummin (Array<octave_idx_type>& index,
-                              int dim = -1) const;
+                              int dim = -1, bool includenan = false) const;
 
   FloatComplexNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
diff --git a/liboctave/array/fMatrix.cc b/liboctave/array/fMatrix.cc
--- a/liboctave/array/fMatrix.cc
+++ b/liboctave/array/fMatrix.cc
@@ -3030,28 +3030,148 @@ operator * (const FloatMatrix& a, const 
     return T (nr, nc);
 
 FloatMatrix
-min (float d, const FloatMatrix& m)
+min (float d, const FloatMatrix& m, bool includenan)
 {
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
   EMPTY_RETURN_CHECK (FloatMatrix);
 
+  if (octave::math::isnan (d) && includenan)
+    return FloatMatrix (nr, nc, d);
+
+  FloatMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = d <= m(i, j) ? d : m(i, j);
+          }
+    }
+  else
+    {
+      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 (d, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+FloatMatrix
+min (const FloatMatrix& m, float d, bool includenan)
+{
+  if (! octave::math::isnan (d))
+    return min (d, m, includenan);
+
+  if (! includenan)
+    return m;
+
+  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 (d, m(i, j));
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : d;
       }
 
   return result;
 }
 
 FloatMatrix
-min (const FloatMatrix& m, float d)
+min (const FloatMatrix& a, const FloatMatrix& b, bool includenan)
 {
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.columns ();
+
+  if (nr != b.rows () || nc != b.columns ())
+    (*current_liboctave_error_handler)
+      ("two-arg min requires same size arguments");
+
+  EMPTY_RETURN_CHECK (FloatMatrix);
+
+  FloatMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) <= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
+    }
+
+  return result;
+}
+
+FloatMatrix
+max (float d, const FloatMatrix& m, bool includenan)
+{
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (FloatMatrix);
+
+  if (octave::math::isnan (d) && includenan)
+    return FloatMatrix (nr, nc, d);
+
+  FloatMatrix result (nr, nc);
+
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = d >= m(i, j) ? d : m(i, j);
+          }
+    }
+  else
+    {
+      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 (d, m(i, j));
+          }
+    }
+
+  return result;
+}
+
+FloatMatrix
+max (const FloatMatrix& m, float d, bool includenan)
+{
+  if (! octave::math::isnan (d))
+    return max (d, m, includenan);
+
+  if (! includenan)
+    return m;
+
   octave_idx_type nr = m.rows ();
   octave_idx_type nc = m.columns ();
 
@@ -3063,96 +3183,46 @@ min (const FloatMatrix& m, float d)
     for (octave_idx_type i = 0; i < nr; i++)
       {
         octave_quit ();
-        result(i, j) = octave::math::min (m(i, j), d);
+        result(i, j) = octave::math::isnan (m(i, j)) ? m(i, j) : d;
       }
 
   return result;
 }
 
 FloatMatrix
-min (const FloatMatrix& a, const FloatMatrix& b)
+max (const FloatMatrix& a, const FloatMatrix& b, bool includenan)
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.columns ();
 
   if (nr != b.rows () || nc != b.columns ())
     (*current_liboctave_error_handler)
-      ("two-arg min requires same size arguments");
+      ("two-arg max requires same size arguments");
 
   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 (a(i, j), b(i, j));
-      }
-
-  return result;
-}
-
-FloatMatrix
-max (float d, const FloatMatrix& m)
-{
-  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 (d, m(i, j));
-      }
-
-  return result;
-}
-
-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;
-}
-
-FloatMatrix
-max (const FloatMatrix& a, const FloatMatrix& b)
-{
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.columns ();
-
-  if (nr != b.rows () || nc != b.columns ())
-    (*current_liboctave_error_handler)
-      ("two-arg max requires same size arguments");
-
-  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 (a(i, j), b(i, j));
-      }
+  if (includenan)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            octave_quit ();
+            result(i, j) = a(i, j) >= b(i, j) ? a(i, j)
+                           : (octave::math::isnan (a(i, j)) ? a(i, j)
+                              : b(i, j));
+          }
+    }
+  else
+    {
+      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 (a(i, j), b(i, j));
+          }
+    }
 
   return result;
 }
diff --git a/liboctave/array/fMatrix.h b/liboctave/array/fMatrix.h
--- a/liboctave/array/fMatrix.h
+++ b/liboctave/array/fMatrix.h
@@ -391,13 +391,19 @@ extern OCTAVE_API FloatMatrix xgemm (con
 extern OCTAVE_API FloatMatrix operator * (const FloatMatrix& a,
                                           const FloatMatrix& b);
 
-extern OCTAVE_API FloatMatrix min (float d, const FloatMatrix& m);
-extern OCTAVE_API FloatMatrix min (const FloatMatrix& m, float d);
-extern OCTAVE_API FloatMatrix min (const FloatMatrix& a, const FloatMatrix& b);
+extern OCTAVE_API FloatMatrix min (float d, const FloatMatrix& m,
+                                   bool includenan = false);
+extern OCTAVE_API FloatMatrix min (const FloatMatrix& m, float d,
+                                   bool includenan = false);
+extern OCTAVE_API FloatMatrix min (const FloatMatrix& a, const FloatMatrix& b,
+                                   bool includenan = false);
 
-extern OCTAVE_API FloatMatrix max (float d, const FloatMatrix& m);
-extern OCTAVE_API FloatMatrix max (const FloatMatrix& m, float d);
-extern OCTAVE_API FloatMatrix max (const FloatMatrix& a, const FloatMatrix& b);
+extern OCTAVE_API FloatMatrix max (float d, const FloatMatrix& m,
+                                   bool includenan = false);
+extern OCTAVE_API FloatMatrix max (const FloatMatrix& m, float d,
+                                   bool includenan = false);
+extern OCTAVE_API FloatMatrix max (const FloatMatrix& a, const FloatMatrix& b,
+                                   bool includenan = false);
 
 extern OCTAVE_API FloatMatrix linspace (const FloatColumnVector& x1,
                                         const FloatColumnVector& x2,
diff --git a/liboctave/array/fNDArray.cc b/liboctave/array/fNDArray.cc
--- a/liboctave/array/fNDArray.cc
+++ b/liboctave/array/fNDArray.cc
@@ -635,51 +635,59 @@ FloatNDArray::sumsq (int dim) const
 }
 
 FloatNDArray
-FloatNDArray::max (int dim) const
+FloatNDArray::max (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<float> (*this, dim, mx_inline_max);
+  return do_mx_minmax_op<float> (*this, dim, mx_inline_max, includenan);
 }
 
 FloatNDArray
-FloatNDArray::max (Array<octave_idx_type>& idx_arg, int dim) const
+FloatNDArray::max (Array<octave_idx_type>& idx_arg,
+                   int dim, bool includenan) const
 {
-  return do_mx_minmax_op<float> (*this, idx_arg, dim, mx_inline_max);
+  return do_mx_minmax_op<float> (*this, idx_arg, dim, mx_inline_max,
+                                 includenan);
 }
 
 FloatNDArray
-FloatNDArray::min (int dim) const
+FloatNDArray::min (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<float> (*this, dim, mx_inline_min);
+  return do_mx_minmax_op<float> (*this, dim, mx_inline_min, includenan);
 }
 
 FloatNDArray
-FloatNDArray::min (Array<octave_idx_type>& idx_arg, int dim) const
+FloatNDArray::min (Array<octave_idx_type>& idx_arg,
+                   int dim, bool includenan) const
 {
-  return do_mx_minmax_op<float> (*this, idx_arg, dim, mx_inline_min);
+  return do_mx_minmax_op<float> (*this, idx_arg, dim, mx_inline_min,
+                                 includenan);
 }
 
 FloatNDArray
-FloatNDArray::cummax (int dim) const
+FloatNDArray::cummax (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<float> (*this, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<float> (*this, dim, mx_inline_cummax, includenan);
 }
 
 FloatNDArray
-FloatNDArray::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+FloatNDArray::cummax (Array<octave_idx_type>& idx_arg,
+                      int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<float> (*this, idx_arg, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<float> (*this, idx_arg, dim, mx_inline_cummax,
+                                    includenan);
 }
 
 FloatNDArray
-FloatNDArray::cummin (int dim) const
+FloatNDArray::cummin (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<float> (*this, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<float> (*this, dim, mx_inline_cummin, includenan);
 }
 
 FloatNDArray
-FloatNDArray::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+FloatNDArray::cummin (Array<octave_idx_type>& idx_arg,
+                      int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<float> (*this, idx_arg, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<float> (*this, idx_arg, dim, mx_inline_cummin,
+                                    includenan);
 }
 
 FloatNDArray
diff --git a/liboctave/array/fNDArray.h b/liboctave/array/fNDArray.h
--- a/liboctave/array/fNDArray.h
+++ b/liboctave/array/fNDArray.h
@@ -99,15 +99,19 @@ public:
   charNDArray concat (const charNDArray& rb,
                       const Array<octave_idx_type>& ra_idx);
 
-  FloatNDArray max (int dim = -1) const;
-  FloatNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
-  FloatNDArray min (int dim = -1) const;
-  FloatNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatNDArray max (int dim = -1, bool includenan = false) const;
+  FloatNDArray max (Array<octave_idx_type>& index,
+                    int dim = -1, bool includenan = false) const;
+  FloatNDArray min (int dim = -1, bool includenan = false) const;
+  FloatNDArray min (Array<octave_idx_type>& index,
+                    int dim = -1, bool includenan = false) const;
 
-  FloatNDArray cummax (int dim = -1) const;
-  FloatNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
-  FloatNDArray cummin (int dim = -1) const;
-  FloatNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
+  FloatNDArray cummax (int dim = -1, bool includenan = false) const;
+  FloatNDArray cummax (Array<octave_idx_type>& index,
+                       int dim = -1, bool includenan = false) const;
+  FloatNDArray cummin (int dim = -1, bool includenan = false) const;
+  FloatNDArray cummin (Array<octave_idx_type>& index,
+                       int dim = -1, bool includenan = false) const;
 
   FloatNDArray diff (octave_idx_type order = 1, int dim = -1) const;
 
diff --git a/liboctave/array/intNDArray.cc b/liboctave/array/intNDArray.cc
--- a/liboctave/array/intNDArray.cc
+++ b/liboctave/array/intNDArray.cc
@@ -239,58 +239,64 @@ intNDArray<T>::cumsum (int dim) const
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::max (int dim) const
+intNDArray<T>::max (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<T> (*this, dim, mx_inline_max);
+  return do_mx_minmax_op<T> (*this, dim, mx_inline_max, includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::max (Array<octave_idx_type>& idx_arg, int dim) const
+intNDArray<T>::max (Array<octave_idx_type>& idx_arg,
+                    int dim, bool includenan) const
 {
-  return do_mx_minmax_op<T> (*this, idx_arg, dim, mx_inline_max);
+  return do_mx_minmax_op<T> (*this, idx_arg, dim, mx_inline_max, includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::min (int dim) const
+intNDArray<T>::min (int dim, bool includenan) const
 {
-  return do_mx_minmax_op<T> (*this, dim, mx_inline_min);
+  return do_mx_minmax_op<T> (*this, dim, mx_inline_min, includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::min (Array<octave_idx_type>& idx_arg, int dim) const
+intNDArray<T>::min (Array<octave_idx_type>& idx_arg,
+                    int dim, bool includenan) const
 {
-  return do_mx_minmax_op<T> (*this, idx_arg, dim, mx_inline_min);
+  return do_mx_minmax_op<T> (*this, idx_arg, dim, mx_inline_min, includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::cummax (int dim) const
+intNDArray<T>::cummax (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<T> (*this, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<T> (*this, dim, mx_inline_cummax, includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::cummax (Array<octave_idx_type>& idx_arg, int dim) const
+intNDArray<T>::cummax (Array<octave_idx_type>& idx_arg,
+                       int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<T> (*this, idx_arg, dim, mx_inline_cummax);
+  return do_mx_cumminmax_op<T> (*this, idx_arg, dim, mx_inline_cummax,
+                                includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::cummin (int dim) const
+intNDArray<T>::cummin (int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<T> (*this, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<T> (*this, dim, mx_inline_cummin, includenan);
 }
 
 template <typename T>
 intNDArray<T>
-intNDArray<T>::cummin (Array<octave_idx_type>& idx_arg, int dim) const
+intNDArray<T>::cummin (Array<octave_idx_type>& idx_arg,
+                       int dim, bool includenan) const
 {
-  return do_mx_cumminmax_op<T> (*this, idx_arg, dim, mx_inline_cummin);
+  return do_mx_cumminmax_op<T> (*this, idx_arg, dim, mx_inline_cummin,
+                                includenan);
 }
 
 template <typename T>
diff --git a/liboctave/array/intNDArray.h b/liboctave/array/intNDArray.h
--- a/liboctave/array/intNDArray.h
+++ b/liboctave/array/intNDArray.h
@@ -81,15 +81,19 @@ public:
   boolNDArray all (int dim = -1) const;
   boolNDArray any (int dim = -1) const;
 
-  intNDArray max (int dim = -1) const;
-  intNDArray max (Array<octave_idx_type>& index, int dim = -1) const;
-  intNDArray min (int dim = -1) const;
-  intNDArray min (Array<octave_idx_type>& index, int dim = -1) const;
+  intNDArray max (int dim = -1, bool includenan = false) const;
+  intNDArray max (Array<octave_idx_type>& index,
+                  int dim = -1, bool includenan = false) const;
+  intNDArray min (int dim = -1, bool includenan = false) const;
+  intNDArray min (Array<octave_idx_type>& index,
+                  int dim = -1, bool includenan = false) const;
 
-  intNDArray cummax (int dim = -1) const;
-  intNDArray cummax (Array<octave_idx_type>& index, int dim = -1) const;
-  intNDArray cummin (int dim = -1) const;
-  intNDArray cummin (Array<octave_idx_type>& index, int dim = -1) const;
+  intNDArray cummax (int dim = -1, bool includenan = false) const;
+  intNDArray cummax (Array<octave_idx_type>& index,
+                     int dim = -1, bool includenan = false) const;
+  intNDArray cummin (int dim = -1, bool includenan = false) const;
+  intNDArray cummin (Array<octave_idx_type>& index,
+                     int dim = -1, bool includenan = false) const;
 
   intNDArray prod (int dim) const;
   intNDArray sum (int dim) const;
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, <=)
@@ -408,6 +399,60 @@ DEFMINMAXSPEC (double, mx_inline_xmax, >
 DEFMINMAXSPEC (float, mx_inline_xmin, <=)
 DEFMINMAXSPEC (float, mx_inline_xmax, >=)
 
+template <typename T>
+inline void
+mx_inline_xmin_inan (size_t n, T *r, const T *x, const T *y) throw ()
+{
+  for (size_t i = 0; i < n; i++)
+    r[i] = x[i] <= y[i] ? x[i] : (octave::math::isnan (x[i]) ? x[i] : y[i]);
+}
+
+template <typename T>
+inline void
+mx_inline_xmin_inan (size_t n, T *r, const T *x, T y) throw ()
+{
+  for (size_t i = 0; i < n; i++)
+    r[i] = x[i] <= y ? x[i] : (octave::math::isnan (x[i]) ? x[i] : y);
+}
+
+template <typename T>
+inline void
+mx_inline_xmin_inan (size_t n, T *r, T x, const T *y) throw ()
+{
+  if (octave::math::isnan (x))
+    std::fill_n (r, n, x);
+  else
+    for (size_t i = 0; i < n; i++)
+      r[i] = x <= y[i] ? x : y[i];
+}
+
+template <typename T>
+inline void
+mx_inline_xmax_inan (size_t n, T *r, const T *x, const T *y) throw ()
+{
+  for (size_t i = 0; i < n; i++)
+    r[i] = x[i] >= y[i] ? x[i] : (octave::math::isnan (x[i]) ? x[i] : y[i]);
+}
+
+template <typename T>
+inline void
+mx_inline_xmax_inan (size_t n, T *r, const T *x, T y) throw ()
+{
+  for (size_t i = 0; i < n; i++)
+    r[i] = x[i] >= y ? x[i] : (octave::math::isnan (x[i]) ? x[i] : y);
+}
+
+template <typename T>
+inline void
+mx_inline_xmax_inan (size_t n, T *r, T x, const T *y) throw ()
+{
+  if (octave::math::isnan (x))
+    std::fill_n (r, n, x);
+  else
+    for (size_t i = 0; i < n; i++)
+      r[i] = x >= y[i] ? x : y[i];
+}
+
 // Let the compiler decide which pow to use, whichever best matches the
 // arguments provided.
 
@@ -934,46 +979,89 @@ OP_CUM_FCNN (mx_inline_cumcount, bool, T
 
 #define OP_MINMAX_FCN(F, OP)                                            \
   template <typename T>                                                 \
-  void F (const T *v, T *r, octave_idx_type n)                          \
+  void F (const T *v, T *r, octave_idx_type n, bool includenan)         \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
     T tmp = v[0];                                                       \
     octave_idx_type i = 1;                                              \
-    if (octave::math::isnan (tmp))                                      \
+    if (includenan)                                                     \
       {                                                                 \
-        for (; i < n && octave::math::isnan (v[i]); i++) ;              \
-        if (i < n)                                                      \
-          tmp = v[i];                                                   \
+        if (! octave::math::isnan (tmp))                                \
+          {                                                             \
+            for (; i < n; i++)                                          \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    tmp = v[i];                                         \
+                    break;                                              \
+                  }                                                     \
+                if (v[i] OP tmp)                                        \
+                  tmp = v[i];                                           \
+              }                                                         \
+          }                                                             \
       }                                                                 \
-    for (; i < n; i++)                                                  \
-      if (v[i] OP tmp)                                                  \
-        tmp = v[i];                                                     \
+    else                                                                \
+      {                                                                 \
+        if (octave::math::isnan (tmp))                                  \
+          {                                                             \
+            for (; i < n && octave::math::isnan (v[i]); i++) ;          \
+            if (i < n)                                                  \
+              tmp = v[i];                                               \
+          }                                                             \
+        for (; i < n; i++)                                              \
+          if (v[i] OP tmp)                                              \
+            tmp = v[i];                                                 \
+      }                                                                 \
     *r = tmp;                                                           \
   }                                                                     \
   template <typename T>                                                 \
-  void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n)     \
+  void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n,     \
+          bool includenan)                                              \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
     T tmp = v[0];                                                       \
     octave_idx_type tmpi = 0;                                           \
     octave_idx_type i = 1;                                              \
-    if (octave::math::isnan (tmp))                                      \
+    if (includenan)                                                     \
       {                                                                 \
-        for (; i < n && octave::math::isnan (v[i]); i++) ;              \
-        if (i < n)                                                      \
+        if (! octave::math::isnan (tmp))                                \
           {                                                             \
-            tmp = v[i];                                                 \
-            tmpi = i;                                                   \
+            for (; i < n; i++)                                          \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    tmp = v[i];                                         \
+                    tmpi = i;                                           \
+                    break;                                              \
+                  }                                                     \
+                if (v[i] OP tmp)                                        \
+                  {                                                     \
+                    tmp = v[i];                                         \
+                    tmpi = i;                                           \
+                  }                                                     \
+              }                                                         \
           }                                                             \
       }                                                                 \
-    for (; i < n; i++)                                                  \
-      if (v[i] OP tmp)                                                  \
-        {                                                               \
-          tmp = v[i];                                                   \
-          tmpi = i;                                                     \
-        }                                                               \
+    else                                                                \
+      {                                                                 \
+        if (octave::math::isnan (tmp))                                  \
+          {                                                             \
+            for (; i < n && octave::math::isnan (v[i]); i++) ;          \
+            if (i < n)                                                  \
+              {                                                         \
+                tmp = v[i];                                             \
+                tmpi = i;                                               \
+              }                                                         \
+          }                                                             \
+        for (; i < n; i++)                                              \
+          if (v[i] OP tmp)                                              \
+            {                                                           \
+              tmp = v[i];                                               \
+              tmpi = i;                                                 \
+            }                                                           \
+      }                                                                 \
     *r = tmp;                                                           \
     *ri = tmpi;                                                         \
   }
@@ -988,86 +1076,170 @@ OP_MINMAX_FCN (mx_inline_max, >)
 #define OP_MINMAX_FCN2(F, OP)                                           \
   template <typename T>                                                 \
   inline void                                                           \
-  F (const T *v, T *r, octave_idx_type m, octave_idx_type n)            \
+  F (const T *v, T *r, octave_idx_type m, octave_idx_type n,            \
+     bool includenan)                                                   \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
-    bool nan = false;                                                   \
-    octave_idx_type j = 0;                                              \
-    for (octave_idx_type i = 0; i < m; i++)                             \
+    if (includenan)                                                     \
       {                                                                 \
-        r[i] = v[i];                                                    \
-        if (octave::math::isnan (v[i]))                                 \
-          nan = true;                                                   \
-      }                                                                 \
-    j++;                                                                \
-    v += m;                                                             \
-    while (nan && j < n)                                                \
-      {                                                                 \
-        nan = false;                                                    \
+        octave_idx_type prev = 0;                                       \
+        OCTAVE_LOCAL_BUFFER (octave_idx_type, next, m + 1);             \
+        next[0] = 0;                                                    \
         for (octave_idx_type i = 0; i < m; i++)                         \
           {                                                             \
+            r[i] = v[i];                                                \
+            if (! octave::math::isnan (v[i]))                           \
+              prev = i + 1;                                             \
+            next[prev] = i + 1;                                         \
+          }                                                             \
+        v += m;                                                         \
+        for (octave_idx_type j = 1; j < n; j++)                         \
+          {                                                             \
+            prev = 0;                                                   \
+            for (octave_idx_type i = next[0]; i < m; i = next[i + 1])   \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    r[i] = v[i];                                        \
+                    next[prev] = next[i + 1];                           \
+                  }                                                     \
+                else                                                    \
+                  {                                                     \
+                    if (v[i] OP r[i])                                   \
+                      r[i] = v[i];                                      \
+                    prev = i + 1;                                       \
+                  }                                                     \
+              }                                                         \
+            v += m;                                                     \
+          }                                                             \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        bool nan = false;                                               \
+        octave_idx_type j = 0;                                          \
+        for (octave_idx_type i = 0; i < m; i++)                         \
+          {                                                             \
+            r[i] = v[i];                                                \
             if (octave::math::isnan (v[i]))                             \
               nan = true;                                               \
-            else if (octave::math::isnan (r[i]) || v[i] OP r[i])        \
-              r[i] = v[i];                                              \
           }                                                             \
         j++;                                                            \
         v += m;                                                         \
-      }                                                                 \
-    while (j < n)                                                       \
-      {                                                                 \
-        for (octave_idx_type i = 0; i < m; i++)                         \
-          if (v[i] OP r[i])                                             \
-            r[i] = v[i];                                                \
-        j++;                                                            \
-        v += m;                                                         \
+        while (nan && j < n)                                            \
+          {                                                             \
+            nan = false;                                                \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    if (octave::math::isnan (r[i]))                     \
+                      nan = true;                                       \
+                  }                                                     \
+                else if (octave::math::isnan (r[i]) || v[i] OP r[i])    \
+                  r[i] = v[i];                                          \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+          }                                                             \
+        while (j < n)                                                   \
+          {                                                             \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              if (v[i] OP r[i])                                         \
+                r[i] = v[i];                                            \
+            j++;                                                        \
+            v += m;                                                     \
+          }                                                             \
       }                                                                 \
   }                                                                     \
   template <typename T>                                                 \
   inline void                                                           \
   F (const T *v, T *r, octave_idx_type *ri,                             \
-     octave_idx_type m, octave_idx_type n)                              \
+     octave_idx_type m, octave_idx_type n, bool includenan)             \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
-    bool nan = false;                                                   \
-    octave_idx_type j = 0;                                              \
-    for (octave_idx_type i = 0; i < m; i++)                             \
+    if (includenan)                                                     \
       {                                                                 \
-        r[i] = v[i];                                                    \
-        ri[i] = j;                                                      \
-        if (octave::math::isnan (v[i]))                                 \
-          nan = true;                                                   \
-      }                                                                 \
-    j++;                                                                \
-    v += m;                                                             \
-    while (nan && j < n)                                                \
-      {                                                                 \
-        nan = false;                                                    \
+        octave_idx_type prev = 0;                                       \
+        OCTAVE_LOCAL_BUFFER (octave_idx_type, next, m + 1);             \
+        next[0] = 0;                                                    \
         for (octave_idx_type i = 0; i < m; i++)                         \
           {                                                             \
+            r[i] = v[i];                                                \
+            ri[i] = 0;                                                  \
+            if (! octave::math::isnan (v[i]))                           \
+              prev = i + 1;                                             \
+            next[prev] = i + 1;                                         \
+          }                                                             \
+        v += m;                                                         \
+        for (octave_idx_type j = 1; j < n; j++)                         \
+          {                                                             \
+            prev = 0;                                                   \
+            for (octave_idx_type i = next[0]; i < m; i = next[i + 1])   \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    r[i] = v[i];                                        \
+                    ri[i] = j;                                          \
+                    next[prev] = next[i + 1];                           \
+                  }                                                     \
+                else                                                    \
+                  {                                                     \
+                    if (v[i] OP r[i])                                   \
+                      {                                                 \
+                        r[i] = v[i];                                    \
+                        ri[i] = j;                                      \
+                      }                                                 \
+                    prev = i + 1;                                       \
+                  }                                                     \
+              }                                                         \
+            v += m;                                                     \
+          }                                                             \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        bool nan = false;                                               \
+        octave_idx_type j = 0;                                          \
+        for (octave_idx_type i = 0; i < m; i++)                         \
+          {                                                             \
+            r[i] = v[i];                                                \
+            ri[i] = j;                                                  \
             if (octave::math::isnan (v[i]))                             \
               nan = true;                                               \
-            else if (octave::math::isnan (r[i]) || v[i] OP r[i])        \
-              {                                                         \
-                r[i] = v[i];                                            \
-                ri[i] = j;                                              \
-              }                                                         \
           }                                                             \
         j++;                                                            \
         v += m;                                                         \
-      }                                                                 \
-    while (j < n)                                                       \
-      {                                                                 \
-        for (octave_idx_type i = 0; i < m; i++)                         \
-          if (v[i] OP r[i])                                             \
-            {                                                           \
-              r[i] = v[i];                                              \
-              ri[i] = j;                                                \
-            }                                                           \
-        j++;                                                            \
-        v += m;                                                         \
+        while (nan && j < n)                                            \
+          {                                                             \
+            nan = false;                                                \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    if (octave::math::isnan (r[i]))                     \
+                      nan = true;                                       \
+                  }                                                     \
+                else if (octave::math::isnan (r[i]) || v[i] OP r[i])    \
+                  {                                                     \
+                    r[i] = v[i];                                        \
+                    ri[i] = j;                                          \
+                  }                                                     \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+          }                                                             \
+        while (j < n)                                                   \
+          {                                                             \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              if (v[i] OP r[i])                                         \
+                {                                                       \
+                  r[i] = v[i];                                          \
+                  ri[i] = j;                                            \
+                }                                                       \
+            j++;                                                        \
+            v += m;                                                     \
+          }                                                             \
       }                                                                 \
   }
 
@@ -1078,7 +1250,8 @@ OP_MINMAX_FCN2 (mx_inline_max, >)
   template <typename T>                                         \
   inline void                                                   \
   F (const T *v, T *r, octave_idx_type l,                       \
-     octave_idx_type n, octave_idx_type u)                      \
+     octave_idx_type n, octave_idx_type u,                      \
+     bool includenan)                                           \
   {                                                             \
     if (! n)                                                    \
       return;                                                   \
@@ -1086,7 +1259,7 @@ OP_MINMAX_FCN2 (mx_inline_max, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, n);                                        \
+            F (v, r, n, includenan);                            \
             v += n;                                             \
             r++;                                                \
           }                                                     \
@@ -1095,7 +1268,7 @@ OP_MINMAX_FCN2 (mx_inline_max, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, l, n);                                     \
+            F (v, r, l, n, includenan);                         \
             v += l*n;                                           \
             r += l;                                             \
           }                                                     \
@@ -1104,14 +1277,15 @@ OP_MINMAX_FCN2 (mx_inline_max, >)
   template <typename T>                                         \
   inline void                                                   \
   F (const T *v, T *r, octave_idx_type *ri,                     \
-     octave_idx_type l, octave_idx_type n, octave_idx_type u)   \
+     octave_idx_type l, octave_idx_type n, octave_idx_type u,   \
+     bool includenan)                                           \
   {                                                             \
     if (! n) return;                                            \
     if (l == 1)                                                 \
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, ri, n);                                    \
+            F (v, r, ri, n, includenan);                        \
             v += n;                                             \
             r++;                                                \
             ri++;                                               \
@@ -1121,7 +1295,7 @@ OP_MINMAX_FCN2 (mx_inline_max, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, ri, l, n);                                 \
+            F (v, r, ri, l, n, includenan);                     \
             v += l*n;                                           \
             r += l;                                             \
             ri += l;                                            \
@@ -1134,7 +1308,7 @@ OP_MINMAX_FCNN (mx_inline_max)
 
 #define OP_CUMMINMAX_FCN(F, OP)                                         \
   template <typename T>                                                 \
-  void F (const T *v, T *r, octave_idx_type n)                          \
+  void F (const T *v, T *r, octave_idx_type n, bool includenan)         \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
@@ -1143,24 +1317,54 @@ OP_MINMAX_FCNN (mx_inline_max)
     octave_idx_type j = 0;                                              \
     if (octave::math::isnan (tmp))                                      \
       {                                                                 \
+        if (includenan)                                                 \
+          {                                                             \
+            for (; j < n; j++)                                          \
+              r[j] = tmp;                                               \
+            return;                                                     \
+          }                                                             \
         for (; i < n && octave::math::isnan (v[i]); i++) ;              \
         for (; j < i; j++)                                              \
           r[j] = tmp;                                                   \
         if (i < n)                                                      \
           tmp = v[i];                                                   \
       }                                                                 \
-    for (; i < n; i++)                                                  \
-      if (v[i] OP tmp)                                                  \
-        {                                                               \
-          for (; j < i; j++)                                            \
-            r[j] = tmp;                                                 \
-          tmp = v[i];                                                   \
-        }                                                               \
+    if (includenan)                                                     \
+      {                                                                 \
+        for (; i < n; i++)                                              \
+          {                                                             \
+            if (v[i] OP tmp)                                            \
+              {                                                         \
+                for (; j < i; j++)                                      \
+                  r[j] = tmp;                                           \
+                tmp = v[i];                                             \
+              }                                                         \
+            else if (octave::math::isnan (v[i]))                        \
+              {                                                         \
+                for (; j < i; j++)                                      \
+                  r[j] = tmp;                                           \
+                for (; j < n; j++)                                      \
+                  r[j] = v[i];                                          \
+                return;                                                 \
+              }                                                         \
+          }                                                             \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        for (; i < n; i++)                                              \
+          if (v[i] OP tmp)                                              \
+            {                                                           \
+              for (; j < i; j++)                                        \
+                r[j] = tmp;                                             \
+              tmp = v[i];                                               \
+            }                                                           \
+      }                                                                 \
     for (; j < i; j++)                                                  \
       r[j] = tmp;                                                       \
   }                                                                     \
   template <typename T>                                                 \
-  void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n)     \
+  void F (const T *v, T *r, octave_idx_type *ri, octave_idx_type n,     \
+          bool includenan)                                              \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
@@ -1170,6 +1374,15 @@ OP_MINMAX_FCNN (mx_inline_max)
     octave_idx_type j = 0;                                              \
     if (octave::math::isnan (tmp))                                      \
       {                                                                 \
+        if (includenan)                                                 \
+          {                                                             \
+            for (; j < n; j++)                                          \
+              {                                                         \
+                r[j] = tmp;                                             \
+                ri[j] = 0;                                              \
+              }                                                         \
+            return;                                                     \
+          }                                                             \
         for (; i < n && octave::math::isnan (v[i]); i++) ;              \
         for (; j < i; j++)                                              \
           {                                                             \
@@ -1182,17 +1395,50 @@ OP_MINMAX_FCNN (mx_inline_max)
             tmpi = i;                                                   \
           }                                                             \
       }                                                                 \
-    for (; i < n; i++)                                                  \
-      if (v[i] OP tmp)                                                  \
-        {                                                               \
-          for (; j < i; j++)                                            \
+    if (includenan)                                                     \
+      {                                                                 \
+        for (; i < n; i++)                                              \
+          {                                                             \
+            if (v[i] OP tmp)                                            \
+              {                                                         \
+                for (; j < i; j++)                                      \
+                  {                                                     \
+                    r[j] = tmp;                                         \
+                    ri[j] = tmpi;                                       \
+                  }                                                     \
+                tmp = v[i];                                             \
+                tmpi = i;                                               \
+              }                                                         \
+            else if (octave::math::isnan (v[i]))                        \
+              {                                                         \
+                for (; j < i; j++)                                      \
+                  {                                                     \
+                    r[j] = tmp;                                         \
+                    ri[j] = tmpi;                                       \
+                  }                                                     \
+                for (; j < n; j++)                                      \
+                  {                                                     \
+                    r[j] = v[i];                                        \
+                    ri[j] = i;                                          \
+                  }                                                     \
+                return;                                                 \
+              }                                                         \
+          }                                                             \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        for (; i < n; i++)                                              \
+          if (v[i] OP tmp)                                              \
             {                                                           \
-              r[j] = tmp;                                               \
-              ri[j] = tmpi;                                             \
+              for (; j < i; j++)                                        \
+                {                                                       \
+                  r[j] = tmp;                                           \
+                  ri[j] = tmpi;                                         \
+                }                                                       \
+              tmp = v[i];                                               \
+              tmpi = i;                                                 \
             }                                                           \
-          tmp = v[i];                                                   \
-          tmpi = i;                                                     \
-        }                                                               \
+      }                                                                 \
     for (; j < i; j++)                                                  \
       {                                                                 \
         r[j] = tmp;                                                     \
@@ -1210,100 +1456,103 @@ OP_CUMMINMAX_FCN (mx_inline_cummax, >)
 #define OP_CUMMINMAX_FCN2(F, OP)                                        \
   template <typename T>                                                 \
   inline void                                                           \
-  F (const T *v, T *r, octave_idx_type m, octave_idx_type n)            \
+  F (const T *v, T *r, octave_idx_type m, octave_idx_type n,            \
+     bool includenan)                                                   \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
-    bool nan = false;                                                   \
     const T *r0;                                                        \
     octave_idx_type j = 0;                                              \
-    for (octave_idx_type i = 0; i < m; i++)                             \
+    if (includenan)                                                     \
       {                                                                 \
-        r[i] = v[i];                                                    \
-        if (octave::math::isnan (v[i]))                                 \
-          nan = true;                                                   \
+        for (octave_idx_type i = 0; i < m; i++)                         \
+          r[i] = v[i];                                                  \
+        j++;                                                            \
+        v += m;                                                         \
+        r0 = r;                                                         \
+        r += m;                                                         \
+        while (j < n)                                                   \
+          {                                                             \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (v[i] OP r0[i] || (octave::math::isnan (v[i])        \
+                                      && ! octave::math::isnan (r0[i])))\
+                  r[i] = v[i];                                          \
+                else                                                    \
+                  r[i] = r0[i];                                         \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+            r0 = r;                                                     \
+            r += m;                                                     \
+          }                                                             \
       }                                                                 \
-    j++;                                                                \
-    v += m;                                                             \
-    r0 = r;                                                             \
-    r += m;                                                             \
-    while (nan && j < n)                                                \
+    else                                                                \
       {                                                                 \
-        nan = false;                                                    \
+        bool nan = false;                                               \
         for (octave_idx_type i = 0; i < m; i++)                         \
           {                                                             \
+            r[i] = v[i];                                                \
             if (octave::math::isnan (v[i]))                             \
-              {                                                         \
-                r[i] = r0[i];                                           \
-                nan = true;                                             \
-              }                                                         \
-            else if (octave::math::isnan (r0[i]) || v[i] OP r0[i])      \
-              r[i] = v[i];                                              \
-            else                                                        \
-              r[i] = r0[i];                                             \
+              nan = true;                                               \
           }                                                             \
         j++;                                                            \
         v += m;                                                         \
         r0 = r;                                                         \
         r += m;                                                         \
-      }                                                                 \
-    while (j < n)                                                       \
-      {                                                                 \
-        for (octave_idx_type i = 0; i < m; i++)                         \
-          if (v[i] OP r0[i])                                            \
-            r[i] = v[i];                                                \
-          else                                                          \
-            r[i] = r0[i];                                               \
-        j++;                                                            \
-        v += m;                                                         \
-        r0 = r;                                                         \
-        r += m;                                                         \
+        while (nan && j < n)                                            \
+          {                                                             \
+            nan = false;                                                \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    r[i] = r0[i];                                       \
+                    if (octave::math::isnan (r[i]))                     \
+                      nan = true;                                       \
+                  }                                                     \
+                else if (octave::math::isnan (r0[i]) || v[i] OP r0[i])  \
+                  r[i] = v[i];                                          \
+                else                                                    \
+                  r[i] = r0[i];                                         \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+            r0 = r;                                                     \
+            r += m;                                                     \
+          }                                                             \
+        while (j < n)                                                   \
+          {                                                             \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (v[i] OP r0[i])                                      \
+                  r[i] = v[i];                                          \
+                else                                                    \
+                  r[i] = r0[i];                                         \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+            r0 = r;                                                     \
+            r += m;                                                     \
+          }                                                             \
       }                                                                 \
   }                                                                     \
   template <typename T>                                                 \
   inline void                                                           \
   F (const T *v, T *r, octave_idx_type *ri,                             \
-     octave_idx_type m, octave_idx_type n)                              \
+     octave_idx_type m, octave_idx_type n, bool includenan)             \
   {                                                                     \
     if (! n)                                                            \
       return;                                                           \
-    bool nan = false;                                                   \
     const T *r0;                                                        \
     const octave_idx_type *r0i;                                         \
     octave_idx_type j = 0;                                              \
-    for (octave_idx_type i = 0; i < m; i++)                             \
+    if (includenan)                                                     \
       {                                                                 \
-        r[i] = v[i]; ri[i] = 0;                                         \
-        if (octave::math::isnan (v[i]))                                 \
-          nan = true;                                                   \
-      }                                                                 \
-    j++;                                                                \
-    v += m;                                                             \
-    r0 = r;                                                             \
-    r += m;                                                             \
-    r0i = ri;                                                           \
-    ri += m;                                                            \
-    while (nan && j < n)                                                \
-      {                                                                 \
-        nan = false;                                                    \
         for (octave_idx_type i = 0; i < m; i++)                         \
           {                                                             \
-            if (octave::math::isnan (v[i]))                             \
-              {                                                         \
-                r[i] = r0[i];                                           \
-                ri[i] = r0i[i];                                         \
-                nan = true;                                             \
-              }                                                         \
-            else if (octave::math::isnan (r0[i]) || v[i] OP r0[i])      \
-              {                                                         \
-                r[i] = v[i];                                            \
-                ri[i] = j;                                              \
-              }                                                         \
-            else                                                        \
-              {                                                         \
-                r[i] = r0[i];                                           \
-                ri[i] = r0i[i];                                         \
-              }                                                         \
+            r[i] = v[i];                                                \
+            ri[i] = 0;                                                  \
           }                                                             \
         j++;                                                            \
         v += m;                                                         \
@@ -1311,26 +1560,97 @@ OP_CUMMINMAX_FCN (mx_inline_cummax, >)
         r += m;                                                         \
         r0i = ri;                                                       \
         ri += m;                                                        \
+        while (j < n)                                                   \
+          {                                                             \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (v[i] OP r0[i] || (octave::math::isnan (v[i])        \
+                                      && ! octave::math::isnan (r0[i])))\
+                  {                                                     \
+                    r[i] = v[i];                                        \
+                    ri[i] = j;                                          \
+                  }                                                     \
+                else                                                    \
+                  {                                                     \
+                    r[i] = r0[i];                                       \
+                    ri[i] = r0i[i];                                     \
+                  }                                                     \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+            r0 = r;                                                     \
+            r += m;                                                     \
+            r0i = ri;                                                   \
+            ri += m;                                                    \
+          }                                                             \
       }                                                                 \
-    while (j < n)                                                       \
+    else                                                                \
       {                                                                 \
+        bool nan = false;                                               \
         for (octave_idx_type i = 0; i < m; i++)                         \
-          if (v[i] OP r0[i])                                            \
-            {                                                           \
-              r[i] = v[i];                                              \
-              ri[i] = j;                                                \
-            }                                                           \
-          else                                                          \
-            {                                                           \
-              r[i] = r0[i];                                             \
-              ri[i] = r0i[i];                                           \
-            }                                                           \
+          {                                                             \
+            r[i] = v[i]; ri[i] = 0;                                     \
+            if (octave::math::isnan (v[i]))                             \
+              nan = true;                                               \
+          }                                                             \
         j++;                                                            \
         v += m;                                                         \
         r0 = r;                                                         \
         r += m;                                                         \
         r0i = ri;                                                       \
         ri += m;                                                        \
+        while (nan && j < n)                                            \
+          {                                                             \
+            nan = false;                                                \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (octave::math::isnan (v[i]))                         \
+                  {                                                     \
+                    r[i] = r0[i];                                       \
+                    ri[i] = r0i[i];                                     \
+                    if (octave::math::isnan (r[i]))                     \
+                      nan = true;                                       \
+                  }                                                     \
+                else if (octave::math::isnan (r0[i]) || v[i] OP r0[i])  \
+                  {                                                     \
+                    r[i] = v[i];                                        \
+                    ri[i] = j;                                          \
+                  }                                                     \
+                else                                                    \
+                  {                                                     \
+                    r[i] = r0[i];                                       \
+                    ri[i] = r0i[i];                                     \
+                  }                                                     \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+            r0 = r;                                                     \
+            r += m;                                                     \
+            r0i = ri;                                                   \
+            ri += m;                                                    \
+          }                                                             \
+        while (j < n)                                                   \
+          {                                                             \
+            for (octave_idx_type i = 0; i < m; i++)                     \
+              {                                                         \
+                if (v[i] OP r0[i])                                      \
+                  {                                                     \
+                    r[i] = v[i];                                        \
+                    ri[i] = j;                                          \
+                  }                                                     \
+                else                                                    \
+                  {                                                     \
+                    r[i] = r0[i];                                       \
+                    ri[i] = r0i[i];                                     \
+                  }                                                     \
+              }                                                         \
+            j++;                                                        \
+            v += m;                                                     \
+            r0 = r;                                                     \
+            r += m;                                                     \
+            r0i = ri;                                                   \
+            ri += m;                                                    \
+          }                                                             \
       }                                                                 \
   }
 
@@ -1341,7 +1661,7 @@ OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
   template <typename T>                                         \
   inline void                                                   \
   F (const T *v, T *r, octave_idx_type l,                       \
-     octave_idx_type n, octave_idx_type u)                      \
+     octave_idx_type n, octave_idx_type u, bool includenan)     \
   {                                                             \
     if (! n)                                                    \
       return;                                                   \
@@ -1349,7 +1669,7 @@ OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, n);                                        \
+            F (v, r, n, includenan);                            \
             v += n;                                             \
             r += n;                                             \
           }                                                     \
@@ -1358,7 +1678,7 @@ OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, l, n);                                     \
+            F (v, r, l, n, includenan);                         \
             v += l*n;                                           \
             r += l*n;                                           \
           }                                                     \
@@ -1367,7 +1687,8 @@ OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
   template <typename T>                                         \
   inline void                                                   \
   F (const T *v, T *r, octave_idx_type *ri,                     \
-     octave_idx_type l, octave_idx_type n, octave_idx_type u)   \
+     octave_idx_type l, octave_idx_type n, octave_idx_type u,   \
+     bool includenan)                                           \
   {                                                             \
     if (! n)                                                    \
       return;                                                   \
@@ -1375,7 +1696,7 @@ OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, ri, n);                                    \
+            F (v, r, ri, n, includenan);                        \
             v += n;                                             \
             r += n;                                             \
             ri += n;                                            \
@@ -1385,7 +1706,7 @@ OP_CUMMINMAX_FCN2 (mx_inline_cummax, >)
       {                                                         \
         for (octave_idx_type i = 0; i < u; i++)                 \
           {                                                     \
-            F (v, r, ri, l, n);                                 \
+            F (v, r, ri, l, n, includenan);                     \
             v += l*n;                                           \
             r += l*n;                                           \
             ri += l*n;                                          \
@@ -1580,7 +1901,8 @@ template <typename R>
 inline Array<R>
 do_mx_minmax_op (const Array<R>& src, int dim,
                  void (*mx_minmax_op) (const R *, R *, octave_idx_type,
-                                       octave_idx_type, octave_idx_type))
+                                       octave_idx_type, octave_idx_type, bool),
+                 bool includenan = false)
 {
   octave_idx_type l, n, u;
   dim_vector dims = src.dims ();
@@ -1591,7 +1913,7 @@ do_mx_minmax_op (const Array<R>& src, in
   dims.chop_trailing_singletons ();
 
   Array<R> ret (dims);
-  mx_minmax_op (src.data (), ret.fortran_vec (), l, n, u);
+  mx_minmax_op (src.data (), ret.fortran_vec (), l, n, u, includenan);
 
   return ret;
 }
@@ -1600,7 +1922,9 @@ template <typename R>
 inline Array<R>
 do_mx_minmax_op (const Array<R>& src, Array<octave_idx_type>& idx, int dim,
                  void (*mx_minmax_op) (const R *, R *, octave_idx_type *,
-                                       octave_idx_type, octave_idx_type, octave_idx_type))
+                                       octave_idx_type, octave_idx_type,
+                                       octave_idx_type, bool),
+                 bool includenan = false)
 {
   octave_idx_type l, n, u;
   dim_vector dims = src.dims ();
@@ -1614,7 +1938,7 @@ do_mx_minmax_op (const Array<R>& src, Ar
   if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
 
   mx_minmax_op (src.data (), ret.fortran_vec (), idx.fortran_vec (),
-                l, n, u);
+                l, n, u, includenan);
 
   return ret;
 }
@@ -1623,14 +1947,16 @@ template <typename R>
 inline Array<R>
 do_mx_cumminmax_op (const Array<R>& src, int dim,
                     void (*mx_cumminmax_op) (const R *, R *, octave_idx_type,
-                                             octave_idx_type, octave_idx_type))
+                                             octave_idx_type, octave_idx_type,
+                                             bool),
+                    bool includenan = false)
 {
   octave_idx_type l, n, u;
   dim_vector dims = src.dims ();
   get_extent_triplet (dims, dim, l, n, u);
 
   Array<R> ret (dims);
-  mx_cumminmax_op (src.data (), ret.fortran_vec (), l, n, u);
+  mx_cumminmax_op (src.data (), ret.fortran_vec (), l, n, u, includenan);
 
   return ret;
 }
@@ -1639,7 +1965,9 @@ template <typename R>
 inline Array<R>
 do_mx_cumminmax_op (const Array<R>& src, Array<octave_idx_type>& idx, int dim,
                     void (*mx_cumminmax_op) (const R *, R *, octave_idx_type *,
-                                             octave_idx_type, octave_idx_type, octave_idx_type))
+                                             octave_idx_type, octave_idx_type,
+                                             octave_idx_type, bool),
+                    bool includenan = false)
 {
   octave_idx_type l, n, u;
   dim_vector dims = src.dims ();
@@ -1649,7 +1977,7 @@ do_mx_cumminmax_op (const Array<R>& src,
   if (idx.dims () != dims) idx = Array<octave_idx_type> (dims);
 
   mx_cumminmax_op (src.data (), ret.fortran_vec (), idx.fortran_vec (),
-                   l, n, u);
+                   l, n, u, includenan);
 
   return ret;
 }
diff --git a/liboctave/operators/mx-op-decl.h b/liboctave/operators/mx-op-decl.h
--- a/liboctave/operators/mx-op-decl.h
+++ b/liboctave/operators/mx-op-decl.h
@@ -286,12 +286,12 @@ class boolNDArray;
 // scalar by N-D array min/max ops
 
 #define MINMAX_DECLS(T, S, API)                 \
-  extern API T min (S d, const T& m);           \
-  extern API T min (const T& m, S d);           \
-  extern API T min (const T& a, const T& b);    \
-  extern API T max (S d, const T& m);           \
-  extern API T max (const T& m, S d);           \
-  extern API T max (const T& a, const T& b);
+  extern API T min (S d, const T& m, bool includenan = false);           \
+  extern API T min (const T& m, S d, bool includenan = false);           \
+  extern API T min (const T& a, const T& b, bool includenan = false);    \
+  extern API T max (S d, const T& m, bool includenan = false);           \
+  extern API T max (const T& m, S d, bool includenan = false);           \
+  extern API T max (const T& a, const T& b, bool includenan = false);
 
 // permutation matrix by matrix ops and vice versa
 
diff --git a/liboctave/operators/mx-op-defs.h b/liboctave/operators/mx-op-defs.h
--- a/liboctave/operators/mx-op-defs.h
+++ b/liboctave/operators/mx-op-defs.h
@@ -562,23 +562,48 @@ along with Octave; see the file COPYING.
 
 #define SND_MINMAX_FCN(FCN, OP, T, S)                                   \
   T                                                                     \
-  FCN (S d, const T& m)                                                 \
+  FCN (S d, const T& m, bool includenan)                                \
   {                                                                     \
-    return do_sm_binary_op<T::element_type, S, T::element_type> (d, m, mx_inline_x##FCN); \
+    void (*op) (size_t, T::element_type *, S,                           \
+                const T::element_type *) throw ();                      \
+                                                                        \
+    if (includenan)                                                     \
+      op = mx_inline_x##FCN##_inan;                                     \
+    else                                                                \
+      op = mx_inline_x##FCN;                                            \
+                                                                        \
+    return do_sm_binary_op<T::element_type, S, T::element_type> (d, m,  \
+                                                                 op);   \
   }
 
 #define NDS_MINMAX_FCN(FCN, OP, T, S)                                   \
   T                                                                     \
-  FCN (const T& m, S d)                                                 \
+  FCN (const T& m, S d, bool includenan)                                \
   {                                                                     \
-    return do_ms_binary_op<T::element_type, T::element_type, S> (m, d, mx_inline_x##FCN); \
+    void (*op) (size_t, T::element_type *,                              \
+                const T::element_type *, S) throw ();                   \
+                                                                        \
+    if (includenan)                                                     \
+      op = mx_inline_x##FCN##_inan;                                     \
+    else                                                                \
+      op = mx_inline_x##FCN;                                            \
+                                                                        \
+    return do_ms_binary_op<T::element_type, T::element_type, S> (m, d,  \
+                                                                 op);   \
   }
 
 #define NDND_MINMAX_FCN(FCN, OP, T, S)                                  \
   T                                                                     \
-  FCN (const T& a, const T& b)                                          \
+  FCN (const T& a, const T& b, bool includenan)                         \
   {                                                                     \
-    return do_mm_binary_op<T::element_type, T::element_type, T::element_type> (a, b, mx_inline_x##FCN, mx_inline_x##FCN, mx_inline_x##FCN, #FCN); \
+    if (includenan)                                                     \
+      return do_mm_binary_op<T::element_type, T::element_type, T::element_type> (a, b, mx_inline_x##FCN##_inan, \
+                                                                                 mx_inline_x##FCN##_inan, \
+                                                                                 mx_inline_x##FCN##_inan, #FCN); \
+                                                                        \
+    return do_mm_binary_op<T::element_type, T::element_type, T::element_type> (a, b, mx_inline_x##FCN, \
+                                                                               mx_inline_x##FCN, \
+                                                                               mx_inline_x##FCN, #FCN); \
   }
 
 #define MINMAX_FCNS(T, S)                       \
