diff -r 89a8a90ac850 libinterp/corefcn/find.cc
--- a/libinterp/corefcn/find.cc	Tue Mar 18 11:06:13 2014 +0100
+++ b/libinterp/corefcn/find.cc	Wed Mar 19 18:27:37 2014 -0600
@@ -31,306 +31,329 @@
 #include "gripes.h"
 #include "oct-obj.h"
 
+template<typename T, T v>
+struct integral_constant {};
+
+template<typename T>
+inline T range_start(T start, T, integral_constant<int, 1>)
+{ return start; }
+
+template<typename T>
+inline T range_start(T, T end, integral_constant<int, -1>)
+{ return end - 1; }
+
+template<typename T>
+inline bool finished(T i, T, T end, integral_constant<int, 1>)
+{ return i >= end; }
+
+template<typename T>
+inline bool finished(T i, T start, T, integral_constant<int, -1>)
+{ return i < start; }
+
+template<typename T>
+class array_nz_iter {
+private:
+  const Array<T>& mat;
+  octave_idx_type idx;
+  octave_idx_type row;
+  octave_idx_type col;
+  template<int direction>
+  void take_step(integral_constant<int,direction> dir)
+    {
+      idx += direction;
+      row += direction;
+      if (finished (row, 0, mat.rows (), dir))
+        {
+          row = range_start (0, mat.rows (), dir);
+          col += direction;
+        }
+    }
+  template<int direction>
+  void find_nz(integral_constant<int,direction> dir)
+    {
+      while(!iter_finished (dir) && mat(idx) == T(0))
+        {
+          take_step(dir);
+        }
+    }
+public:
+  template <int direction>
+  array_nz_iter(const Array<T>& arg_mat, integral_constant<int, direction> dir):mat(arg_mat)
+    {
+      row = range_start (0, mat.rows (), dir);
+      //The 'int' argument to numel is undocumented,
+      //but it seems like this is what it does. Can I
+      //use it safely?
+      col = range_start (0, mat.dims().numel(1), dir);
+      idx = range_start (0, mat.numel (), dir);
+      find_nz (dir);
+    }
+  template <int direction>
+  bool iter_finished(integral_constant<int, direction> dir) const
+    {
+      return finished<octave_idx_type>(idx, 0, mat.numel (), dir);
+    }
+  template <int direction>
+  void next(integral_constant<int, direction> dir)
+    {
+      take_step (dir);
+      find_nz (dir);
+    }
+  T operator*() const
+    {
+      return mat(idx);
+    }
+  octave_idx_type get_row() const { return row; }
+  octave_idx_type get_col() const { return col; }
+  double get_idx() const { return idx; }
+};
+
+template<typename T>
+class sparse_nz_iter {
+private:
+  const Sparse<T>& mat;
+  octave_idx_type nz_idx;
+  octave_idx_type col;
+  template <int direction>
+  void catchup_col(integral_constant<int, direction> dir)
+    {
+      while (!iter_finished (dir)
+          && finished (nz_idx, mat.cidx()[col], mat.cidx()[col + 1], dir))
+        col += direction;
+    }
+public:
+  template<int direction>
+  sparse_nz_iter(const Sparse<T>& arg_mat,
+      integral_constant<int, direction> dir) :
+      mat(arg_mat)
+    {
+      nz_idx = range_start(0, mat.nnz (), dir);
+      col = range_start(0, mat.cols (), dir);
+      catchup_col (dir);
+    }
+  template<int direction>
+  bool iter_finished(integral_constant<int, direction> dir) const
+    {
+      return finished<octave_idx_type>(nz_idx, 0, mat.nnz (), dir);
+    }
+  template<int direction>
+  void next(integral_constant<int, direction> dir)
+    {
+      nz_idx += direction;
+      catchup_col (dir);
+    }
+  T operator*() const
+    {
+      return mat.data()[nz_idx];
+    }
+  octave_idx_type get_row() const { return mat.ridx()[nz_idx]; }
+  octave_idx_type get_col() const { return col; }
+  double get_idx() const { return get_col () * double(mat.rows ()) + get_row (); }
+};
+
+class perm_nz_iter
+{
+private:
+  const PermMatrix& mat;
+  octave_idx_type col;
+public:
+  template<int direction>
+  perm_nz_iter(const PermMatrix& arg_mat, integral_constant<int, direction> dir) :
+      mat(arg_mat)
+    {
+      assert(mat.is_col_perm ());
+      col = range_start(0, mat.cols (), dir);
+    }
+  template<int direction>
+  bool iter_finished(integral_constant<int, direction> dir) const
+    {
+      return finished<octave_idx_type>(col, 0, mat.cols (), dir);
+    }
+  template<int direction>
+  void next(integral_constant<int, direction>)
+    {
+      col += direction;
+    }
+  double operator*() const
+    {
+      return 1;
+    }
+  octave_idx_type get_row() const { return mat.fortran_vec()[col]; }
+  octave_idx_type get_col() const { return col; }
+  double get_idx() const { return get_col() * double(mat.rows ()) + get_row (); }
+};
+
+//Surprisingly, this isn't already in the
+//code somewhere.  Perhaps it should be moved
+//to some utility file somewhere.
+template<typename T>
+inline T min(T a, T b)
+{
+  return a < b ? a : b;
+}
+
+template<typename ItType, typename M>
+octave_idx_type
+min_with_nnz(const M& nda, octave_idx_type n_to_find)
+{
+  if(n_to_find == 0)
+    return 0;
+  octave_idx_type i = 0;
+  integral_constant<int, 1> dir;
+  for (ItType iter(nda, dir); !iter.iter_finished (dir); iter.next (dir))
+    {
+      ++i;
+      if (n_to_find >= 0 && i >= n_to_find)
+        return n_to_find;
+    }
+  return i;
+}
+
+template<typename ItType, typename T>
+octave_idx_type
+min_with_nnz(const Sparse<T>& nda, octave_idx_type n_to_find)
+{
+  if(n_to_find < 0)
+    return nda.nnz ();
+  else
+    return min<octave_idx_type> (nda.nnz (), n_to_find);
+}
+
+template<typename ItType>
+octave_idx_type
+min_with_nnz(const PermMatrix& nda, octave_idx_type n_to_find)
+{
+  if(n_to_find < 0)
+    return nda.cols ();
+  else
+    return min<octave_idx_type> (nda.cols (), n_to_find);
+}
+
 // 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.
-
-template <typename T>
+template <int nargout, int direction, typename T, typename M, typename ItType>
 octave_value_list
-find_nonzero_elem_idx (const Array<T>& nda, int nargout,
-                       octave_idx_type n_to_find, int direction)
+find_nonzero_elem_idx_iterized(const M& nda, octave_idx_type n_to_find)
 {
-  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
-
-  Array<octave_idx_type> idx;
-  if (n_to_find >= 0)
-    idx = nda.find (n_to_find, direction == -1);
-  else
-    idx = nda.find ();
-
-  // The maximum element is always at the end.
-  octave_idx_type iext = idx.is_empty () ? 0 : idx.xelem (idx.numel () - 1) + 1;
-
-  switch (nargout)
+  n_to_find = min_with_nnz<ItType> (nda, n_to_find);
+  dim_vector res_dims;
+  //Absurd size rules for Matlab compatibility
+  if (n_to_find == 0
+      && ((nda.rows() == 0 && nda.cols() == 0) || nda.dims().all_ones()))
     {
-    default:
-    case 3:
-      retval(2) = Array<T> (nda.index (idx_vector (idx)));
-      // Fall through!
-
-    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++)
-          {
-            jdx.xelem (i) = idx.xelem (i) / nr;
-            idx.xelem (i) %= nr;
-          }
-        iext = -1;
-        retval(1) = idx_vector (jdx, -1);
-      }
-      // Fall through!
-
-    case 1:
-    case 0:
-      retval(0) = idx_vector (idx, iext);
-      break;
+      res_dims = dim_vector(0, 0);
     }
-
-  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)
+  else if(nda.rows() == 1)
     {
-      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;
-            }
-        }
+      res_dims = dim_vector(1, n_to_find);
     }
   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;
-            }
-        }
+      res_dims = dim_vector(n_to_find, 1);
     }
 
-  count = (n_to_find > v.cidx (end_nc) - v.cidx (start_nc) ?
-           v.cidx (end_nc) - v.cidx (start_nc) : n_to_find);
+  static const int retsize = min<int>(nargout, 2);
+  Array<double> retval[retsize];
+  for(int i = 0; i < retsize; ++i)
+    {
+      retval[i] = Array<double> (res_dims);
+    }
+  Array<T> vals;
+  if(nargout >= 3)
+    vals = Array<T> (res_dims);
 
-  // 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...
+  integral_constant<int,direction> dir;
+  octave_idx_type i = range_start (0, n_to_find, dir);
 
-  octave_idx_type result_nr = count;
-  octave_idx_type result_nc = 1;
+  for (ItType iter (nda, dir);
+      !(finished (i, 0, n_to_find, dir) || iter.iter_finished (dir)); iter.next (dir))
+    {
+      switch(nargout)
+        {
+        case 3:
+          vals(i) = *iter;
+          //Fall through!
+        case 2:
+          retval[1](i) = iter.get_col() + 1;
+          retval[0](i) = iter.get_row() + 1;
+          break;
+        case 1:
+          retval[0](i) = iter.get_idx() + 1;
+          break;
+        }
+      i += direction;
+    }
+  octave_value_list res (nargout);
+  for(int j = 0; j < retsize; ++j)
+    {
+      res(j) = retval[j];
+    }
+  if (nargout >= 3)
+    {
+      res(2) = vals;
+    }
+  return res;
+}
 
-  bool scalar_arg = false;
+template<int nargout, int direction, typename T>
+octave_value_list find_nonzero_elem_idx_for(const Array<T>& nda,
+    octave_idx_type n_to_find)
+{
+  return find_nonzero_elem_idx_iterized<nargout, direction, T, Array<T>,
+      array_nz_iter<T> > (nda, n_to_find);
+}
 
-  if (v.rows () == 1)
+template<int nargout, int direction, typename T>
+octave_value_list find_nonzero_elem_idx_for(const Sparse<T>& nda,
+    octave_idx_type n_to_find)
+{
+  return find_nonzero_elem_idx_iterized<nargout, direction, T, Sparse<T>,
+      sparse_nz_iter<T> > (nda, n_to_find);
+}
+
+template<int nargout, int direction>
+octave_value_list find_nonzero_elem_idx_for(const PermMatrix& nda,
+    octave_idx_type n_to_find)
+{
+  return find_nonzero_elem_idx_iterized<nargout, direction, double, PermMatrix,
+      perm_nz_iter> (nda.is_col_perm () ? nda : nda.flip_repr (), n_to_find);
+}
+
+
+template<int nargout, typename M>
+octave_value_list
+find_nonzero_elem_idx_for_nargout(const M& nda,
+    octave_idx_type n_to_find, int direction)
+{
+  if (direction > 0)
+    return find_nonzero_elem_idx_for<nargout, 1> (nda, n_to_find);
+  else if(direction < 0)
+    return find_nonzero_elem_idx_for<nargout, -1> (nda, n_to_find);
+  else
     {
-      result_nr = 1;
-      result_nc = count;
+      assert(false);
+      return octave_value_list ();
+    }
+}
 
-      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));
-    }
-
+template<typename M>
+octave_value_list find_nonzero_elem_idx(const M& nda, int nargout,
+    octave_idx_type n_to_find, int direction)
+{
   switch (nargout)
     {
     case 0:
     case 1:
-      retval(0) = idx;
-      break;
-
-    case 5:
-      retval(4) = nc;
-      // Fall through
-
-    case 4:
-      retval(3) = nr;
-      // Fall through
-
+      return find_nonzero_elem_idx_for_nargout<1> (nda, n_to_find, direction);
+    case 2:
+      return find_nonzero_elem_idx_for_nargout<2> (nda, n_to_find, direction);
     case 3:
-      retval(2) = val;
-      // Fall through!
-
-    case 2:
-      retval(1) = j_idx;
-      retval(0) = i_idx;
-      break;
-
+      return find_nonzero_elem_idx_for_nargout<3> (nda, n_to_find, direction);
     default:
-      panic_impossible ();
-      break;
+      assert(false);
+      return octave_value_list();
     }
-
-  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;
-      n_to_find = nc;
-      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 ();
-      if (v.is_col_perm ())
-        {
-          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
-        {
-          for (octave_idx_type k = 0; k < count; k++)
-            {
-              OCTAVE_QUIT;
-              const octave_idx_type i = start_nc + k;
-              const octave_idx_type j = p[i];
-              // Scatter into the index arrays according to
-              // j adjusted by the start point.
-              const octave_idx_type koff = j - start_nc;
-              i_idx(koff) = static_cast<double> (1+i);
-              j_idx(koff) = static_cast<double> (1+j);
-              idx(koff) = 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,
@@ -447,6 +470,15 @@
 
   octave_value arg = args(0);
 
+  if (nargout <= 0)
+    {
+      nargout = 1;
+    }
+  else if(nargout > 3)
+    {
+      nargout = 3;
+    }
+
   if (arg.is_bool_type ())
     {
       if (arg.is_sparse_type ())
diff -r 89a8a90ac850 liboctave/array/PermMatrix.cc
--- a/liboctave/array/PermMatrix.cc	Tue Mar 18 11:06:13 2014 +0100
+++ b/liboctave/array/PermMatrix.cc	Wed Mar 19 18:27:37 2014 -0600
@@ -216,3 +216,15 @@
 
   return r;
 }
+
+PermMatrix PermMatrix::flip_repr() const
+{
+  octave_idx_type n = dim2();
+  PermMatrix r(n);
+  r._colp = !_colp;
+  for (octave_idx_type i = 0; i < n; ++i)
+    {
+      r.xelem(xelem(i)) = i;
+    }
+  return r;
+}
diff -r 89a8a90ac850 liboctave/array/PermMatrix.h
--- a/liboctave/array/PermMatrix.h	Tue Mar 18 11:06:13 2014 +0100
+++ b/liboctave/array/PermMatrix.h	Wed Mar 19 18:27:37 2014 -0600
@@ -122,6 +122,8 @@
 
   static PermMatrix eye (octave_idx_type n);
 
+  PermMatrix flip_repr() const;
+
 private:
   bool _colp;
 };
diff -r 89a8a90ac850 test/build-sparse-tests.sh
--- a/test/build-sparse-tests.sh	Tue Mar 18 11:06:13 2014 +0100
+++ b/test/build-sparse-tests.sh	Wed Mar 19 18:27:37 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)));
