# HG changeset patch
# User David Spies <dnspies@gmail.com>
# Date 1401003429 21600
#      Sun May 25 01:37:09 2014 -0600
# Node ID df0221588028f7c21875245a4483343c98ae1730
# Parent  330cb3efc6baaff4e3033fa976952636f644bb01
Implemented a generic version of the find method which in turn uses templated iterators to handle the differences between matrix types

diff -r 330cb3efc6ba -r df0221588028 libinterp/corefcn/find.cc
--- a/libinterp/corefcn/find.cc	Sat May 24 14:14:21 2014 -0600
+++ b/libinterp/corefcn/find.cc	Sun May 25 01:37:09 2014 -0600
@@ -20,300 +20,183 @@
 
 */
 
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include <cassert>
 
-#include "quit.h"
+#include "find.h"
 
 #include "defun.h"
 #include "error.h"
 #include "gripes.h"
 #include "oct-obj.h"
 
-// Find at most N_TO_FIND nonzero elements in NDA.  Search forward if
-// DIRECTION is 1, backward if it is -1.  NARGOUT is the number of
-// output arguments.  If N_TO_FIND is -1, find all nonzero elements.
+/*
+ * res_type_of is a generic type used for storing the result of
+ * a find operation.  The way in which this result is stored will
+ * vary based on whether the number of requested return values
+ * is 1, 2, or 3.
+ * Each instantiation of res_type_of must support a couple different
+ * operations.  It is constructed with a dim_vector which indicates
+ * the size of the return vectors (the size of all return values for
+ * find is the same).  It supports an add() operation which
+ * the generic "find" method will call to add an element to the return
+ * values.  "add" takes an index and an nz-iterator for the matrix type
+ * being searched (from nz-iterators.h)
+ * Finally it supports get_list which returns an octave_value_list of the
+ * return values (1,2, or 3 depending on the nargout template argument)
+ */
 
-template <typename T>
-octave_value_list
-find_nonzero_elem_idx (const Array<T>& nda, int nargout,
-                       octave_idx_type n_to_find, int direction)
-{
-  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
+template<int nargout, typename T>
+struct res_type_of;
 
-  Array<octave_idx_type> idx;
-  if (n_to_find >= 0)
-    idx = nda.find (n_to_find, direction == -1);
-  else
-    idx = nda.find ();
+template<typename T>
+struct res_type_of<1, T>
+  {
+  res_type_of()
+    {
+    }
+  res_type_of(const dim_vector& nnz) :
+      res (nnz)
+    {
+    }
+  Array<double> res;
+  template<typename iter_t>
+  void add(octave_idx_type place, iter_t iter)
+    {
+      res.xelem (place) = iter.interp_idx ();
+    }
+  octave_value_list get_list()
+    {
+      return octave_value_list (octave_value (res));
+    }
+  };
 
-  // The maximum element is always at the end.
-  octave_idx_type iext = idx.is_empty () ? 0 : idx.xelem (idx.numel () - 1) + 1;
+template<typename T>
+struct res_type_of<2, T>
+  {
+  res_type_of()
+    {
+    }
+  res_type_of(const dim_vector& nnz) :
+      rescol (nnz), resrow (nnz)
+    {
+    }
+  Array<double> rescol;
+  Array<double> resrow;
+  template<typename iter_t>
+  void add(octave_idx_type place, iter_t iter)
+    {
+      rescol.xelem (place) = to_interp_idx(iter.col ());
+      resrow.xelem (place) = to_interp_idx(iter.row ());
+    }
+  octave_value_list get_list()
+    {
+      octave_value_list res (2);
+      res.xelem (0) = resrow;
+      res.xelem (1) = rescol;
+      return res;
+    }
+  };
 
-  switch (nargout)
+template<typename T>
+struct res_type_of<3, T>
+  {
+  res_type_of()
     {
-    default:
-    case 3:
-      retval(2) = Array<T> (nda.index (idx_vector (idx)));
-      // Fall through!
+    }
+  res_type_of(const dim_vector& nnz) :
+      rescol (nnz), resrow (nnz), elems (nnz)
+    {
+    }
+  Array<double> rescol;
+  Array<double> resrow;
+  Array<T> elems;
+  template<typename iter_t>
+  void add(octave_idx_type place, iter_t iter)
+    {
+      rescol.xelem (place) = to_interp_idx(iter.col ());
+      resrow.xelem (place) = to_interp_idx(iter.row ());
+      elems.xelem (place) = iter.data ();
+    }
+  octave_value_list get_list()
+    {
+      octave_value_list res (3);
+      res.xelem (0) = resrow;
+      res.xelem (1) = rescol;
+      res.xelem (2) = elems;
+      return res;
+    }
+  };
 
-    case 2:
-      {
-        Array<octave_idx_type> jdx (idx.dims ());
-        octave_idx_type n = idx.length ();
-        octave_idx_type nr = nda.rows ();
-        for (octave_idx_type i = 0; i < n; i++)
+
+/*
+ * Calls the find_to_iter template method in "find.h" which
+ * in turn will fill in resvec of type res_type_of (see above).
+ * This is generic enough to work for any matrix type S
+ */
+template<int direction, int nargout, typename S>
+octave_value_list find_call(const S& v,
+    octave_idx_type n_to_find)
+  {
+    typedef typename S::element_type T;
+    typedef res_type_of<nargout, T> R;
+
+    R resvec;
+    dir_constant<direction> dirc;
+    find_to_iter (dirc, v, n_to_find, resvec);
+
+    return resvec.get_list ();
+  }
+
+/*
+ * A large annoying plethora of switch statements necessary to
+ * transform normal arguments into template arguments. These must
+ * be template arguments in order to construct the proper res_type_of
+ */
+template<typename M>
+octave_value_list find_templated(const M& v, int nargout,
+		octave_idx_type n_to_find, int direction)
+    {
+      switch (nargout)
+        {
+      case 1:
+        switch (direction)
           {
-            jdx.xelem (i) = idx.xelem (i) / nr;
-            idx.xelem (i) %= nr;
+        case -1:
+          return find_call<-1, 1> (v, n_to_find);
+        case 1:
+          return find_call<1, 1> (v, n_to_find);
+        default:
+          assert(false); //Direction can only be 1 or -1
           }
-        iext = -1;
-        retval(1) = idx_vector (jdx, -1);
-      }
-      // Fall through!
-
-    case 1:
-    case 0:
-      retval(0) = idx_vector (idx, iext);
-      break;
+        break;
+      case 2:
+        switch (direction)
+          {
+        case -1:
+          return find_call<-1, 2> (v, n_to_find);
+        case 1:
+          return find_call<1, 2> (v, n_to_find);
+        default:
+          assert(false); //Direction can only be 1 or -1
+          }
+        break;
+      case 3:
+        switch (direction)
+          {
+        case -1:
+          return find_call<-1, 3> (v, n_to_find);
+        case 1:
+          return find_call<1, 3> (v, n_to_find);
+        default:
+          assert(false); //Direction can only be 1 or -1
+          }
+        break;
+      default:
+        assert(false);
+        }
+      return octave_value_list ();
     }
 
-  return retval;
-}
-
-template <typename T>
-octave_value_list
-find_nonzero_elem_idx (const Sparse<T>& v, int nargout,
-                       octave_idx_type n_to_find, int direction)
-{
-  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
-
-
-  octave_idx_type nc = v.cols ();
-  octave_idx_type nr = v.rows ();
-  octave_idx_type nz = v.nnz ();
-
-  // Search in the default range.
-  octave_idx_type start_nc = -1;
-  octave_idx_type end_nc = -1;
-  octave_idx_type count;
-
-  // Search for the range to search
-  if (n_to_find < 0)
-    {
-      start_nc = 0;
-      end_nc = nc;
-      n_to_find = nz;
-      count = nz;
-    }
-  else if (direction > 0)
-    {
-      for (octave_idx_type j = 0; j < nc; j++)
-        {
-          OCTAVE_QUIT;
-          if (v.cidx (j) == 0 && v.cidx (j+1) != 0)
-            start_nc = j;
-          if (v.cidx (j+1) >= n_to_find)
-            {
-              end_nc = j + 1;
-              break;
-            }
-        }
-    }
-  else
-    {
-      for (octave_idx_type j = nc; j > 0; j--)
-        {
-          OCTAVE_QUIT;
-          if (v.cidx (j) == nz && v.cidx (j-1) != nz)
-            end_nc = j;
-          if (nz - v.cidx (j-1) >= n_to_find)
-            {
-              start_nc = j - 1;
-              break;
-            }
-        }
-    }
-
-  count = (n_to_find > v.cidx (end_nc) - v.cidx (start_nc) ?
-           v.cidx (end_nc) - v.cidx (start_nc) : n_to_find);
-
-  // If the original argument was a row vector, force a row vector of
-  // the overall indices to be returned.  But see below for scalar
-  // case...
-
-  octave_idx_type result_nr = count;
-  octave_idx_type result_nc = 1;
-
-  bool scalar_arg = false;
-
-  if (v.rows () == 1)
-    {
-      result_nr = 1;
-      result_nc = count;
-
-      scalar_arg = (v.columns () == 1);
-    }
-
-  Matrix idx (result_nr, result_nc);
-
-  Matrix i_idx (result_nr, result_nc);
-  Matrix j_idx (result_nr, result_nc);
-
-  Array<T> val (dim_vector (result_nr, result_nc));
-
-  if (count > 0)
-    {
-      // Search for elements to return.  Only search the region where
-      // there are elements to be found using the count that we want
-      // to find.
-      for (octave_idx_type j = start_nc, cx = 0; j < end_nc; j++)
-        for (octave_idx_type i = v.cidx (j); i < v.cidx (j+1); i++)
-          {
-            OCTAVE_QUIT;
-            if (direction < 0 && i < nz - count)
-              continue;
-            i_idx(cx) = static_cast<double> (v.ridx (i) + 1);
-            j_idx(cx) = static_cast<double> (j + 1);
-            idx(cx) = j * nr + v.ridx (i) + 1;
-            val(cx) = v.data(i);
-            cx++;
-            if (cx == count)
-              break;
-          }
-    }
-  else if (scalar_arg)
-    {
-      idx.resize (0, 0);
-
-      i_idx.resize (0, 0);
-      j_idx.resize (0, 0);
-
-      val.resize (dim_vector (0, 0));
-    }
-
-  switch (nargout)
-    {
-    case 0:
-    case 1:
-      retval(0) = idx;
-      break;
-
-    case 5:
-      retval(4) = nc;
-      // Fall through
-
-    case 4:
-      retval(3) = nr;
-      // Fall through
-
-    case 3:
-      retval(2) = val;
-      // Fall through!
-
-    case 2:
-      retval(1) = j_idx;
-      retval(0) = i_idx;
-      break;
-
-    default:
-      panic_impossible ();
-      break;
-    }
-
-  return retval;
-}
-
-octave_value_list
-find_nonzero_elem_idx (const PermMatrix& v, int nargout,
-                       octave_idx_type n_to_find, int direction)
-{
-  // There are far fewer special cases to handle for a PermMatrix.
-  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
-
-  octave_idx_type nc = v.cols ();
-  octave_idx_type start_nc, count;
-
-  // Determine the range to search.
-  if (n_to_find < 0 || n_to_find >= nc)
-    {
-      start_nc = 0;
-      count = nc;
-    }
-  else if (direction > 0)
-    {
-      start_nc = 0;
-      count = n_to_find;
-    }
-  else
-    {
-      start_nc = nc - n_to_find;
-      count = n_to_find;
-    }
-
-  bool scalar_arg = (v.rows () == 1 && v.cols () == 1);
-
-  Matrix idx (count, 1);
-  Matrix i_idx (count, 1);
-  Matrix j_idx (count, 1);
-  // Every value is 1.
-  Array<double> val (dim_vector (count, 1), 1.0);
-
-  if (count > 0)
-    {
-      const octave_idx_type* p = v.data ();
-	  for (octave_idx_type k = 0; k < count; k++)
-		{
-		  OCTAVE_QUIT;
-		  const octave_idx_type j = start_nc + k;
-		  const octave_idx_type i = p[j];
-		  i_idx(k) = static_cast<double> (1+i);
-		  j_idx(k) = static_cast<double> (1+j);
-		  idx(k) = j * nc + i + 1;
-		}
-    }
-  else if (scalar_arg)
-    {
-      // Same odd compatibility case as the other overrides.
-      idx.resize (0, 0);
-      i_idx.resize (0, 0);
-      j_idx.resize (0, 0);
-      val.resize (dim_vector (0, 0));
-    }
-
-  switch (nargout)
-    {
-    case 0:
-    case 1:
-      retval(0) = idx;
-      break;
-
-    case 5:
-      retval(4) = nc;
-      // Fall through
-
-    case 4:
-      retval(3) = nc;
-      // Fall through
-
-    case 3:
-      retval(2) = val;
-      // Fall through!
-
-    case 2:
-      retval(1) = j_idx;
-      retval(0) = i_idx;
-      break;
-
-    default:
-      panic_impossible ();
-      break;
-    }
-
-  return retval;
-}
-
 DEFUN (find, args, nargout,
        "-*- texinfo -*-\n\
 @deftypefn  {Built-in Function} {@var{idx} =} find (@var{x})\n\
@@ -435,7 +318,7 @@
           SparseBoolMatrix v = arg.sparse_bool_matrix_value ();
 
           if (! error_state)
-            retval = find_nonzero_elem_idx (v, nargout,
+            retval = find_templated (v, nargout,
                                             n_to_find, direction);
         }
       else if (nargout <= 1 && n_to_find == -1 && direction == 1)
@@ -449,7 +332,7 @@
           boolNDArray v = arg.bool_array_value ();
 
           if (! error_state)
-            retval = find_nonzero_elem_idx (v, nargout,
+            retval = find_templated (v, nargout,
                                             n_to_find, direction);
         }
     }
@@ -461,7 +344,7 @@
           INTT ## NDArray v = arg.INTT ## _array_value (); \
           \
           if (! error_state) \
-            retval = find_nonzero_elem_idx (v, nargout, \
+            retval = find_templated (v, nargout, \
                                             n_to_find, direction);\
         }
 
@@ -485,7 +368,7 @@
           SparseMatrix v = arg.sparse_matrix_value ();
 
           if (! error_state)
-            retval = find_nonzero_elem_idx (v, nargout,
+            retval = find_templated (v, nargout,
                                             n_to_find, direction);
         }
       else if (arg.is_complex_type ())
@@ -493,7 +376,7 @@
           SparseComplexMatrix v = arg.sparse_complex_matrix_value ();
 
           if (! error_state)
-            retval = find_nonzero_elem_idx (v, nargout,
+            retval = find_templated (v, nargout,
                                             n_to_find, direction);
         }
       else
@@ -504,14 +387,14 @@
       PermMatrix P = arg.perm_matrix_value ();
 
       if (! error_state)
-        retval = find_nonzero_elem_idx (P, nargout, n_to_find, direction);
+        retval = find_templated (P, nargout, n_to_find, direction);
     }
   else if (arg.is_string ())
     {
       charNDArray chnda = arg.char_array_value ();
 
       if (! error_state)
-        retval = find_nonzero_elem_idx (chnda, nargout, n_to_find, direction);
+        retval = find_templated (chnda, nargout, n_to_find, direction);
     }
   else if (arg.is_single_type ())
     {
@@ -520,7 +403,7 @@
           FloatNDArray nda = arg.float_array_value ();
 
           if (! error_state)
-            retval = find_nonzero_elem_idx (nda, nargout, n_to_find,
+            retval = find_templated (nda, nargout, n_to_find,
                                             direction);
         }
       else if (arg.is_complex_type ())
@@ -528,7 +411,7 @@
           FloatComplexNDArray cnda = arg.float_complex_array_value ();
 
           if (! error_state)
-            retval = find_nonzero_elem_idx (cnda, nargout, n_to_find,
+            retval = find_templated (cnda, nargout, n_to_find,
                                             direction);
         }
     }
@@ -537,14 +420,14 @@
       NDArray nda = arg.array_value ();
 
       if (! error_state)
-        retval = find_nonzero_elem_idx (nda, nargout, n_to_find, direction);
+        retval = find_templated (nda, nargout, n_to_find, direction);
     }
   else if (arg.is_complex_type ())
     {
       ComplexNDArray cnda = arg.complex_array_value ();
 
       if (! error_state)
-        retval = find_nonzero_elem_idx (cnda, nargout, n_to_find, direction);
+        retval = find_templated (cnda, nargout, n_to_find, direction);
     }
   else
     gripe_wrong_type_arg ("find", arg);
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/Array.cc
--- a/liboctave/array/Array.cc	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/Array.cc	Sun May 25 01:37:09 2014 -0600
@@ -2242,87 +2242,6 @@
 }
 
 template <class T>
-Array<octave_idx_type>
-Array<T>::find (octave_idx_type n, bool backward) const
-{
-  Array<octave_idx_type> retval;
-  const T *src = data ();
-  octave_idx_type nel = nelem ();
-  const T zero = T ();
-  if (n < 0 || n >= nel)
-    {
-      // We want all elements, which means we'll almost surely need
-      // to resize. So count first, then allocate array of exact size.
-      octave_idx_type cnt = 0;
-      for (octave_idx_type i = 0; i < nel; i++)
-        cnt += src[i] != zero;
-
-      retval.clear (cnt, 1);
-      octave_idx_type *dest = retval.fortran_vec ();
-      for (octave_idx_type i = 0; i < nel; i++)
-        if (src[i] != zero) *dest++ = i;
-    }
-  else
-    {
-      // We want a fixed max number of elements, usually small. So be
-      // optimistic, alloc the array in advance, and then resize if
-      // needed.
-      retval.clear (n, 1);
-      if (backward)
-        {
-          // Do the search as a series of successive single-element searches.
-          octave_idx_type k = 0;
-          octave_idx_type l = nel - 1;
-          for (; k < n; k++)
-            {
-              for (; l >= 0 && src[l] == zero; l--) ;
-              if (l >= 0)
-                retval(k) = l--;
-              else
-                break;
-            }
-          if (k < n)
-            retval.resize2 (k, 1);
-          octave_idx_type *rdata = retval.fortran_vec ();
-          std::reverse (rdata, rdata + k);
-        }
-      else
-        {
-          // Do the search as a series of successive single-element searches.
-          octave_idx_type k = 0;
-          octave_idx_type l = 0;
-          for (; k < n; k++)
-            {
-              for (; l != nel && src[l] == zero; l++) ;
-              if (l != nel)
-                retval(k) = l++;
-              else
-                break;
-            }
-          if (k < n)
-            retval.resize2 (k, 1);
-        }
-    }
-
-  // Fixup return dimensions, for Matlab compatibility.
-  // find (zeros (0,0)) -> zeros (0,0)
-  // find (zeros (1,0)) -> zeros (1,0)
-  // find (zeros (0,1)) -> zeros (0,1)
-  // find (zeros (0,X)) -> zeros (0,1)
-  // find (zeros (1,1)) -> zeros (0,0) !!!! WHY?
-  // find (zeros (0,1,0)) -> zeros (0,0)
-  // find (zeros (0,1,0,1)) -> zeros (0,0) etc
-
-  if ((numel () == 1 && retval.is_empty ())
-      || (rows () == 0 && dims ().numel (1) == 0))
-    retval.dimensions = dim_vector ();
-  else if (rows () == 1 && ndims () == 2)
-    retval.dimensions = dim_vector (1, retval.length ());
-
-  return retval;
-}
-
-template <class T>
 Array<T>
 Array<T>::nth_element (const idx_vector& n, int dim) const
 {
@@ -2516,9 +2435,6 @@
 template <> octave_idx_type \
 Array<T>::nnz (void) const\
 { return 0; } \
-template <> Array<octave_idx_type> \
-Array<T>::find (octave_idx_type, bool) const\
-{ return Array<octave_idx_type> (); } \
  \
 template <> Array<T>  \
 Array<T>::nth_element (const idx_vector&, int) const { return Array<T> (); } \
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/Array.h
--- a/liboctave/array/Array.h	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/Array.h	Sun May 25 01:37:09 2014 -0600
@@ -40,6 +40,7 @@
 #include "quit.h"
 #include "oct-mem.h"
 #include "oct-refcount.h"
+#include "nz-iterators-decl.h"
 
 // One dimensional array class.  Handles the reference counting for
 // all the derived classes.
@@ -123,6 +124,8 @@
 
   typedef T element_type;
 
+  typedef array_iterator<T> iter_type;
+
   typedef typename ref_param<T>::type crefT;
 
   typedef bool (*compare_fcn_type) (typename ref_param<T>::type,
@@ -611,11 +614,6 @@
   // Count nonzero elements.
   octave_idx_type nnz (void) const;
 
-  // Find indices of (at most n) nonzero elements. If n is specified, backward
-  // specifies search from backward.
-  Array<octave_idx_type> find (octave_idx_type n = -1,
-                               bool backward = false) const;
-
   // Returns the n-th element in increasing order, using the same ordering as
   // used for sort. n can either be a scalar index or a contiguous range.
   Array<T> nth_element (const idx_vector& n, int dim = 0) const;
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/DiagArray2.h
--- a/liboctave/array/DiagArray2.h	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/DiagArray2.h	Sun May 25 01:37:09 2014 -0600
@@ -30,6 +30,7 @@
 #include <cstdlib>
 
 #include "Array.h"
+#include "nz-iterators-decl.h"
 
 // Array<T> is inherited privately so that some methods, like index, don't
 // produce unexpected results.
@@ -45,6 +46,8 @@
 
   using typename Array<T>::element_type;
 
+  typedef diag_iterator<T> iter_type;
+
   DiagArray2 (void)
     : Array<T> (), d1 (0), d2 (0) { }
 
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/PermMatrix.h
--- a/liboctave/array/PermMatrix.h	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/PermMatrix.h	Sun May 25 01:37:09 2014 -0600
@@ -25,6 +25,7 @@
 
 #include "Array.h"
 #include "mx-defs.h"
+#include "nz-iterators-decl.h"
 
 // Array<T> is inherited privately so that some methods, like index, don't
 // produce unexpected results.
@@ -37,6 +38,9 @@
 
 public:
 
+  typedef bool element_type;
+  typedef perm_iterator iter_type;
+
   PermMatrix (void) : Array<octave_idx_type> () { }
 
   PermMatrix (octave_idx_type n);
@@ -64,6 +68,7 @@
   { return perm_length (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
   octave_idx_type numel (void) const { return nelem (); }
+  octave_idx_type nnz (void) const {return perm_length (); }
 
   size_t byte_size (void) const
   { return Array<octave_idx_type>::byte_size (); }
@@ -105,6 +110,9 @@
   friend OCTAVE_API PermMatrix operator *(const PermMatrix& a,
                                           const PermMatrix& b);
 
+  octave_idx_type perm_elem(octave_idx_type i) const
+    { return Array<octave_idx_type>::elem (i);}
+
   const octave_idx_type *data (void) const
   { return Array<octave_idx_type>::data (); }
 
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/Sparse.h
--- a/liboctave/array/Sparse.h	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/Sparse.h	Sun May 25 01:37:09 2014 -0600
@@ -36,6 +36,7 @@
 #include "dim-vector.h"
 #include "lo-error.h"
 #include "lo-utils.h"
+#include "nz-iterators-decl.h"
 
 #include "oct-sort.h"
 #include "oct-mem.h"
@@ -54,6 +55,8 @@
 
   typedef T element_type;
 
+  typedef sparse_iterator<T> iter_type;
+
 protected:
   //--------------------------------------------------------------------
   // The real representation of all Sparse arrays.
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/dim-vector.h
--- a/liboctave/array/dim-vector.h	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/dim-vector.h	Sun May 25 01:37:09 2014 -0600
@@ -339,7 +339,9 @@
 
   // Return the number of elements that a matrix with this dimension
   // vector would have, NOT the number of dimensions (elements in the
-  // dimension vector).
+  // dimension vector).  If given an argument n, returns the number of
+  // elements an array with dimensions dim[n..end] would have
+  // (ex. numel(1) returns the number of columns)
 
   octave_idx_type numel (int n = 0) const
   {
diff -r 330cb3efc6ba -r df0221588028 liboctave/array/idx-vector.cc
--- a/liboctave/array/idx-vector.cc	Sat May 24 14:14:21 2014 -0600
+++ b/liboctave/array/idx-vector.cc	Sun May 25 01:37:09 2014 -0600
@@ -39,6 +39,7 @@
 #include "oct-locbuf.h"
 #include "lo-error.h"
 #include "lo-mappers.h"
+#include "find.h"
 
 static void
 gripe_invalid_range (void)
@@ -753,7 +754,7 @@
 idx_vector::idx_mask_rep::as_array (void)
 {
   if (aowner)
-    return aowner->find ().reshape (orig_dims);
+    return find (*aowner).reshape (orig_dims);
   else
     {
       Array<bool> retval (orig_dims);
diff -r 330cb3efc6ba -r df0221588028 liboctave/util/find.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/util/find.h	Sun May 25 01:37:09 2014 -0600
@@ -0,0 +1,151 @@
+/*
+ * find.h
+ *
+ *  Created on: May 20, 2014
+ *      Author: dspies
+ */
+
+#ifndef FIND_H_
+#define FIND_H_
+
+/*
+
+ Copyright (C) 1996-2013 John W. Eaton
+
+ This file is part of Octave.
+
+ Octave is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 3 of the License, or (at your
+ option) any later version.
+
+ Octave is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Octave; see the file COPYING.  If not, see
+ <http://www.gnu.org/licenses/>.
+
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "quit.h"
+
+#include "interp-idx.h"
+#include "min-with-nnz.h"
+#include "nz-iterators.h"
+#include "step-dir.h"
+
+/*
+ * This struct mostly matches the signature of res_type_of (see find.cc)
+ * except without a get_list function.  Internal calls to "find" from
+ * liboctave should call the "find" method which in turn constructs one
+ * of these as a value to be filled.
+ */
+struct array_res_1 {
+  array_res_1()
+    {
+    }
+  array_res_1(const dim_vector& nnz) :
+      res (nnz)
+    {
+    }
+  Array<octave_idx_type> res;
+  template<typename T>
+  void add(octave_idx_type place, array_iterator<T> iter)
+    {
+      res.xelem (place) = iter.flat_idx ();
+    }
+};
+
+/*
+ * This is the generic "find" method for all matrix types.  All calls to
+ * "find" get routed here one way or another as this method contains logic
+ * for determining the proper output dimensions and handling forward/reversed,
+ * n_to_find etc.
+ * It would be unwise to try and duplicate this logic elsewhere as it's fairly
+ * nuanced and unintuitive (but Matlab compatible).
+ */
+template<int direction, typename S, typename R>
+void find_to_iter(dir_constant<direction> dirc, const S& v,
+    octave_idx_type n_to_find, R& resvec)
+  {
+    typedef typename S::iter_type iter_t;
+
+    octave_idx_type numres;
+    if (n_to_find == -1)
+      numres = v.nnz ();
+    else
+      numres = min_with_nnz (v, n_to_find);
+
+    const dim_vector& dv = v.dims ();
+    octave_idx_type col_len = dv.elem (0);
+    octave_idx_type cols = dv.numel (1);
+
+    // Fixup return dimensions, for Matlab compatibility.
+    // find (zeros (0,0)) -> zeros (0,0)
+    // find (zeros (1,0)) -> zeros (1,0)
+    // find (zeros (0,1)) -> zeros (0,1)
+    // find (zeros (0,X)) -> zeros (0,1)
+    // find (zeros (1,1)) -> zeros (0,0) !!!! WHY?
+    // find (zeros (0,1,0)) -> zeros (0,0)
+    // find (zeros (0,1,0,1)) -> zeros (0,0) etc
+
+    dim_vector res_dims;
+
+    if ((cols == 1 && col_len == 1 && numres == 0)
+        || (col_len == 0 && cols == 0))
+      res_dims = dim_vector (0, 0);
+    else if (col_len == 1)
+      res_dims = dim_vector (1, numres);
+    else
+      res_dims = dim_vector (numres, 1);
+
+    iter_t iterator (v);
+    resvec = R (res_dims);
+    octave_idx_type count = choose_dir (dirc, numres - 1, 0);
+    octave_idx_type end = choose_dir (dirc, 0, numres);
+
+    for (iterator.begin (dirc);
+        !iterator.finished (dirc) && cmp_ge_lt (dirc, count, end);
+        iterator.step (dirc), count += direction)
+      {
+        resvec.add (count, iterator);
+      }
+  }
+
+/*
+ * The method call used for internal liboctave calls to find.  It constructs an
+ * array_res_1 and fills that rather than a res_type_of (see "find.cc"). Note
+ * that this method only supports full arrays because octave_idx_type will
+ * generally overflow if used to represent the nonzero indices of sparse arrays.
+ * Instead consider using the sparse-iterator class from nz-iterators.h to
+ * avoid this problem
+ */
+template<typename T>
+Array<octave_idx_type> find(const Array<T>& v, octave_idx_type n_to_find = -1,
+    int dir = 1)
+  {
+    array_res_1 resvec;
+    switch(dir){
+    case 1:
+      dir_constant<1> dirc_forward;
+      find_to_iter (dirc_forward, v, n_to_find, resvec);
+      break;
+    case -1:
+      dir_constant<-1> dirc_back;
+      find_to_iter (dirc_back, v, n_to_find, resvec);
+      break;
+    default:
+      liboctave_fatal("Unknown direction: %d; Must be 1 or -1", dir);
+    }
+
+    return resvec.res;
+  }
+
+#endif /* FIND_H_ */
diff -r 330cb3efc6ba -r df0221588028 liboctave/util/interp-idx.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/util/interp-idx.h	Sun May 25 01:37:09 2014 -0600
@@ -0,0 +1,29 @@
+/*
+ * interp_idx.h
+ *
+ *  Created on: May 22, 2014
+ *      Author: dspies
+ */
+
+#ifndef INTERP_IDX_H_
+#define INTERP_IDX_H_
+
+// Simple method for converting between C++ octave_idx_type and
+// Octave data-types (convert to double and add one).
+inline double to_interp_idx(octave_idx_type idx)
+  {
+    return idx + 1.L;
+  }
+
+// Simple method for taking a row-column pair together with the matrix
+// height and returning the corresponding index as an octave data-type
+// (note that for large heights, there's a risk of losing precision.
+// This method will not overflow or throw a bad alloc, it will simply
+// choose the nearest possible double-value to the proper index).
+inline double to_interp_idx(octave_idx_type row, octave_idx_type col,
+    octave_idx_type height)
+  {
+    return col * static_cast<double> (height) + row + 1;
+  }
+
+#endif /* INTERP_IDX_H_ */
diff -r 330cb3efc6ba -r df0221588028 liboctave/util/min-with-nnz.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/util/min-with-nnz.h	Sun May 25 01:37:09 2014 -0600
@@ -0,0 +1,48 @@
+/*
+ * minWithNNZ.h
+ *
+ *  Created on: May 20, 2014
+ *      Author: dspies
+ */
+
+#ifndef MINWITHNNZ_H_
+#define MINWITHNNZ_H_
+
+#include "Sparse.h"
+#include "step-dir.h"
+
+/*
+ * Generic efficient methods for finding the equivalent of
+ * min(arr.nnz(), minwith).  In case arr.nnz() happens to run
+ * in O(arr.numel()) time, this operation will short-circuit and return
+ * as soon as it counts minwith nonzero elements.
+ * But if arr has a constant-time nnz() operation, then this will use
+ * that instead and run in constant-time.
+ */
+
+template<typename M>
+octave_idx_type min_with_nnz(M arr, octave_idx_type minwith)
+  {
+    octave_idx_type count;
+    typename M::iter_type iter (arr);
+    dir_constant<1> dirc;
+    for (iter.begin (dirc), count = 0; !iter.finished (dirc) && count < minwith;
+        iter.step (dirc), ++count)
+      ;
+    return count;
+  }
+
+template<typename T>
+octave_idx_type min_with_nnz(Sparse<T> arr, octave_idx_type minwith)
+  {
+    return std::min (arr.nnz (), minwith);
+  }
+
+/*
+ * FIXME: I want to add one for PermMatrix, but for some reason C++ doesn't
+ * like it when you overload templated functions with non-templated ones.
+ * Except that PermMatrix is not templated.  Anyone know how to handle that?
+ * Should I add a dummy template parameter? (that seems like quite an ugly hack)
+ */
+
+#endif /* MINWITHNNZ_H_ */
diff -r 330cb3efc6ba -r df0221588028 liboctave/util/nz-iterators-decl.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/util/nz-iterators-decl.h	Sun May 25 01:37:09 2014 -0600
@@ -0,0 +1,32 @@
+/*
+ * nz_iterators_decl.h
+ *
+ *  Created on: May 18, 2014
+ *      Author: dspies
+ */
+
+#ifndef NZ_ITERATORS_DECL_H_
+#define NZ_ITERATORS_DECL_H_
+
+/*
+ * Contains just the class declarations for nz-iterators.h
+ * (this is all the base matrix classes Array, Sparse, DiagArray2,
+ * and PermMatrix need to compile, so changes to nz-iterators.h
+ * will not cause those classes to need to be recompiled unless the
+ * class declarations are changed in some way)
+ *
+ * If you wish to use these iterators, you must #include nz-iterators.h
+ */
+
+template<typename T>
+class array_iterator;
+
+template<typename T>
+class sparse_iterator;
+
+template<typename T>
+class diag_iterator;
+
+class perm_iterator;
+
+#endif /* NZ_ITERATORS_DECL_H_ */
diff -r 330cb3efc6ba -r df0221588028 liboctave/util/nz-iterators.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/util/nz-iterators.h	Sun May 25 01:37:09 2014 -0600
@@ -0,0 +1,360 @@
+/*
+ * nz_iterators.h
+ *
+ *  Created on: May 18, 2014
+ *      Author: dspies
+ */
+
+#ifndef NZ_ITERATORS_H_
+#define NZ_ITERATORS_H_
+
+#include "interp-idx.h"
+#include "nz-iterators-decl.h"
+#include "oct-inttypes.h"
+#include "Array.h"
+#include "DiagArray2.h"
+#include "PermMatrix.h"
+#include "Sparse.h"
+#include "step-dir.h"
+
+/*
+ * This file contains generic column-major iterators over
+ * the nonzero elements of any array or matrix.  If you have a matrix mat
+ * of type M, you can construct the proper iterator type using
+ * M::iter_type iter(mat) and iter will iterate efficiently (forwards
+ * or backwards) over the nonzero elements of mat.
+ *
+ * The parameter T is the element-type except for PermMatrix where
+ * the element type is always bool
+ *
+ * Use a dir_constant to indicate which direction you intend to iterate.
+ * (see step-dir.h).  begin() resets the iterator to the beginning or
+ * end of the matrix (for dir_constant<1> and <-1> respectively).
+ * finished(dirc) indicates whether the iterators has finished traversing
+ * the nonzero elements.  step(dirc) steps from one element to the next.
+ *
+ * You can, for instance, use a for-loop as follows:
+ *
+ typedef M::iter_type iter_t;
+ typedef M::element_type T;
+
+ iter_t iter(mat);
+ dir_constant<1> dirc;
+
+ for(iter.begin (dirc); !iter.finished (dirc); iter.step (dirc))
+ {
+   octave_idx_type row = iter.row();
+   octave_idx_type col = iter.col();
+   double doub_index = iter.interp_index ();
+   T elem = iter.data();
+   // ... Do something with these
+ }
+
+ * Note that array_iter for indexing over full matrices also includes
+ * a iter.flat_index () method which returns an octave_idx_type.
+ *
+ * The other iterators to not have a flat_index() method because they
+ * risk overflowing octave_idx_type.  It is recommended you take care
+ * to implement your function in a way that accounts for this problem.
+ *
+ * FIXME: I'd like to add in these
+ * default no-parameter versions of
+ * begin() and step() to each of the
+ * classes.  But the C++ compiler complains
+ * because apparently I'm not allowed to overload
+ * templated methods with non-templated ones.  Any
+ * ideas for work-arounds?
+
+ #define INCLUDE_DEFAULT_STEPS \
+  void begin() \
+    { \
+      dir_constant<1> dirc; \
+      begin (dirc); \
+    } \
+  void step() \
+    { \
+      dir_constant<1> dirc; \
+      step (dirc); \
+    } \
+  bool finished() const \
+    { \
+      dir_constant<1> dirc; \
+      return finished (dirc); \
+    }
+ */
+
+// A generic method for checking if some element of a matrix with
+// element type T is zero.
+template<typename T>
+bool is_zero(T t)
+  {
+    return t == static_cast<T> (0);
+  }
+
+/*
+ * An iterator over full arrays.  When the number of dimensions exceeds
+ * 2, calls to iter.col() may exceed mat.cols() up to mat.dims().numel(1)
+ *
+ * This mimics the behavior of the "find" method (both in Octave and Matlab)
+ * on many-dimensional matrices.
+ */
+template<typename T>
+class array_iterator
+  {
+private:
+  const Array<T>& mat;
+
+  //Actual total number of columns = mat.dims().numel(1)
+  //can be different from length of row dimension
+  const octave_idx_type totcols;
+
+  octave_idx_type coli;
+  octave_idx_type rowj;
+  octave_idx_type my_idx;
+
+  template<int dir>
+  void step_once(dir_constant<dir> dirc)
+    {
+      my_idx += dir;
+      rowj += dir;
+      octave_idx_type cmpwith = choose_dir (dirc, 0, mat.rows ());
+      if (cmp_lt_ge (dirc, rowj, cmpwith))
+        {
+          rowj = choose_dir (dirc, mat.rows () - 1, 0);
+          coli += dir;
+        }
+    }
+
+public:
+  array_iterator(const Array<T>& arg_mat) :
+      mat (arg_mat), totcols (arg_mat.dims ().numel (1))
+    {
+      dir_constant<1> dirc;
+      begin (dirc);
+    }
+
+  template<int dir>
+  void begin(dir_constant<dir> dirc)
+    {
+      coli = choose_dir (dirc, totcols, 0);
+      rowj = choose_dir (dirc, 0, -1);
+      my_idx = choose_dir (dirc, mat.numel (), -1);
+      step (dirc);
+    }
+
+  octave_idx_type col() const
+    {
+      return coli;
+    }
+  octave_idx_type row() const
+    {
+      return rowj;
+    }
+  double interp_idx() const
+    {
+      return to_interp_idx (my_idx);
+    }
+  octave_idx_type flat_idx() const
+    {
+      return my_idx;
+    }
+  T data() const
+    {
+      return mat.elem (my_idx);
+    }
+
+  template<int dir>
+  void step(dir_constant<dir> dirc)
+    {
+      do
+        {
+          step_once (dirc);
+        }
+      while (!finished (dirc) && is_zero (data ()));
+    }
+  template<int dir>
+  bool finished(dir_constant<dir> dirc) const
+    {
+      octave_idx_type cmpwith = choose_dir (dirc, 0, totcols);
+      return cmp_lt_ge (dirc, coli, cmpwith);
+    }
+  };
+
+template<typename T>
+class sparse_iterator
+  {
+private:
+  const Sparse<T>& mat;
+  octave_idx_type coli;
+  octave_idx_type my_idx;
+
+  template<int dir>
+  octave_idx_type test_col_idx(dir_constant<dir>)
+    {
+      return mat.cidx (coli + std::max (dir, 0));
+    }
+
+  template<int dir>
+  void adjust_col(dir_constant<dir> dirc)
+    {
+      while (!finished (dirc) && cmp_lt_ge (dirc, my_idx, test_col_idx (dirc)))
+        coli += dir;
+    }
+
+public:
+  sparse_iterator(const Sparse<T>& arg_mat) :
+      mat (arg_mat)
+    {
+      dir_constant<1> dirc;
+      begin (dirc);
+    }
+
+  template<int dir>
+  void begin(dir_constant<dir> dirc)
+    {
+      coli = choose_dir (dirc, mat.cols () - 1, 0);
+      my_idx = choose_dir (dirc, mat.nnz () - 1, 0);
+      adjust_col (dirc);
+    }
+
+  double interp_idx() const
+    {
+      return to_interp_idx (row (), col (), mat.rows ());
+    }
+  octave_idx_type col() const
+    {
+      return coli;
+    }
+  octave_idx_type row() const
+    {
+      return mat.ridx (my_idx);
+    }
+  T data() const
+    {
+      return mat.data (my_idx);
+    }
+  template<int dir>
+  void step(dir_constant<dir> dirc)
+    {
+      my_idx += dir;
+      adjust_col (dirc);
+    }
+  template<int dir>
+  bool finished(dir_constant<dir> dirc) const
+    {
+      octave_idx_type cmpwith = choose_dir (dirc, 0, mat.cols ());
+      return cmp_lt_ge (dirc, coli, cmpwith);
+    }
+  };
+
+template<typename T>
+class diag_iterator
+  {
+private:
+  const DiagArray2<T>& mat;
+  octave_idx_type my_idx;
+
+public:
+  diag_iterator(const DiagArray2<T>& arg_mat) :
+      mat (arg_mat)
+    {
+      dir_constant<1> dirc;
+      begin (dirc);
+    }
+
+  void begin()
+    {
+      dir_constant<1> dirc;
+      begin (dirc);
+    }
+
+  template<int dir>
+  void begin(dir_constant<dir> dirc)
+    {
+      my_idx = choose_dir (dirc, mat.diag_length (), -1);
+      step (dirc);
+    }
+
+  double interp_idx() const
+    {
+      return to_interp_idx (row (), col (), mat.rows ());
+    }
+  octave_idx_type col() const
+    {
+      return my_idx;
+    }
+  octave_idx_type row() const
+    {
+      return my_idx;
+    }
+  T data() const
+    {
+      return mat.dgelem (my_idx);
+    }
+  template<int dir>
+  void step(dir_constant<dir> dirc)
+    {
+      do
+        {
+          my_idx += dir;
+        }
+      while (!finished (dirc) && is_zero (data ()));
+    }
+  template<int dir>
+  bool finished(dir_constant<dir> dirc) const
+    {
+      octave_idx_type cmpwith = choose_dir (dirc, 0, mat.diag_length ());
+      return cmp_lt_ge (dirc, my_idx, cmpwith);
+    }
+  };
+
+class perm_iterator
+  {
+private:
+  const PermMatrix& mat;
+  octave_idx_type my_idx;
+
+public:
+  perm_iterator(const PermMatrix& arg_mat) :
+      mat (arg_mat)
+    {
+      dir_constant<1> dirc;
+      begin (dirc);
+    }
+
+  template<int dir>
+  void begin(dir_constant<dir> dirc)
+    {
+      my_idx = choose_dir (dirc, mat.cols () - 1, 0);
+    }
+
+  octave_idx_type interp_idx() const
+    {
+      return to_interp_idx (row (), col (), mat.rows ());
+    }
+  octave_idx_type col() const
+    {
+      return my_idx;
+    }
+  octave_idx_type row() const
+    {
+      return mat.perm_elem (my_idx);
+    }
+  bool data() const
+    {
+      return true;
+    }
+  template<int dir>
+  void step(dir_constant<dir>)
+    {
+      my_idx += dir;
+    }
+  template<int dir>
+  bool finished(dir_constant<dir> dirc) const
+    {
+      octave_idx_type cmpwith = choose_dir (dirc, 0, mat.rows ());
+      return cmp_lt_ge (dirc, my_idx, cmpwith);
+    }
+  };
+
+#endif /* NZ_ITERATORS_H_ */
diff -r 330cb3efc6ba -r df0221588028 liboctave/util/step-dir.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/util/step-dir.h	Sun May 25 01:37:09 2014 -0600
@@ -0,0 +1,147 @@
+/*
+ * step_dir.h
+ *
+ *  Created on: May 18, 2014
+ *      Author: dspies
+ */
+
+#ifndef STEP_DIR_H_
+#define STEP_DIR_H_
+
+#include <cassert>
+
+/*
+ * Simple generic parameterized functions for stepping "forward" or "backward"
+ * 1 indicates "forward"
+ * -1 indicates "backward"
+ *
+ * Do not call these with any other template argument besides 1 or -1
+ */
+
+/*
+ * An empty struct parameterized with a direction.  By using this type as an
+ * argument to all the functions in this file, the compiler can automatically
+ * to infer the template parameter.  This avoids those awkward instances where
+ * C++ doesn't allow you to specify the template arguments of the method
+ * you want to call explicitly.  In those cases, type inference is the only way
+ * to say which one you want.
+ */
+template<int dir>
+struct dir_constant{};
+
+/*
+ * Returns either ifn1 or if1 depending on whether dir is -1 or 1 respectively
+ */
+template<int dir>
+inline octave_idx_type choose_dir(dir_constant<dir>, octave_idx_type ifn1, octave_idx_type ifp1)
+  {
+    if (dir == -1)
+      return ifn1;
+    else if (dir == 1)
+      return ifp1;
+    else
+      assert(false);
+  }
+
+/*
+ * The following eight functions compare a with b according to the
+ * corresponding comparison operation for the passed dir_constant.
+ *
+ * For instance, cmp_lt_gt<-1> returns a < b
+ * and
+ * cmp_lt_gt<1> returns a > b
+ *
+ * This is handy when the structure of a method is identical for
+ * "forwards" and "backwards" cases except for the comparison method to
+ * be used.
+ */
+
+template<int dir>
+inline bool cmp_lt_gt(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a < b;
+    else if (dir == 1)
+      return a > b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_lt_ge(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a < b;
+    else if (dir == 1)
+      return a >= b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_le_gt(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a <= b;
+    else if (dir == 1)
+      return a > b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_le_ge(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a <= b;
+    else if (dir == 1)
+      return a >= b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_gt_lt(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a > b;
+    else if (dir == 1)
+      return a < b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_gt_le(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a > b;
+    else if (dir == 1)
+      return a <= b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_ge_lt(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a >= b;
+    else if (dir == 1)
+      return a < b;
+    else
+      assert(false);
+  }
+
+template<int dir>
+inline bool cmp_ge_le(dir_constant<dir>, octave_idx_type a, octave_idx_type b)
+  {
+    if (dir == -1)
+      return a >= b;
+    else if (dir == 1)
+      return a <= b;
+    else
+      assert(false);
+  }
+
+#endif /* STEP_DIR_H_ */
diff -r 330cb3efc6ba -r df0221588028 test/build-sparse-tests.sh
--- a/test/build-sparse-tests.sh	Sat May 24 14:14:21 2014 -0600
+++ b/test/build-sparse-tests.sh	Sun May 25 01:37:09 2014 -0600
@@ -569,7 +569,8 @@
 %!assert (as==af)
 %!assert (af==as)
 %!test
-%! [ii,jj,vv,nr,nc] = find (as);
+%! [ii,jj,vv] = find (as);
+%! [nr,nc] = size (as);
 %! assert (af, full (sparse (ii,jj,vv,nr,nc)));
 %!assert (nnz (as), sum (af(:)!=0))
 %!assert (nnz (as), nnz (af))
@@ -598,7 +599,8 @@
 %! x = sparse (i,j,v,m,n);
 %! assert (x, as);
 %!test
-%! [i,j,v,m,n] = find (as);
+%! [i,j,v] = find (as);
+%! [m,n] = size (as);
 %! x = sparse (i,j,v,m,n);
 %! assert (x, as);
 %!assert (issparse (horzcat (as,as)));
# HG changeset patch
# User David Spies <dnspies@gmail.com>
# Date 1401042017 21600
#      Sun May 25 12:20:17 2014 -0600
# Node ID 3c3ca5ddd4293f6ebeb93db9e583b3276b10e333
# Parent  df0221588028f7c21875245a4483343c98ae1730
Added unit test for bug #42421 to show since the previous commit fixes it

diff -r df0221588028 -r 3c3ca5ddd429 libinterp/corefcn/find.cc
--- a/libinterp/corefcn/find.cc	Sun May 25 01:37:09 2014 -0600
+++ b/libinterp/corefcn/find.cc	Sun May 25 12:20:17 2014 -0600
@@ -483,5 +483,11 @@
 %!assert (find ([2 0 1 0 5 0], Inf), [1, 3, 5])
 %!assert (find ([2 0 1 0 5 0], Inf, "last"), [1, 3, 5])
 
+%!test
+%! x = sparse(100000, 30000);
+%! x(end, end) = 1;
+%! i = find(x);
+%! assert (i == 3e09);
+
 %!error find ()
 */
