# HG changeset patch
# User dbateman
# Date 1201805951 0
# Node ID 402168152bb930dde56494826c026576c9125787
# Parent  3c999b2b5de83d8196a5adcd5c8f69b3f0263d33
[project @ 2008-01-31 18:59:09 by dbateman]

diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-C.cc
--- a/liboctave/Array-C.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-C.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -31,6 +31,65 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+static double
+xabs (const Complex& x)
+{
+  return (xisinf (x.real ()) || xisinf (x.imag ())) ? octave_Inf : abs (x);
+}
+
+static bool
+operator < (const Complex& a, const Complex& b)
+{
+  return (xisnan (b) || (xabs (a) < xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) < arg (b))));
+}
+
+static bool
+operator > (const Complex& a, const Complex& b)
+{
+  return (xisnan (a) || (xabs (a) > xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) > arg (b))));
+}
+
+template <>
+bool
+ascending_compare (Complex a, Complex b)
+{
+  return (xisnan (b) || (xabs (a) < xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) < arg (b))));
+}
+
+template <>
+bool
+ascending_compare (vec_index<Complex> *a, vec_index<Complex> *b)
+{
+  return (xisnan (b->vec)
+	  || (xabs (a->vec) < xabs (b->vec))
+	  || ((xabs (a->vec) == xabs (b->vec))
+	      && (arg (a->vec) < arg (b->vec))));
+}
+
+template <>
+bool
+descending_compare (Complex a, Complex b)
+{
+  return (xisnan (a) || (xabs (a) > xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) > arg (b))));
+}
+
+template <>
+bool
+descending_compare (vec_index<Complex> *a, vec_index<Complex> *b)
+{
+  return (xisnan (a->vec)
+	  || (xabs (a->vec) > xabs (b->vec))
+	  || ((xabs (a->vec) == xabs (b->vec))
+	      && (arg (a->vec) > arg (b->vec))));
+}
+
+INSTANTIATE_ARRAY_SORT (Complex);
 
 INSTANTIATE_ARRAY_AND_ASSIGN (Complex, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-b.cc
--- a/liboctave/Array-b.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-b.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -29,6 +29,9 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+INSTANTIATE_ARRAY_SORT (bool);
 
 INSTANTIATE_ARRAY_AND_ASSIGN (bool, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-ch.cc
--- a/liboctave/Array-ch.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-ch.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -29,6 +29,9 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+INSTANTIATE_ARRAY_SORT (char);
 
 INSTANTIATE_ARRAY_AND_ASSIGN (char, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-d.cc
--- a/liboctave/Array-d.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-d.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -29,6 +29,352 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+#if defined (HAVE_IEEE754_DATA_FORMAT)
+
+static inline uint64_t
+FloatFlip (uint64_t f)
+{
+  uint64_t mask
+    = -static_cast<int64_t>(f >> 63) | 0x8000000000000000ULL;
+
+  return f ^ mask;
+}
+
+static inline uint64_t
+IFloatFlip (uint64_t f)
+{
+  uint64_t mask = ((f >> 63) - 1) | 0x8000000000000000ULL;
+
+  return f ^ mask;
+}
+
+template <>
+bool
+ascending_compare (double a, double b)
+{
+  return (xisnan (b) || (a < b));
+}
+
+template <>
+bool
+ascending_compare (vec_index<double> *a, vec_index<double> *b)
+{
+  return (xisnan (b->vec) || (a->vec < b->vec));
+}
+
+template <>
+bool
+descending_compare (double a, double b)
+{
+  return (xisnan (a) || (a > b));
+}
+
+template <>
+bool
+descending_compare (vec_index<double> *a, vec_index<double> *b)
+{
+  return (xisnan (b->vec) || (a->vec > b->vec));
+}
+
+INSTANTIATE_ARRAY_SORT (uint64_t);
+
+template <>
+Array<double>
+Array<double>::sort (octave_idx_type dim, sortmode mode) const
+{
+  Array<double> m = *this;
+
+  dim_vector dv = m.dims ();
+
+  if (m.length () < 1)
+    return m;
+
+  octave_idx_type ns = dv(dim);
+  octave_idx_type iter = dv.numel () / ns;
+  octave_idx_type stride = 1;
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  double *v = m.fortran_vec ();
+
+  uint64_t *p = reinterpret_cast<uint64_t *> (v);
+
+  octave_sort<uint64_t> lsort;
+
+  if (mode == ASCENDING)
+    lsort.set_compare (ascending_compare);
+  else if (mode == DESCENDING)
+    lsort.set_compare (descending_compare);
+
+  if (stride == 1)
+    {
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	  // Flip the data in the vector so that int compares on
+	  // IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    p[i] = FloatFlip (p[i]);
+	      
+	  lsort.sort (p, ns);
+
+	  // Flip the data out of the vector so that int compares
+	  // on IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    p[i] = IFloatFlip (p[i]);
+
+	  // There are two representations of NaN.  One will be
+	  // sorted to the beginning of the vector and the other
+	  // to the end.  If it will be sorted incorrectly, fix
+	  // things up.
+
+	  if (lo_ieee_signbit (octave_NaN))
+	    {
+	      if (mode == UNDEFINED || mode == ASCENDING)
+		{
+		  octave_idx_type i = 0;
+		  double *vtmp = reinterpret_cast<double *> (p);
+		  while (xisnan (vtmp[i++]) && i < ns);
+		  for (octave_idx_type l = 0; l < ns - i + 1; l++)
+		    vtmp[l] = vtmp[l+i-1];
+		  for (octave_idx_type l = ns - i + 1; l < ns; l++)
+		    vtmp[l] = octave_NaN;
+		}
+	      else
+		{
+		  octave_idx_type i = ns;
+		  double *vtmp = reinterpret_cast<double *> (p);
+		  while (xisnan (vtmp[--i]) && i > 0);
+		  for (octave_idx_type l = i; l >= 0; l--)
+		    vtmp[l-i+ns-1] = vtmp[l];
+		  for (octave_idx_type l = 0; l < ns - i - 1; l++)
+		    vtmp[l] = octave_NaN;
+		}
+	    }
+
+	  p += ns;
+	}
+    }
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (uint64_t, vi, ns);
+
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	  octave_idx_type offset = j;
+	  octave_idx_type offset2 = 0;
+	  while (offset >= stride)
+	    {
+	      offset -= stride;
+	      offset2++;
+	    }
+	  offset += offset2 * stride * ns;
+
+	  // Flip the data in the vector so that int compares on
+	  // IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    vi[i] = FloatFlip (p[i*stride + offset]);
+
+	  lsort.sort (vi, ns);
+
+	  // Flip the data out of the vector so that int compares
+	  // on IEEE754 give the correct ordering.
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    p[i*stride + offset] = IFloatFlip (vi[i]);
+	      
+	  // There are two representations of NaN. One will be
+	  // sorted to the beginning of the vector and the other
+	  // to the end. If it will be sorted to the beginning,
+	  // fix things up.
+
+	  if (lo_ieee_signbit (octave_NaN))
+	    {
+	      if (mode == UNDEFINED || mode == ASCENDING)
+		{
+		   octave_idx_type i = 0;
+		  while (xisnan (v[i++*stride + offset]) && i < ns);
+		  for (octave_idx_type l = 0; l < ns - i + 1; l++)
+		    v[l*stride + offset] = v[(l+i-1)*stride + offset];
+		  for (octave_idx_type l = ns - i + 1; l < ns; l++)
+		    v[l*stride + offset] = octave_NaN;
+		}
+	      else
+		{
+		   octave_idx_type i = ns;
+		  while (xisnan (v[--i*stride + offset]) && i > 0);
+		  for (octave_idx_type l = i; l >= 0; l--)
+		    v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
+		  for (octave_idx_type l = 0; l < ns - i - 1; l++)
+		    v[l*stride + offset] = octave_NaN;
+		}
+	    }
+	}
+    }
+
+  return m;
+}
+
+template <>
+Array<double>
+Array<double>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim, 
+		     sortmode mode) const
+{
+  Array<double> m = *this;
+
+  dim_vector dv = m.dims ();
+
+  if (m.length () < 1)
+    {
+      sidx = Array<octave_idx_type> (dv);
+      return m;
+    }
+
+  octave_idx_type ns = dv(dim);
+  octave_idx_type iter = dv.numel () / ns;
+  octave_idx_type stride = 1;
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  double *v = m.fortran_vec ();
+
+  uint64_t *p = reinterpret_cast<uint64_t *> (v);
+
+  octave_sort<vec_index<uint64_t> *> indexed_sort;
+
+  if (mode == ASCENDING)
+    indexed_sort.set_compare (ascending_compare);
+  else if (mode == DESCENDING)
+    indexed_sort.set_compare (descending_compare);
+
+  OCTAVE_LOCAL_BUFFER (vec_index<uint64_t> *, vi, ns);
+  OCTAVE_LOCAL_BUFFER (vec_index<uint64_t>, vix, ns);
+  
+  for (octave_idx_type i = 0; i < ns; i++)
+    vi[i] = &vix[i];
+
+  sidx = Array<octave_idx_type> (dv);
+      
+  for (octave_idx_type j = 0; j < iter; j++)
+    {
+      octave_idx_type offset = j;
+      octave_idx_type offset2 = 0;
+      while (offset >= stride)
+	{
+	  offset -= stride;
+	  offset2++;
+	}
+      offset += offset2 * stride * ns;
+
+      // Flip the data in the vector so that int compares on
+      // IEEE754 give the correct ordering.
+
+      for (octave_idx_type i = 0; i < ns; i++)
+	{
+	  vi[i]->vec = FloatFlip (p[i*stride + offset]);
+	  vi[i]->indx = i;
+	}
+
+      indexed_sort.sort (vi, ns);
+
+      // Flip the data out of the vector so that int compares on
+      // IEEE754 give the correct ordering
+
+      for (octave_idx_type i = 0; i < ns; i++)
+	{
+	  p[i*stride + offset] = IFloatFlip (vi[i]->vec);
+	  sidx(i*stride + offset) = vi[i]->indx;
+	}
+
+      // There are two representations of NaN.  One will be sorted
+      // to the beginning of the vector and the other to the end.
+      // If it will be sorted to the beginning, fix things up.
+
+      if (lo_ieee_signbit (octave_NaN))
+	{
+	  if (mode == UNDEFINED || mode == ASCENDING)
+	    {
+	      octave_idx_type i = 0;
+	      while (xisnan (v[i++*stride+offset]) && i < ns);
+	      OCTAVE_LOCAL_BUFFER (double, itmp, i - 1);
+	      for (octave_idx_type l = 0; l < i -1; l++)
+		itmp[l] = sidx(l*stride + offset);
+	      for (octave_idx_type l = 0; l < ns - i + 1; l++)
+		{
+		  v[l*stride + offset] = v[(l+i-1)*stride + offset];
+		  sidx(l*stride + offset) = sidx((l+i-1)*stride + offset);
+		}
+	      for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++)
+		{
+		  v[l*stride + offset] = octave_NaN;
+		  sidx(l*stride + offset) = 
+		    static_cast<octave_idx_type>(itmp[k]);
+		}
+	    }
+	  else 
+	    {
+	      octave_idx_type i = ns;
+	      while (xisnan (v[--i*stride+offset]) && i > 0);
+	      OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1);
+	      for (octave_idx_type l = 0; l < ns - i -1; l++)
+		itmp[l] = sidx((l+i+1)*stride + offset);
+	      for (octave_idx_type l = i; l >= 0; l--)
+		{
+		  v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
+		  sidx((l-i+ns-1)*stride + offset) = sidx(l*stride + offset);
+		}
+	      for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++)
+		{
+		  v[l*stride + offset] = octave_NaN;
+		  sidx(l*stride + offset) = 
+		    static_cast<octave_idx_type>(itmp[k]);
+		}
+	    }
+	}
+    }
+
+  return m;
+}
+
+#else
+
+template <>
+bool
+Array<double>::ascending_compare (double a, double b) const
+{
+  return (xisnan (b) || (a < b));
+}
+
+template <>
+bool
+Array<double>::ascending_compare (vec_index<double> *a, 
+				  vec_index<double> *b) const
+{
+  return (xisnan (b->vec) || (a->vec < b->vec));
+}
+
+template <>
+bool
+Array<double>::descending_compare (double a, double b) const
+{
+  return (xisnan (a) || (a > b));
+}
+
+template <>
+bool
+Array<double>::descending_compare (vec_index<double> *a, 
+				   vec_index<double> *b) const
+{
+  return (xisnan (b->vec) || (a->vec > b->vec));
+}
+
+INSTANTIATE_ARRAY_SORT (double);
+
+#endif
 
 INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-i.cc
--- a/liboctave/Array-i.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-i.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -31,6 +31,10 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+INSTANTIATE_ARRAY_SORT (int);
+INSTANTIATE_ARRAY_SORT (long);
 
 INSTANTIATE_ARRAY_AND_ASSIGN (int, OCTAVE_API);
 INSTANTIATE_ARRAY_AND_ASSIGN (long, OCTAVE_API);
@@ -38,11 +42,21 @@
 INSTANTIATE_ARRAY_ASSIGN (int, short, OCTAVE_API);
 INSTANTIATE_ARRAY_ASSIGN (int, char, OCTAVE_API);
 
+INSTANTIATE_ARRAY_SORT (octave_int8);
+INSTANTIATE_ARRAY_SORT (octave_int16);
+INSTANTIATE_ARRAY_SORT (octave_int32);
+INSTANTIATE_ARRAY_SORT (octave_int64);
+
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_int8, OCTAVE_API);
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_int16, OCTAVE_API);
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_int32, OCTAVE_API);
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_int64, OCTAVE_API);
 
+INSTANTIATE_ARRAY_SORT (octave_uint8);
+INSTANTIATE_ARRAY_SORT (octave_uint16);
+INSTANTIATE_ARRAY_SORT (octave_uint32);
+INSTANTIATE_ARRAY_SORT (octave_uint64);
+
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_uint8, OCTAVE_API);
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_uint16, OCTAVE_API);
 INSTANTIATE_ARRAY_AND_ASSIGN (octave_uint32, OCTAVE_API);
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-idx-vec.cc
--- a/liboctave/Array-idx-vec.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-idx-vec.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -31,6 +31,8 @@
 #include "Array.h"
 #include "Array.cc"
 
+NO_INSTANTIATE_ARRAY_SORT (idx_vector);
+
 INSTANTIATE_ARRAY (idx_vector, OCTAVE_API);
 
 /*
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-s.cc
--- a/liboctave/Array-s.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-s.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -29,6 +29,9 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
+
+INSTANTIATE_ARRAY_SORT (short);
 
 INSTANTIATE_ARRAY_AND_ASSIGN (short, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array-str.cc
--- a/liboctave/Array-str.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array-str.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -28,9 +28,12 @@
 
 #include "Array.h"
 #include "Array.cc"
+#include "oct-sort.cc"
 
 #include <string>
 
+INSTANTIATE_ARRAY_SORT (std::string);
+
 INSTANTIATE_ARRAY (std::string, OCTAVE_API);
 
 /*
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array.cc
--- a/liboctave/Array.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -2447,6 +2447,189 @@
   return retval;
 }
 
+template <class T>
+bool 
+ascending_compare (T a, T b)
+{
+  return (a < b);
+}
+
+template <class T>
+bool 
+descending_compare (T a, T b)
+{
+  return (a > b);
+}
+
+template <class T>
+bool 
+ascending_compare (vec_index<T> *a, vec_index<T> *b)
+{
+  return (a->vec < b->vec);
+}
+
+template <class T>
+bool 
+descending_compare (vec_index<T> *a, vec_index<T> *b)
+{
+  return (a->vec > b->vec);
+}
+
+template <class T>
+Array<T>
+Array<T>::sort (octave_idx_type dim, sortmode mode) const
+{
+  Array<T> m = *this;
+
+  dim_vector dv = m.dims ();
+
+  if (m.length () < 1)
+    return m;
+
+  octave_idx_type ns = dv(dim);
+  octave_idx_type iter = dv.numel () / ns;
+  octave_idx_type stride = 1;
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  T *v = m.fortran_vec ();
+  octave_sort<T> lsort;
+  
+  if (mode == ASCENDING) 
+    lsort.set_compare (ascending_compare);
+  else if (mode == DESCENDING)
+    lsort.set_compare (descending_compare);
+
+  if (stride == 1)
+    {
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	  lsort.sort (v, ns);
+	  v += ns;
+	}
+    }
+  else
+    {
+      // Don't use OCTAVE_LOCAL_BUFFER here as it doesn't work with bool
+      // on some compilers.
+      Array<T> vi (ns);
+      T *pvi = vi.fortran_vec ();
+
+      for (octave_idx_type j = 0; j < iter; j++) 
+	{
+	   octave_idx_type offset = j;
+	   octave_idx_type offset2 = 0;
+	  while (offset >= stride)
+	    {
+	      offset -= stride;
+	      offset2++;
+	    }
+	  offset += offset2 * stride * ns;
+	  
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    pvi[i] = v[i*stride + offset];
+
+	  lsort.sort (pvi, ns);
+	      
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    v[i*stride + offset] = pvi[i];
+	}
+    }
+
+  return m;
+}
+
+template <class T>
+Array<T>
+Array<T>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim, 
+		sortmode mode) const
+{
+  Array<T> m = *this;
+
+  dim_vector dv = m.dims ();
+
+  if (m.length () < 1)
+    {
+      sidx = Array<octave_idx_type> (dv);
+      return m;
+    }
+
+  octave_idx_type ns = dv(dim);
+  octave_idx_type iter = dv.numel () / ns;
+  octave_idx_type stride = 1;
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  T *v = m.fortran_vec ();
+  octave_sort<vec_index<T> *> indexed_sort;
+
+  if (mode == ASCENDING) 
+    indexed_sort.set_compare (ascending_compare);
+  else if (mode == DESCENDING)
+    indexed_sort.set_compare (descending_compare);
+
+  OCTAVE_LOCAL_BUFFER (vec_index<T> *, vi, ns);
+  OCTAVE_LOCAL_BUFFER (vec_index<T>, vix, ns);
+
+  for (octave_idx_type i = 0; i < ns; i++)
+    vi[i] = &vix[i];
+
+  sidx = Array<octave_idx_type> (dv);
+      
+  if (stride == 1)
+    {
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	   octave_idx_type offset = j * ns;
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    {
+	      vi[i]->vec = v[i];
+	      vi[i]->indx = i;
+	    }
+
+	  indexed_sort.sort (vi, ns);
+
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    {
+	      v[i] = vi[i]->vec;
+	      sidx(i + offset) = vi[i]->indx;
+	    }
+	  v += ns;
+	}
+    }
+  else
+    {
+      for (octave_idx_type j = 0; j < iter; j++)
+	{
+	  octave_idx_type offset = j;
+	  octave_idx_type offset2 = 0;
+	  while (offset >= stride)
+	    {
+	      offset -= stride;
+	      offset2++;
+	    }
+	  offset += offset2 * stride * ns;
+	      
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    {
+	      vi[i]->vec = v[i*stride + offset];
+	      vi[i]->indx = i;
+	    }
+
+	  indexed_sort.sort (vi, ns);
+	      
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    {
+	      v[i*stride+offset] = vi[i]->vec;
+	      sidx(i*stride+offset) = vi[i]->indx;
+	    }
+	}
+    }
+
+  return m;
+}
+
 // FIXME -- this is a mess.
 
 template <class LT, class RT>
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array.h
--- a/liboctave/Array.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array.h	Thu Jan 31 18:59:11 2008 +0000
@@ -33,6 +33,8 @@
 #include "dim-vector.h"
 #include "lo-utils.h"
 
+#include "oct-sort.h"
+
 class idx_vector;
 
 // One dimensional array class.  Handles the reference counting for
@@ -543,6 +545,10 @@
   // Unsafe.  This function exists to support the MEX interface.
   // You should not use it anywhere else.
   void *mex_get_data (void) const { return const_cast<T *> (data ()); }
+
+  Array<T> sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const;
+  Array<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		 sortmode mode = UNDEFINED) const;
 };
 
 // NOTE: these functions should be friends of the Array<T> class and
@@ -588,6 +594,25 @@
   INSTANTIATE_ARRAY (T, API); \
   INSTANTIATE_ARRAY_ASSIGN (T, T, API)
 
+#define INSTANTIATE_ARRAY_SORT(T) \
+  template class octave_sort<T>; \
+  template class vec_index<T>; \
+  template class octave_sort<vec_index<T> *>;
+
+#define NO_INSTANTIATE_ARRAY_SORT(T) \
+  template class vec_index<T>; \
+  template <> bool ascending_compare (T, T) { return true; } \
+  template <> bool ascending_compare (vec_index<T> *, vec_index<T> *) \
+    { return true; } \
+  template <> bool descending_compare (T, T) { return true; } \
+  template <> bool descending_compare (vec_index<T> *, vec_index<T> *) \
+    { return true; } \
+  template <> Array<T> Array<T>::sort \
+    (octave_idx_type, sortmode) const { return *this; } \
+  template <> Array<T> Array<T>::sort (Array<octave_idx_type> &sidx, \
+    octave_idx_type, sortmode) const \
+    { sidx = Array<octave_idx_type> (); return *this; }
+
 #endif
 
 /*
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array2.h
--- a/liboctave/Array2.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array2.h	Thu Jan 31 18:59:11 2008 +0000
@@ -115,6 +115,19 @@
       Array<T> tmp = Array<T>::index (i, j, resize_ok, rfv);
       return Array2<T> (tmp, tmp.rows (), tmp.columns ());
     }
+
+  Array2<T> sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+    {
+      Array<T> tmp = Array<T>::sort (dim, mode);
+      return Array2<T> (tmp, tmp.rows (), tmp.columns ());
+    }
+
+  Array2<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		 sortmode mode = UNDEFINED) const
+    {
+      Array<T> tmp = Array<T>::sort (sidx, dim, mode);
+      return Array2<T> (tmp, tmp.rows (), tmp.columns ());
+    }
 };
 
 #endif
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Array3.h
--- a/liboctave/Array3.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Array3.h	Thu Jan 31 18:59:11 2008 +0000
@@ -75,6 +75,19 @@
 
   void resize (octave_idx_type r, octave_idx_type c, octave_idx_type p, const T& val)
     { this->resize_and_fill (r, c, p, val); }
+
+  Array3<T> sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+    {
+      Array<T> tmp = Array<T>::sort (dim, mode);
+      return Array3<T> (tmp, tmp.rows (), tmp.columns (), tmp.pages ());
+    }
+
+  Array3<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		 sortmode mode = UNDEFINED) const
+    {
+      Array<T> tmp = Array<T>::sort (sidx, dim, mode);
+      return Array3<T> (tmp, tmp.rows (), tmp.columns (), tmp.pages ());
+    }
 };
 
 #endif
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/ArrayN.h
--- a/liboctave/ArrayN.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/ArrayN.h	Thu Jan 31 18:59:11 2008 +0000
@@ -135,6 +135,19 @@
       Array<T> tmp = Array<T>::index (ra_idx, resize_ok, rfv);
       return ArrayN<T> (tmp, tmp.dims ());
     }
+
+  ArrayN<T> sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+    {
+      Array<T> tmp = Array<T>::sort (dim, mode);
+      return ArrayN<T> (tmp, tmp.dims ());
+    }
+
+  ArrayN<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		 sortmode mode = UNDEFINED) const
+    {
+      Array<T> tmp = Array<T>::sort (sidx, dim, mode);
+      return ArrayN<T> (tmp, tmp.dims ());
+    }
 };
 
 template <class T>
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/CSparse.cc
--- a/liboctave/CSparse.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/CSparse.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -46,8 +46,6 @@
 #include "SparseCmplxCHOL.h"
 #include "SparseCmplxQR.h"
 
-#include "oct-sort.h"
-
 // Define whether to use a basic QR solver or one that uses a Dulmange
 // Mendelsohn factorization to seperate the problem into under-determined,
 // well-determined and over-determined parts and solves them seperately
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/ChangeLog	Thu Jan 31 18:59:11 2008 +0000
@@ -1,3 +1,71 @@
+2008-01-31  David Bateman  <dbateman@free.fr>
+
+	* oct-sort.cc: conversion of int to octave_idx_type where needed
+	for 64-bit builds.
+	(IFLT): Allow IFLT macro to be overridden.
+	* oct-sort.h: conversion of int to octave_idx_type where needed
+	for 64-bit builds.
+	(enum sortmode): Type of sort to perform.
+	(vec_index): Simple class to aid in indexed sorts.
+	
+	* Array.h ( Array<T> sort (octave_idx_type, sortmode) const,
+	Array<T> sort (Array<octave_idx_type> &, octave_idx_type,
+	sortmode) const): Array sorting methods.
+	(INSTANTIATE_ARRAY_SORT, NO_INSTANTIATE_ARRAY_SORT): Macros to
+	instantiate the array sorting methods.
+	* Array.cc (ascending_compare, descending_compare): New template
+	functions for generic sort comparison.
+  	( Array<T> Array<T>::sort (octave_idx_type, sortmode) const,
+	Array<T> Array<T>::sort (Array<octave_idx_type> &, octave_idx_type,
+	sortmode) const): Array sorting functions based of octave_sort
+	class.
+	* Array-C.cc: Instantiate the complex array sort methods. 
+	(IFLT): New macro to override the one in the
+	octave_sort class to avoid need for Complex < and > operators.
+	(static double xabs (const Complex&)): Complex abs function
+	avoiding std::abs(Inf) returning NaN with some compilers.
+	(ascending_compare, descending compare): override template
+	functions for complex comparison.
+	* Array-d.cc: Instantiate the double array sort methods. 
+	(Array<double> Array<double>::sort (octave_idx_type, 
+	sortmode) const, Array<double> Array<double>::sort 
+	(Array<octave_idx_type> &, octave_idx_type, sortmode) const): 
+	Array sorting functions based of octave_sort using uint64 sorting
+	on IEE754 doubles, for speed and correct sorting of Inf and NaN.
+	(ascending_compare, descending compare): override template
+	functions for double and uint64 comparison.
+	* Array-b.cc, Array-ch.cc, Array-i.cc, Array-s.cc, Array-str.cc: 
+	Instantiate the array sort methods.
+	* Array-idx-vec.cc: Null instantiation of array sort methods.
+	* Array2.h, Array3.h, ArrayN.h (sort): 2, 3 and N-dimensional
+	versions of the sort methods based on Array<T>::sort.
+
+	* CSparse.cc, dSparse.cc: Remove inclusion of octa-sort.h.
+	* Sparse.h ( Sparse<T> sort (octave_idx_type, sortmode) const,
+	Sparse<T> sort (Array<octave_idx_type> &, octave_idx_type,
+	sortmode) const): Sparse sorting methods.
+	(INSTANTIATE_ARRAY_SORT): Macro to instantiate the sparse sorting 
+	methods.
+	* Sparse.cc: replace sort with lsort throughout to avoid shadowing
+	of new sort method.
+	(sparse_ascending_compare, sparse_descending_compare): New template
+	functions for generic sort comparison.
+  	( Sparse<T> Sparse<T>::sort (octave_idx_type, sortmode) const,
+	Sparse<T> Sparse<T>::sort (Sparse<octave_idx_type> &, octave_idx_type,
+	sortmode) const): Sparse sorting functions based of octave_sort
+	class.
+	* Sparse-C.cc: Instantiate the complex sparse sort methods. 
+	(IFLT): New macro to override the one in the
+	octave_sort class to avoid need for Complex < and > operators.
+	(static double xabs (const Complex&)): Complex abs function
+	avoiding std::abs(Inf) returning NaN with some compilers.
+	(sparse_ascending_compare, sparse_descending compare): override
+	template functions for complex comparison.
+	* Sparse-d.cc: Instantiate the cdouble sparse sort methods. 
+	(sparse_ascending_compare, sparse_descending compare): override
+	template functions for double comparison.
+	* Array-b.cc: Instantiate the sparse sort methods.
+
 2008-01-25  Jaroslav Hajek  <highegg@gmail.com>
 
 	* idx-vector.h (idx_vector::idx_vector_rep::range_base,
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Sparse-C.cc
--- a/liboctave/Sparse-C.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Sparse-C.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -28,10 +28,71 @@
 // Instantiate Sparse matrix of complex values.
 
 #include "oct-cmplx.h"
-
+#include "lo-mappers.h"
+#include "lo-ieee.h"
 #include "Sparse.h"
 #include "Sparse.cc"
 
+#include "oct-sort.cc"
+
+static double
+xabs (const Complex& x)
+{
+  return (xisinf (x.real ()) || xisinf (x.imag ())) ? octave_Inf : abs (x);
+}
+
+static bool
+operator < (const Complex& a, const Complex& b)
+{
+  return (xisnan (b) || (xabs (a) < xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) < arg (b))));
+}
+
+static bool
+operator > (const Complex& a, const Complex& b)
+{
+  return (xisnan (a) || (xabs (a) > xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) > arg (b))));
+}
+
+template <>
+bool
+sparse_ascending_compare (Complex a, Complex b)
+{
+  return (xisnan (b) || (xabs (a) < xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) < arg (b))));
+}
+
+template <>
+bool
+sparse_ascending_compare (vec_index<Complex> *a, vec_index<Complex> *b)
+{
+  return (xisnan (b->vec)
+	  || (xabs (a->vec) < xabs (b->vec))
+	  || ((xabs (a->vec) == xabs (b->vec))
+	      && (arg (a->vec) < arg (b->vec))));
+}
+
+template <>
+bool
+sparse_descending_compare (Complex a, Complex b)
+{
+  return (xisnan (a) || (xabs (a) > xabs (b))
+	  || ((xabs (a) == xabs (b)) && (arg (a) > arg (b))));
+}
+
+template <>
+bool
+sparse_descending_compare (vec_index<Complex> *a, vec_index<Complex> *b)
+{
+  return (xisnan (a->vec)
+	  || (xabs (a->vec) > xabs (b->vec))
+	  || ((xabs (a->vec) == xabs (b->vec))
+	      && (arg (a->vec) > arg (b->vec))));
+}
+
+INSTANTIATE_SPARSE_SORT (Complex);
+
 INSTANTIATE_SPARSE_AND_ASSIGN (Complex, OCTAVE_API);
 
 INSTANTIATE_SPARSE_ASSIGN (Complex, double, OCTAVE_API);
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Sparse-b.cc
--- a/liboctave/Sparse-b.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Sparse-b.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -29,6 +29,9 @@
 
 #include "Sparse.h"
 #include "Sparse.cc"
+#include "oct-sort.cc"
+
+INSTANTIATE_SPARSE_SORT (bool);
 
 INSTANTIATE_SPARSE_AND_ASSIGN (bool, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Sparse-d.cc
--- a/liboctave/Sparse-d.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Sparse-d.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -27,8 +27,40 @@
 
 // Instantiate Sparse matrix of double values.
 
+#include "lo-mappers.h"
 #include "Sparse.h"
 #include "Sparse.cc"
+#include "oct-sort.cc"
+
+template <>
+bool
+sparse_ascending_compare (double a, double b)
+{
+  return (xisnan (b) || (a < b));
+}
+
+template <>
+bool
+sparse_ascending_compare (vec_index<double> *a, vec_index<double> *b)
+{
+  return (xisnan (b->vec) || (a->vec < b->vec));
+}
+
+template <>
+bool
+sparse_descending_compare (double a, double b)
+{
+  return (xisnan (a) || (a > b));
+}
+
+template <>
+bool
+sparse_descending_compare (vec_index<double> *a, vec_index<double> *b)
+{
+  return (xisnan (b->vec) || (a->vec > b->vec));
+}
+
+INSTANTIATE_SPARSE_SORT (double);
 
 INSTANTIATE_SPARSE_AND_ASSIGN (double, OCTAVE_API);
 
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Sparse.cc
--- a/liboctave/Sparse.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Sparse.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -369,9 +369,9 @@
 	{
 	  OCTAVE_QUIT;
 	  octave_sort<octave_sparse_sort_idxl *> 
-	    sort (octave_sparse_sidxl_comp);
-
-	  sort.sort (sidx, actual_nzmx);
+	    lsort (octave_sparse_sidxl_comp);
+
+	  lsort.sort (sidx, actual_nzmx);
 	  OCTAVE_QUIT;
 
 	  // Now count the unique non-zero values
@@ -487,9 +487,9 @@
 	{
 	  OCTAVE_QUIT;
 	  octave_sort<octave_sparse_sort_idxl *> 
-	    sort (octave_sparse_sidxl_comp);
-
-	  sort.sort (sidx, actual_nzmx);
+	    lsort (octave_sparse_sidxl_comp);
+
+	  lsort.sort (sidx, actual_nzmx);
 	  OCTAVE_QUIT;
 
 	  // Now count the unique non-zero values
@@ -1855,7 +1855,7 @@
 
 	      OCTAVE_LOCAL_BUFFER (octave_idx_type, itmp, 
 				   (nr > nc ? nr : nc));
-	      octave_sort<octave_idx_type> sort;
+	      octave_sort<octave_idx_type> lsort;
 
 	      if (n > nr || m > nc)
 		permutation = false;
@@ -1871,7 +1871,7 @@
 		  // for idx_vector type.
 		  for (octave_idx_type i = 0; i < n; i++)
 		    itmp [i] = idx_i.elem (i);
-		  sort.sort (itmp, n);
+		  lsort.sort (itmp, n);
 		  for (octave_idx_type i = 1; i < n; i++)
 		    if (itmp[i-1] == itmp[i])
 		      {
@@ -1883,7 +1883,7 @@
 		{
 		  for (octave_idx_type i = 0; i < m; i++)
 		    itmp [i] = idx_j.elem (i);
-		  sort.sort (itmp, m);
+		  lsort.sort (itmp, m);
 		  for (octave_idx_type i = 1; i < m; i++)
 		    if (itmp[i-1] == itmp[i])
 		      {
@@ -1920,7 +1920,7 @@
 			      retval.xridx (kk++) = ii;
 			    }
 			}
-		      sort.sort (ri + retval.xcidx (j), kk - retval.xcidx (j));
+		      lsort.sort (ri + retval.xcidx (j), kk - retval.xcidx (j));
 		      for (octave_idx_type p = retval.xcidx (j); p < kk; p++)
 			retval.xdata (p) = X [retval.xridx (p)]; 
 		      retval.xcidx(j+1) = kk;
@@ -2022,7 +2022,7 @@
 				    }
 				}
 			    }
-			  sort.sort (ri + retval.xcidx (j), 
+			  lsort.sort (ri + retval.xcidx (j), 
 				     kk - retval.xcidx (j));
 			  for (octave_idx_type p = retval.xcidx (j); 
 			       p < kk; p++)
@@ -2053,6 +2053,215 @@
   return index (ra_idx (0), ra_idx (1), resize_ok);
 }
 
+// Can't use versions of these in Array.cc due to duplication of the 
+// instantiations for Array<double and Sparse<double>, etc
+template <class T>
+bool 
+sparse_ascending_compare (T a, T b)
+{
+  return (a < b);
+}
+
+template <class T>
+bool 
+sparse_descending_compare (T a, T b)
+{
+  return (a > b);
+}
+
+template <class T>
+bool 
+sparse_ascending_compare (vec_index<T> *a, vec_index<T> *b)
+{
+  return (a->vec < b->vec);
+}
+
+template <class T>
+bool 
+sparse_descending_compare (vec_index<T> *a, vec_index<T> *b)
+{
+  return (a->vec > b->vec);
+}
+
+template <class T>
+Sparse<T>
+Sparse<T>::sort (octave_idx_type dim, sortmode mode) const
+{
+  Sparse<T> m = *this;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  if (m.length () < 1)
+    return m;
+
+  if (dim > 0)
+    {
+      m = m.transpose ();
+      nr = m.rows ();
+      nc = m.columns ();
+    }
+
+  octave_sort<T> lsort;
+
+  if (mode == ASCENDING) 
+    lsort.set_compare (sparse_ascending_compare);
+  else if (mode == DESCENDING)
+    lsort.set_compare (sparse_descending_compare);
+
+  T *v = m.data ();
+  octave_idx_type *mcidx = m.cidx ();
+  octave_idx_type *mridx = m.ridx ();
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    {
+      octave_idx_type ns = mcidx [j + 1] - mcidx [j];
+      lsort.sort (v, ns);
+
+      octave_idx_type i;
+      if (mode == ASCENDING) 
+	{
+	  for (i = 0; i < ns; i++)
+	    if (sparse_ascending_compare (static_cast<T> (0), v [i]))
+	      break;
+	}
+      else
+	{
+	  for (i = 0; i < ns; i++)
+	    if (sparse_descending_compare (static_cast<T> (0), v [i]))
+	      break;
+	}
+      for (octave_idx_type k = 0; k < i; k++)
+	mridx [k] = k;
+      for (octave_idx_type k = i; k < ns; k++)
+	mridx [k] = k - ns + nr; 
+
+      v += ns;
+      mridx += ns;
+    }
+
+  if (dim > 0)
+      m = m.transpose ();
+
+  return m;
+}
+
+template <class T>
+Sparse<T>
+Sparse<T>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim, 
+		 sortmode mode) const
+{
+  Sparse<T> m = *this;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  if (m.length () < 1)
+    {
+      sidx = Array<octave_idx_type> (dim_vector (nr, nc));
+      return m;
+    }
+
+  if (dim > 0)
+    {
+      m = m.transpose ();
+      nr = m.rows ();
+      nc = m.columns ();
+    }
+
+  octave_sort<vec_index<T> *> indexed_sort;
+
+  if (mode == ASCENDING) 
+    indexed_sort.set_compare (sparse_ascending_compare);
+  else if (mode == DESCENDING)
+    indexed_sort.set_compare (sparse_descending_compare);
+
+  T *v = m.data ();
+  octave_idx_type *mcidx = m.cidx ();
+  octave_idx_type *mridx = m.ridx ();
+
+  OCTAVE_LOCAL_BUFFER (vec_index<T> *, vi, nr);
+  OCTAVE_LOCAL_BUFFER (vec_index<T>, vix, nr);
+
+  for (octave_idx_type i = 0; i < nr; i++)
+    vi[i] = &vix[i];
+
+  sidx = Array<octave_idx_type> (dim_vector (nr, nc));
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    {
+      octave_idx_type ns = mcidx [j + 1] - mcidx [j];
+      octave_idx_type offset = j * nr;
+
+      if (ns == 0)
+	{
+	  for (octave_idx_type k = 0; k < nr; k++)
+	    sidx (offset + k) = k;
+	}
+      else
+	{
+	  for (octave_idx_type i = 0; i < ns; i++)
+	    {
+	      vi[i]->vec = v[i];
+	      vi[i]->indx = mridx[i];
+	    }
+
+	  indexed_sort.sort (vi, ns);
+
+	  octave_idx_type i;
+	  if (mode == ASCENDING) 
+	    {
+	      for (i = 0; i < ns; i++)
+		if (sparse_ascending_compare (static_cast<T> (0), 
+					      vi [i] -> vec))
+		  break;
+	    }
+	  else
+	    {
+	      for (i = 0; i < ns; i++)
+		if (sparse_descending_compare (static_cast<T> (0), 
+					       vi [i] -> vec))
+		  break;
+	    }
+
+	  octave_idx_type ii = 0;
+	  octave_idx_type jj = i;
+	  for (octave_idx_type k = 0; k < nr; k++)
+	    {
+	      if (ii < ns && mridx[ii] == k)
+		ii++;
+	      else
+		sidx (offset + jj++) = k;
+	    }
+
+	  for (octave_idx_type k = 0; k < i; k++)
+	    {
+	      v [k] = vi [k] -> vec;
+	      sidx (k + offset) = vi [k] -> indx;
+	      mridx [k] = k;
+	    }
+
+	  for (octave_idx_type k = i; k < ns; k++)
+	    {
+	      v [k] = vi [k] -> vec;
+	      sidx (k - ns + nr + offset) = vi [k] -> indx;
+	      mridx [k] = k - ns + nr; 
+	    }
+
+	  v += ns;
+	  mridx += ns;
+	}
+    }
+
+  if (dim > 0)
+    {
+      m = m.transpose ();
+      sidx = sidx.transpose ();
+    }
+
+  return m;
+}
+
 // FIXME
 // Unfortunately numel can overflow for very large but very sparse matrices.
 // For now just flag an error when this happens.
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/Sparse.h
--- a/liboctave/Sparse.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/Sparse.h	Thu Jan 31 18:59:11 2008 +0000
@@ -35,6 +35,8 @@
 #include "dim-vector.h"
 #include "lo-utils.h"
 
+#include "oct-sort.h"
+
 class idx_vector;
 
 // Two dimensional sparse class.  Handles the reference counting for
@@ -515,6 +517,10 @@
   octave_idx_type *mex_get_ir (void) const { return const_cast<octave_idx_type *> (ridx ()); }
 
   octave_idx_type *mex_get_jc (void) const { return const_cast<octave_idx_type *> (cidx ()); }
+
+  Sparse<T> sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const;
+  Sparse<T> sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		 sortmode mode = UNDEFINED) const;
 };
 
 // NOTE: these functions should be friends of the Sparse<T> class and
@@ -540,6 +546,11 @@
   INSTANTIATE_SPARSE (T, API); \
   INSTANTIATE_SPARSE_ASSIGN (T, T, API)
 
+#define INSTANTIATE_SPARSE_SORT(T) \
+  template class octave_sort<T>; \
+  template class vec_index<T>; \
+  template class octave_sort<vec_index<T> *>;
+
 #endif
 
 /*
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/dSparse.cc
--- a/liboctave/dSparse.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/dSparse.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -47,8 +47,6 @@
 #include "SparsedbleCHOL.h"
 #include "SparseQR.h"
 
-#include "oct-sort.h"
-
 // Define whether to use a basic QR solver or one that uses a Dulmange
 // Mendelsohn factorization to seperate the problem into under-determined,
 // well-determined and over-determined parts and solves them seperately
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/oct-sort.cc
--- a/liboctave/oct-sort.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/oct-sort.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -90,7 +90,9 @@
 #include "quit.h"
 #include "oct-sort.h"
 
+#ifndef IFLT
 #define IFLT(a,b)  if (compare ? compare ((a), (b)) : ((a) < (b)))
+#endif
 
 template <class T>
 octave_sort<T>::octave_sort (void) : compare (0)
@@ -186,10 +188,10 @@
 Returns -1 in case of error.
 */
 template <class T>
-int
+octave_idx_type
 octave_sort<T>::count_run (T *lo, T *hi, int *descending)
 {
-  int n;
+  octave_idx_type n;
 
   *descending = 0;
   ++lo;
@@ -243,12 +245,12 @@
 Returns -1 on error.  See listsort.txt for info on the method.
 */
 template <class T>
-int
-octave_sort<T>::gallop_left (T key, T *a, int n, int hint)
+octave_idx_type
+octave_sort<T>::gallop_left (T key, T *a, octave_idx_type n, octave_idx_type hint)
 {
-  int ofs;
-  int lastofs;
-  int k;
+  octave_idx_type ofs;
+  octave_idx_type lastofs;
+  octave_idx_type k;
 
   a += hint;
   lastofs = 0;
@@ -258,7 +260,7 @@
       /* a[hint] < key -- gallop right, until
        * a[hint + lastofs] < key <= a[hint + ofs]
        */
-      const int maxofs = n - hint;	/* &a[n-1] is highest */
+      const octave_idx_type maxofs = n - hint;	/* &a[n-1] is highest */
       while (ofs < maxofs) 
 	{
 	  IFLT (a[ofs], key)
@@ -282,7 +284,7 @@
       /* key <= a[hint] -- gallop left, until
        * a[hint - ofs] < key <= a[hint - lastofs]
        */
-      const int maxofs = hint + 1;	/* &a[0] is lowest */
+      const octave_idx_type maxofs = hint + 1;	/* &a[0] is lowest */
       while (ofs < maxofs) 
 	{
 	  IFLT (*(a-ofs), key)
@@ -309,7 +311,7 @@
   ++lastofs;
   while (lastofs < ofs) 
     {
-      int m = lastofs + ((ofs - lastofs) >> 1);
+      octave_idx_type m = lastofs + ((ofs - lastofs) >> 1);
 
       IFLT (a[m], key)
 	lastofs = m+1;	/* a[m] < key */
@@ -335,12 +337,12 @@
 written as one routine with yet another "left or right?" flag.
 */
 template <class T>
-int
-octave_sort<T>::gallop_right (T key, T *a, int n, int hint)
+octave_idx_type
+octave_sort<T>::gallop_right (T key, T *a, octave_idx_type n, octave_idx_type hint)
 {
-  int ofs;
-  int lastofs;
-  int k;
+  octave_idx_type ofs;
+  octave_idx_type lastofs;
+  octave_idx_type k;
 
   a += hint;
   lastofs = 0;
@@ -350,7 +352,7 @@
       /* key < a[hint] -- gallop left, until
        * a[hint - ofs] <= key < a[hint - lastofs]
        */
-      const int maxofs = hint + 1;	/* &a[0] is lowest */
+      const octave_idx_type maxofs = hint + 1;	/* &a[0] is lowest */
       while (ofs < maxofs) 
 	{
 	  IFLT (key, *(a-ofs))
@@ -375,7 +377,7 @@
       /* a[hint] <= key -- gallop right, until
        * a[hint + lastofs] <= key < a[hint + ofs]
        */
-      const int maxofs = n - hint;	/* &a[n-1] is highest */
+      const octave_idx_type maxofs = n - hint;	/* &a[n-1] is highest */
       while (ofs < maxofs) 
 	{
 	  IFLT (key, a[ofs])
@@ -401,7 +403,7 @@
   ++lastofs;
   while (lastofs < ofs) 
     {
-      int m = lastofs + ((ofs - lastofs) >> 1);
+      octave_idx_type m = lastofs + ((ofs - lastofs) >> 1);
 
       IFLT (key, a[m])
 	ofs = m;	/* key < a[m] */
@@ -437,11 +439,11 @@
   ms.a = 0;
 }
 
-static inline int
-roundupsize (int n)
+static inline octave_idx_type
+roundupsize (octave_idx_type n)
 {
   unsigned int nbits = 3;
-  unsigned int n2 = static_cast<unsigned int> (n) >> 8;
+  octave_idx_type n2 = static_cast<octave_idx_type> (n) >> 8;
 
   /* Round up:
    * If n <       256, to a multiple of        8.
@@ -479,7 +481,7 @@
  */
 template <class T>
 int
-octave_sort<T>::merge_getmem (int need)
+octave_sort<T>::merge_getmem (octave_idx_type need)
 {
   if (need <= ms.alloced)
     return 0;
@@ -510,12 +512,12 @@
  */
 template <class T>
 int
-octave_sort<T>::merge_lo (T *pa, int na, T *pb, int nb)
+octave_sort<T>::merge_lo (T *pa, octave_idx_type na, T *pb, octave_idx_type nb)
 {
-  int k;
+  octave_idx_type k;
   T *dest;
   int result = -1;	/* guilty until proved innocent */
-  int min_gallop = ms.min_gallop;
+  octave_idx_type min_gallop = ms.min_gallop;
 
   if (MERGE_GETMEM (na) < 0)
     return -1;
@@ -532,8 +534,8 @@
 
   for (;;)
     {
-      int acount = 0;	/* # of times A won in a row */
-      int bcount = 0;	/* # of times B won in a row */
+      octave_idx_type acount = 0;	/* # of times A won in a row */
+      octave_idx_type bcount = 0;	/* # of times B won in a row */
 
       /* Do the straightforward thing until (if ever) one run
        * appears to win consistently.
@@ -647,14 +649,14 @@
  */
 template <class T>
 int
-octave_sort<T>::merge_hi (T *pa, int na, T *pb, int nb)
+octave_sort<T>::merge_hi (T *pa, octave_idx_type na, T *pb, octave_idx_type nb)
 {
-  int k;
+  octave_idx_type k;
   T *dest;
   int result = -1;	/* guilty until proved innocent */
   T *basea;
   T *baseb;
-  int min_gallop = ms.min_gallop;
+  octave_idx_type min_gallop = ms.min_gallop;
 
   if (MERGE_GETMEM (nb) < 0)
     return -1;
@@ -674,8 +676,8 @@
 
   for (;;) 
     {
-      int acount = 0;	/* # of times A won in a row */
-      int bcount = 0;	/* # of times B won in a row */
+      octave_idx_type acount = 0;	/* # of times A won in a row */
+      octave_idx_type bcount = 0;	/* # of times B won in a row */
 
       /* Do the straightforward thing until (if ever) one run
        * appears to win consistently.
@@ -787,11 +789,11 @@
  */
 template <class T>
 int
-octave_sort<T>::merge_at (int i)
+octave_sort<T>::merge_at (octave_idx_type i)
 {
   T *pa, *pb;
-  int na, nb;
-  int k;
+  octave_idx_type na, nb;
+  octave_idx_type k;
 
   pa = ms.pending[i].base;
   na = ms.pending[i].len;
@@ -852,7 +854,7 @@
 
   while (ms.n > 1) 
     {
-      int n = ms.n - 2;
+      octave_idx_type n = ms.n - 2;
       if (n > 0 && p[n-1].len <= p[n].len + p[n+1].len) 
 	{
 	  if (p[n-1].len < p[n+1].len)
@@ -885,7 +887,7 @@
 
   while (ms.n > 1) 
     {
-      int n = ms.n - 2;
+      octave_idx_type n = ms.n - 2;
       if (n > 0 && p[n-1].len < p[n+1].len)
 	--n;
       if (merge_at (n) < 0)
@@ -906,10 +908,10 @@
  * See listsort.txt for more info.
  */
 template <class T>
-int
-octave_sort<T>::merge_compute_minrun (int n)
+octave_idx_type
+octave_sort<T>::merge_compute_minrun (octave_idx_type n)
 {
-  int r = 0;	/* becomes 1 if any 1 bits are shifted off */
+  octave_idx_type r = 0;	/* becomes 1 if any 1 bits are shifted off */
 
   while (n >= 64)
     {
@@ -922,7 +924,7 @@
 
 template <class T>
 void
-octave_sort<T>::sort (T *v, int elements)
+octave_sort<T>::sort (T *v, octave_idx_type elements)
 {
   /* Re-initialize the Mergestate as this might be the second time called */
   ms.n = 0;
@@ -930,18 +932,18 @@
 
   if (elements > 1)
     {
-      int nremaining = elements; 
+      octave_idx_type nremaining = elements; 
       T *lo = v;
       T *hi = v + elements;
 
       /* March over the array once, left to right, finding natural runs,
        * and extending short natural runs to minrun elements.
        */
-      int minrun = merge_compute_minrun (nremaining);
+      octave_idx_type minrun = merge_compute_minrun (nremaining);
       do 
 	{
 	  int descending;
-	  int n;
+	  octave_idx_type n;
 
 	  /* Identify next run. */
 	  n = count_run (lo, hi, &descending);
@@ -952,7 +954,7 @@
 	  /* If short, extend to min(minrun, nremaining). */
 	  if (n < minrun) 
 	    {
-	      const int force = nremaining <= minrun ? nremaining : minrun;
+	      const octave_idx_type force = nremaining <= minrun ? nremaining : minrun;
 	      binarysort (lo, lo + force, lo + n);
 	      n = force;
 	    }
diff -r 3c999b2b5de8 -r 402168152bb9 liboctave/oct-sort.h
--- a/liboctave/oct-sort.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/liboctave/oct-sort.h	Thu Jan 31 18:59:11 2008 +0000
@@ -96,6 +96,9 @@
 // Avoid malloc for small temp arrays.
 #define MERGESTATE_TEMP_SIZE 1024
 
+// Enum for type of sort function
+enum sortmode { UNDEFINED, ASCENDING, DESCENDING };
+
 template <class T>
 class
 octave_sort
@@ -110,7 +113,7 @@
 
   void set_compare (bool (*comp) (T, T)) { compare = comp; }
 
-  void sort (T *v, int elements);
+  void sort (T *v, octave_idx_type elements);
 
 private:
 
@@ -125,7 +128,7 @@
   struct s_slice 
   {
     T *base;
-    int len;
+    octave_idx_type len;
   };
   
   struct MergeState 
@@ -134,12 +137,12 @@
     // initialized to MIN_GALLOP.  merge_lo and merge_hi tend to nudge
     // it higher for random data, and lower for highly structured
     // data.
-    int min_gallop;
+    octave_idx_type min_gallop;
 
     // 'a' is temp storage to help with merges.  It contains room for
     // alloced entries.
     T *a;               // may point to temparray below
-    int alloced;
+    octave_idx_type alloced;
     
     // A stack of n pending runs yet to be merged.  Run #i starts at
     // address base[i] and extends for len[i] elements.  It's always
@@ -149,7 +152,7 @@
     //
     // so we could cut the storage for this, but it's a minor amount,
     // and keeping all the info explicit simplifies the code.
-    int n;
+    octave_idx_type n;
     struct s_slice pending[MAX_MERGE_PENDING];
   };
 
@@ -161,31 +164,39 @@
   
   void binarysort (T *lo, T *hi, T *start);
     
-  int count_run (T *lo, T *hi, int *descending);
+  octave_idx_type count_run (T *lo, T *hi, octave_idx_type *descending);
 
-  int gallop_left (T key, T *a, int n, int hint);
+  octave_idx_type gallop_left (T key, T *a, octave_idx_type n, octave_idx_type hint);
 
-  int gallop_right (T key, T *a, int n, int hint);
+  octave_idx_type gallop_right (T key, T *a, octave_idx_type n, octave_idx_type hint);
 
   void merge_init (void);
 
   void merge_freemem (void);
 
-  int merge_getmem (int need);
+  int merge_getmem (octave_idx_type need);
 
-  int merge_lo (T *pa, int na, T *pb, int nb);
+  int merge_lo (T *pa, octave_idx_type na, T *pb, octave_idx_type nb);
 
-  int merge_hi (T *pa, int na, T *pb, int nb);
+  int merge_hi (T *pa, octave_idx_type na, T *pb, octave_idx_type nb);
 
-  int merge_at (int i);
+  int merge_at (octave_idx_type i);
 
   int merge_collapse (void);
 
   int merge_force_collapse (void);
 
-  int merge_compute_minrun (int n);
+  octave_idx_type merge_compute_minrun (octave_idx_type n);
 };
 
+template <class T>
+class
+vec_index
+{
+public:
+  T vec;
+  octave_idx_type indx;
+};
 #endif
 
 /*
diff -r 3c999b2b5de8 -r 402168152bb9 src/ChangeLog
--- a/src/ChangeLog	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ChangeLog	Thu Jan 31 18:59:11 2008 +0000
@@ -1,3 +1,34 @@
+2008-01-31  David Bateman  <dbateman@free.fr>
+
+	* ov.cc (octave_value::octave_value (const ArrayN<bool>&),
+	octave_value::octave_value (const Sparse<bool>&, const MatrixType &),
+	octave_value::octave_value (const ArrayN<std::streamoff>&)): New 
+	constructors.
+	* ov.h: (octave_value (const ArrayN<bool>&),
+	octave_value (const Sparse<bool>&, const MatrixType &),
+	octave_value (const ArrayN<std::streamoff>&)): Declare them.
+	(octave_value sort (octave_idx_type, sortmode) const, octave_value
+	sort (Array<octave_idx_type> &, octave_idx_type, sortmode) const):
+	octave_value sort method.
+	
+	
+	* ov-base.cc (sort): Base versions of teh octave_value sort methods.
+	* ov-base.h (sort): Declare the octave_value sort methods
+	* ov-base-scalar.h (sort): Simple sort methods for scalars.
+	* ov-base-mat.h, ov-base-sparse.h (sort): Sort methods calling
+	underlying array or sparse sort methods.
+	* ov-str-mat.h (sort): String specific sort methods.
+	
+	* TEMPLATE-INST/Array-tc.cc: Instantiate the array sort methods.
+	* ov-streamoff.h (sort): Sort versions returning and error.
+	* oct-stream.cc, ov-typeinfo.cc, Array-os.cc: Null instantiation
+	of array sort methods.
+	
+	* Makefile.in (DLD_XSRC): Remove sort.cc
+	* DLD-FUNCTIONS/sort.cc: Remove
+	* data.cc (Fdata): New function using octave_value sort methods
+	for the sorting. Add tests.
+
 2008-01-30  Thomas Weber  <thomas.weber.mail@gmail.com>
 
 	* pager.cc (Fmore): Doc fix.
diff -r 3c999b2b5de8 -r 402168152bb9 src/DLD-FUNCTIONS/sort.cc
--- a/src/DLD-FUNCTIONS/sort.cc	Wed Jan 30 09:11:58 2008 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1446 +0,0 @@
-/*
-
-Copyright (C) 2004, 2005, 2006, 2007 David Bateman and John W. Eaton
-Copyright (C) 1996, 1997, 1999, 2000, 2001, 2002, 2003 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 <vector>
-
-#include "lo-mappers.h"
-#include "quit.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "lo-ieee.h"
-#include "data-conv.h"
-#include "ov-cx-mat.h"
-#include "ov-cell.h"
-#include "oct-sort.cc"
-
-enum sortmode { UNDEFINED, ASCENDING, DESCENDING };
-
-template <class T>
-class
-vec_index
-{
-public:
-  T vec;
-  octave_idx_type indx;
-};
-
-template <class T>
-bool 
-ascending_compare (T a, T b)
-{
-  return (a < b);
-}
-
-template <class T>
-bool 
-descending_compare (T a, T b)
-{
-  return (a > b);
-}
-
-template <class T>
-bool 
-ascending_compare (vec_index<T> *a, vec_index<T> *b)
-{
-  return (a->vec < b->vec);
-}
-
-template <class T>
-bool 
-descending_compare (vec_index<T> *a, vec_index<T> *b)
-{
-  return (a->vec > b->vec);
-}
-
-template <class T>
-static octave_value
-mx_sort (ArrayN<T> &m, int dim, sortmode mode = UNDEFINED)
-{
-  octave_value retval;
-
-  dim_vector dv = m.dims ();
-
-  if (m.length () < 1)
-    return ArrayN<T> (dv);
-
-  octave_idx_type ns = dv(dim);
-  octave_idx_type iter = dv.numel () / ns;
-  octave_idx_type stride = 1;
-  for (int i = 0; i < dim; i++)
-    stride *= dv(i);
-
-  T *v = m.fortran_vec ();
-  octave_sort<T> sort;
-
-  if (mode == ASCENDING) 
-    sort.set_compare (ascending_compare);
-  else if (mode == DESCENDING)
-    sort.set_compare (descending_compare);
-
-  if (stride == 1)
-    {
-      for (octave_idx_type j = 0; j < iter; j++)
-	{
-	  sort.sort (v, ns);
-	  v += ns;
-	}
-    }
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (T, vi, ns);
-      for (octave_idx_type j = 0; j < iter; j++) 
-	{
-	   octave_idx_type offset = j;
-	   octave_idx_type offset2 = 0;
-	  while (offset >= stride)
-	    {
-	      offset -= stride;
-	      offset2++;
-	    }
-	  offset += offset2 * stride * ns;
-	  
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    vi[i] = v[i*stride + offset];
-
-	  sort.sort (vi, ns);
-	      
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    v[i*stride + offset] = vi[i];
-	}
-    }
-
-  retval = m;
-
-  return retval;
-}
-
-template <class T>
-static octave_value_list
-mx_sort_indexed (ArrayN<T> &m, int dim, sortmode mode = UNDEFINED)
-{
-  octave_value_list retval;
-
-  dim_vector dv = m.dims ();
-
-  if (m.length () < 1)
-    {
-      retval(1) = NDArray (dv);
-      retval(0) = ArrayN<T> (dv);
-      return retval;
-    }
-
-  octave_idx_type ns = dv(dim);
-  octave_idx_type iter = dv.numel () / ns;
-  octave_idx_type stride = 1;
-  for (int i = 0; i < dim; i++)
-    stride *= dv(i);
-
-  T *v = m.fortran_vec ();
-  octave_sort<vec_index<T> *> indexed_sort;
-
-  if (mode == ASCENDING) 
-    indexed_sort.set_compare (ascending_compare);
-  else if (mode == DESCENDING)
-    indexed_sort.set_compare (descending_compare);
-
-  OCTAVE_LOCAL_BUFFER (vec_index<T> *, vi, ns);
-  OCTAVE_LOCAL_BUFFER (vec_index<T>, vix, ns);
-
-  for (octave_idx_type i = 0; i < ns; i++)
-    vi[i] = &vix[i];
-
-  NDArray idx (dv);
-      
-  if (stride == 1)
-    {
-      for (octave_idx_type j = 0; j < iter; j++)
-	{
-	   octave_idx_type offset = j * ns;
-
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    {
-	      vi[i]->vec = v[i];
-	      vi[i]->indx = i + 1;
-	    }
-
-	  indexed_sort.sort (vi, ns);
-
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    {
-	      v[i] = vi[i]->vec;
-	      idx(i + offset) = vi[i]->indx;
-	    }
-	  v += ns;
-	}
-    }
-  else
-    {
-      for (octave_idx_type j = 0; j < iter; j++)
-	{
-	  octave_idx_type offset = j;
-	  octave_idx_type offset2 = 0;
-	  while (offset >= stride)
-	    {
-	      offset -= stride;
-	      offset2++;
-	    }
-	  offset += offset2 * stride * ns;
-	      
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    {
-	      vi[i]->vec = v[i*stride + offset];
-	      vi[i]->indx = i + 1;
-	    }
-
-	  indexed_sort.sort (vi, ns);
-	      
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    {
-	      v[i*stride+offset] = vi[i]->vec;
-	      idx(i*stride+offset) = vi[i]->indx;
-	    }
-	}
-    }
-
-  retval(1) = idx;
-  retval(0) = octave_value (m);
-
-  return retval;
-}
-
-template <class T>
-static octave_value
-mx_sort_sparse (Sparse<T> &m, int dim, sortmode mode = UNDEFINED)
-{
-  octave_value retval;
-
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  if (m.length () < 1)
-    return Sparse<T> (nr, nc);
-
-  if (dim > 0)
-    {
-      m = m.transpose ();
-      nr = m.rows ();
-      nc = m.columns ();
-    }
-
-  octave_sort<T> sort;
-
-  if (mode == ASCENDING) 
-    sort.set_compare (ascending_compare);
-  else if (mode == DESCENDING)
-    sort.set_compare (descending_compare);
-
-  T *v = m.data ();
-  octave_idx_type *cidx = m.cidx ();
-  octave_idx_type *ridx = m.ridx ();
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      octave_idx_type ns = cidx [j + 1] - cidx [j];
-      sort.sort (v, ns);
-
-      octave_idx_type i;
-      if (mode == ASCENDING) 
-	{
-	  for (i = 0; i < ns; i++)
-	    if (ascending_compare (static_cast<T> (0), v [i]))
-	      break;
-	}
-      else
-	{
-	  for (i = 0; i < ns; i++)
-	    if (descending_compare (static_cast<T> (0), v [i]))
-	      break;
-	}
-      for (octave_idx_type k = 0; k < i; k++)
-	ridx [k] = k;
-      for (octave_idx_type k = i; k < ns; k++)
-	ridx [k] = k - ns + nr; 
-
-      v += ns;
-      ridx += ns;
-    }
-
-  if (dim > 0)
-      m = m.transpose ();
-
-  retval = m;
-
-  return retval;
-}
-
-template <class T>
-static octave_value_list
-mx_sort_sparse_indexed (Sparse<T> &m, int dim, sortmode mode = UNDEFINED)
-{
-  octave_value_list retval;
-
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  if (m.length () < 1)
-    {
-      retval (1) = NDArray (dim_vector (nr, nc));
-      retval (0) = octave_value (SparseMatrix (nr, nc));
-      return retval;
-    }
-
-  if (dim > 0)
-    {
-      m = m.transpose ();
-      nr = m.rows ();
-      nc = m.columns ();
-    }
-
-  octave_sort<vec_index<T> *> indexed_sort;
-
-  if (mode == ASCENDING) 
-    indexed_sort.set_compare (ascending_compare);
-  else if (mode == DESCENDING)
-    indexed_sort.set_compare (descending_compare);
-
-  T *v = m.data ();
-  octave_idx_type *cidx = m.cidx ();
-  octave_idx_type *ridx = m.ridx ();
-
-  OCTAVE_LOCAL_BUFFER (vec_index<T> *, vi, nr);
-  OCTAVE_LOCAL_BUFFER (vec_index<T>, vix, nr);
-
-  for (octave_idx_type i = 0; i < nr; i++)
-    vi[i] = &vix[i];
-
-  Matrix idx (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      octave_idx_type ns = cidx [j + 1] - cidx [j];
-      octave_idx_type offset = j * nr;
-
-      if (ns == 0)
-	{
-	  for (octave_idx_type k = 0; k < nr; k++)
-	    idx (offset + k) = k + 1;
-	}
-      else
-	{
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    {
-	      vi[i]->vec = v[i];
-	      vi[i]->indx = ridx[i] + 1;
-	    }
-
-	  indexed_sort.sort (vi, ns);
-
-	  octave_idx_type i;
-	  if (mode == ASCENDING) 
-	    {
-	      for (i = 0; i < ns; i++)
-		if (ascending_compare (static_cast<T> (0), vi [i] -> vec))
-		  break;
-	    }
-	  else
-	    {
-	      for (i = 0; i < ns; i++)
-		if (descending_compare (static_cast<T> (0), vi [i] -> vec))
-		  break;
-	    }
-
-	  octave_idx_type ii = 0;
-	  octave_idx_type jj = i;
-	  for (octave_idx_type k = 0; k < nr; k++)
-	    {
-	      if (ii < ns && ridx[ii] == k)
-		ii++;
-	      else
-		idx (offset + jj++) = k + 1;
-	    }
-
-	  for (octave_idx_type k = 0; k < i; k++)
-	    {
-	      v [k] = vi [k] -> vec;
-	      idx (k + offset) = vi [k] -> indx;
-	      ridx [k] = k;
-	    }
-
-	  for (octave_idx_type k = i; k < ns; k++)
-	    {
-	      v [k] = vi [k] -> vec;
-	      idx (k - ns + nr + offset) = vi [k] -> indx;
-	      ridx [k] = k - ns + nr; 
-	    }
-
-	  v += ns;
-	  ridx += ns;
-	}
-    }
-
-  if (dim > 0)
-    {
-      m = m.transpose ();
-      idx = idx.transpose ();
-    }
-
-  retval (1) = idx;
-  retval(0) = octave_value (m);
-
-  return retval;
-}
-
-// If we have IEEE 754 data format, then we can use the trick of
-// casting doubles as unsigned eight byte integers, and with a little
-// bit of magic we can automatically sort the NaN's correctly.
-
-#if defined (HAVE_IEEE754_DATA_FORMAT)
-
-static inline uint64_t
-FloatFlip (uint64_t f)
-{
-  uint64_t mask
-    = -static_cast<int64_t>(f >> 63) | 0x8000000000000000ULL;
-
-  return f ^ mask;
-}
-
-static inline uint64_t
-IFloatFlip (uint64_t f)
-{
-  uint64_t mask = ((f >> 63) - 1) | 0x8000000000000000ULL;
-
-  return f ^ mask;
-}
-
-template <>
-bool
-ascending_compare (uint64_t a, 
-		   uint64_t b)
-{
-  return (a < b);
-}
-
-template <>
-bool
-ascending_compare (vec_index<uint64_t> *a, 
-		   vec_index<uint64_t> *b)
-{
-  return (a->vec < b->vec);
-}
-
-template <>
-bool
-descending_compare (uint64_t a, 
-		    uint64_t b)
-{
-  return (a > b);
-}
-
-template <>
-bool
-descending_compare (vec_index<uint64_t> *a, 
-		    vec_index<uint64_t> *b)
-{
-  return (a->vec > b->vec);
-}
-
-template class octave_sort<uint64_t>;
-template class vec_index<uint64_t>;
-template class octave_sort<vec_index<uint64_t> *>;
-
-template <>
-octave_value
-mx_sort (ArrayN<double> &m, int dim, sortmode mode)
-{
-  octave_value retval;
-
-  dim_vector dv = m.dims ();
-
-  if (m.length () < 1)
-    return ArrayN<double> (dv);
-
-  octave_idx_type ns = dv(dim);
-  octave_idx_type iter = dv.numel () / ns;
-  octave_idx_type stride = 1;
-  for (int i = 0; i < dim; i++)
-    stride *= dv(i);
-
-  double *v = m.fortran_vec ();
-
-  uint64_t *p = reinterpret_cast<uint64_t *> (v);
-
-  octave_sort<uint64_t> sort;
-
-  if (mode == ASCENDING)
-    sort.set_compare (ascending_compare);
-  else if (mode == DESCENDING)
-    sort.set_compare (descending_compare);
-
-  if (stride == 1)
-    {
-      for (octave_idx_type j = 0; j < iter; j++)
-	{
-	  // Flip the data in the vector so that int compares on
-	  // IEEE754 give the correct ordering.
-
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    p[i] = FloatFlip (p[i]);
-	      
-	  sort.sort (p, ns);
-
-	  // Flip the data out of the vector so that int compares
-	  // on IEEE754 give the correct ordering.
-
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    p[i] = IFloatFlip (p[i]);
-
-	  // There are two representations of NaN.  One will be
-	  // sorted to the beginning of the vector and the other
-	  // to the end.  If it will be sorted incorrectly, fix
-	  // things up.
-
-	  if (lo_ieee_signbit (octave_NaN))
-	    {
-	      if (mode == UNDEFINED || mode == ASCENDING)
-		{
-		  octave_idx_type i = 0;
-		  double *vtmp = reinterpret_cast<double *> (p);
-		  while (xisnan (vtmp[i++]) && i < ns);
-		  for (octave_idx_type l = 0; l < ns - i + 1; l++)
-		    vtmp[l] = vtmp[l+i-1];
-		  for (octave_idx_type l = ns - i + 1; l < ns; l++)
-		    vtmp[l] = octave_NaN;
-		}
-	      else
-		{
-		  octave_idx_type i = ns;
-		  double *vtmp = reinterpret_cast<double *> (p);
-		  while (xisnan (vtmp[--i]) && i > 0);
-		  for (octave_idx_type l = i; l >= 0; l--)
-		    vtmp[l-i+ns-1] = vtmp[l];
-		  for (octave_idx_type l = 0; l < ns - i - 1; l++)
-		    vtmp[l] = octave_NaN;
-		}
-	    }
-
-	  p += ns;
-	}
-    }
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (uint64_t, vi, ns);
-
-      for (octave_idx_type j = 0; j < iter; j++)
-	{
-	  octave_idx_type offset = j;
-	  octave_idx_type offset2 = 0;
-	  while (offset >= stride)
-	    {
-	      offset -= stride;
-	      offset2++;
-	    }
-	  offset += offset2 * stride * ns;
-
-	  // Flip the data in the vector so that int compares on
-	  // IEEE754 give the correct ordering.
-
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    vi[i] = FloatFlip (p[i*stride + offset]);
-
-	  sort.sort (vi, ns);
-
-	  // Flip the data out of the vector so that int compares
-	  // on IEEE754 give the correct ordering.
-
-	  for (octave_idx_type i = 0; i < ns; i++)
-	    p[i*stride + offset] = IFloatFlip (vi[i]);
-	      
-	  // There are two representations of NaN. One will be
-	  // sorted to the beginning of the vector and the other
-	  // to the end. If it will be sorted to the beginning,
-	  // fix things up.
-
-	  if (lo_ieee_signbit (octave_NaN))
-	    {
-	      if (mode == UNDEFINED || mode == ASCENDING)
-		{
-		   octave_idx_type i = 0;
-		  while (xisnan (v[i++*stride + offset]) && i < ns);
-		  for (octave_idx_type l = 0; l < ns - i + 1; l++)
-		    v[l*stride + offset] = v[(l+i-1)*stride + offset];
-		  for (octave_idx_type l = ns - i + 1; l < ns; l++)
-		    v[l*stride + offset] = octave_NaN;
-		}
-	      else
-		{
-		   octave_idx_type i = ns;
-		  while (xisnan (v[--i*stride + offset]) && i > 0);
-		  for (octave_idx_type l = i; l >= 0; l--)
-		    v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
-		  for (octave_idx_type l = 0; l < ns - i - 1; l++)
-		    v[l*stride + offset] = octave_NaN;
-		}
-	    }
-	}
-    }
-
-  retval = m;
-
-  return retval;
-}
-
-// Should other overloaded functions have their static keywords removed?
-template <>
-octave_value_list
-mx_sort_indexed (ArrayN<double> &m, int dim, sortmode mode)
-{
-  octave_value_list retval;
-
-  dim_vector dv = m.dims ();
-
-  if (m.length () < 1)
-    {
-      retval(1) = NDArray (dv);
-      retval(0) = ArrayN<double> (dv);
-      return retval;
-    }
-
-  octave_idx_type ns = dv(dim);
-  octave_idx_type iter = dv.numel () / ns;
-  octave_idx_type stride = 1;
-  for (int i = 0; i < dim; i++)
-    stride *= dv(i);
-
-  double *v = m.fortran_vec ();
-
-  uint64_t *p = reinterpret_cast<uint64_t *> (v);
-
-  octave_sort<vec_index<uint64_t> *> indexed_sort;
-
-  if (mode == ASCENDING)
-    indexed_sort.set_compare (ascending_compare);
-  else if (mode == DESCENDING)
-    indexed_sort.set_compare (descending_compare);
-
-  OCTAVE_LOCAL_BUFFER (vec_index<uint64_t> *, vi, ns);
-  OCTAVE_LOCAL_BUFFER (vec_index<uint64_t>, vix, ns);
-  
-  for (octave_idx_type i = 0; i < ns; i++)
-    vi[i] = &vix[i];
-
-  NDArray idx (dv);
-      
-  for (octave_idx_type j = 0; j < iter; j++)
-    {
-      octave_idx_type offset = j;
-      octave_idx_type offset2 = 0;
-      while (offset >= stride)
-	{
-	  offset -= stride;
-	  offset2++;
-	}
-      offset += offset2 * stride * ns;
-
-      // Flip the data in the vector so that int compares on
-      // IEEE754 give the correct ordering.
-
-      for (octave_idx_type i = 0; i < ns; i++)
-	{
-	  vi[i]->vec = FloatFlip (p[i*stride + offset]);
-	  vi[i]->indx = i + 1;
-	}
-
-      indexed_sort.sort (vi, ns);
-
-      // Flip the data out of the vector so that int compares on
-      // IEEE754 give the correct ordering
-
-      for (octave_idx_type i = 0; i < ns; i++)
-	{
-	  p[i*stride + offset] = IFloatFlip (vi[i]->vec);
-	  idx(i*stride + offset) = vi[i]->indx;
-	}
-
-      // There are two representations of NaN.  One will be sorted
-      // to the beginning of the vector and the other to the end.
-      // If it will be sorted to the beginning, fix things up.
-
-      if (lo_ieee_signbit (octave_NaN))
-	{
-	  if (mode == UNDEFINED || mode == ASCENDING)
-	    {
-	      octave_idx_type i = 0;
-	      while (xisnan (v[i++*stride+offset]) && i < ns);
-	      OCTAVE_LOCAL_BUFFER (double, itmp, i - 1);
-	      for (octave_idx_type l = 0; l < i -1; l++)
-		itmp[l] = idx(l*stride + offset);
-	      for (octave_idx_type l = 0; l < ns - i + 1; l++)
-		{
-		  v[l*stride + offset] = v[(l+i-1)*stride + offset];
-		  idx(l*stride + offset) = idx((l+i-1)*stride + offset);
-		}
-	      for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++)
-		{
-		  v[l*stride + offset] = octave_NaN;
-		  idx(l*stride + offset) = itmp[k];
-		}
-	    }
-	  else 
-	    {
-	      octave_idx_type i = ns;
-	      while (xisnan (v[--i*stride+offset]) && i > 0);
-	      OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1);
-	      for (octave_idx_type l = 0; l < ns - i -1; l++)
-		itmp[l] = idx((l+i+1)*stride + offset);
-	      for (octave_idx_type l = i; l >= 0; l--)
-		{
-		  v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
-		  idx((l-i+ns-1)*stride + offset) = idx(l*stride + offset);
-		}
-	      for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++)
-		{
-		  v[l*stride + offset] = octave_NaN;
-		  idx(l*stride + offset) = itmp[k];
-		}
-	    }
-	}
-    }
-
-  retval(1) = idx;
-  retval(0) = m;
-
-  return retval;
-}
-
-#else
-
-template <>
-bool
-ascending_compare (double a, double b)
-{
-  return (xisnan (b) || (a < b));
-}
-
-template <>
-bool
-ascending_compare (vec_index<double> *a, vec_index<double> *b)
-{
-  return (xisnan (b->vec) || (a->vec < b->vec));
-}
-
-template <>
-bool
-descending_compare (double a, double b)
-{
-  return (xisnan (a) || (a > b));
-}
-
-template <>
-bool
-descending_compare (vec_index<double> *a, vec_index<double> *b)
-{
-  return (xisnan (a->vec) || (a->vec > b->vec));
-}
-
-template class octave_sort<double>;
-template class vec_index<double>;
-template class octave_sort<vec_index<double> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-static octave_value_list
-mx_sort (ArrayN<double> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<double> &m, int dim, sortmode mode);
-#endif
-#endif
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-static octave_value_list
-mx_sort_sparse (Sparse<double> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_sparse_indexed (Sparse<double> &m, int dim, sortmode mode);
-#endif
-
-// std::abs(Inf) returns NaN!!
-static inline double
-xabs (const Complex& x)
-{
-  return (xisinf (x.real ()) || xisinf (x.imag ())) ? octave_Inf : abs (x);
-}
-
-template <>
-bool
-ascending_compare (Complex a, Complex b)
-{
-  return (xisnan (b) || (xabs (a) < xabs (b))
-	  || ((xabs (a) == xabs (b)) && (arg (a) < arg (b))));
-}
-
-bool
-operator < (const Complex& a, const Complex& b)
-{
-  return (xisnan (b) || (xabs (a) < xabs (b))
-	  || ((xabs (a) == xabs (b)) && (arg (a) < arg (b))));
-}
-
-template <>
-bool
-descending_compare (Complex a, Complex b)
-{
-  return (xisnan (a) || (xabs (a) > xabs (b))
-	  || ((xabs (a) == xabs (b)) && (arg (a) > arg (b))));
-}
-
-bool
-operator > (const Complex& a, const Complex& b)
-{
-  return (xisnan (a) || (xabs (a) > xabs (b))
-	  || ((xabs (a) == xabs (b)) && (arg (a) > arg (b))));
-}
-
-template <>
-bool
-ascending_compare (vec_index<Complex> *a, vec_index<Complex> *b)
-{
-  return (xisnan (b->vec)
-	  || (xabs (a->vec) < xabs (b->vec))
-	  || ((xabs (a->vec) == xabs (b->vec))
-	      && (arg (a->vec) < arg (b->vec))));
-}
-
-template <>
-bool
-descending_compare (vec_index<Complex> *a, vec_index<Complex> *b)
-{
-  return (xisnan (a->vec)
-	  || (xabs (a->vec) > xabs (b->vec))
-	  || ((xabs (a->vec) == xabs (b->vec))
-	      && (arg (a->vec) > arg (b->vec))));
-}
-
-template class octave_sort<Complex>;
-template class vec_index<Complex>;
-template class octave_sort<vec_index<Complex> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-static octave_value_list
-mx_sort (ArrayN<Complex> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<Complex> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_sparse (Sparse<Complex> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_sparse_indexed (Sparse<Complex> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<char>;
-template class vec_index<char>;
-template class octave_sort<vec_index<char> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (char a, char b);
-
-bool
-ascending_compare (vec_index<char> *a, vec_index<char> *b);
-
-bool
-descending_compare (char a, char b);
-
-bool
-descending_compare (vec_index<char> *a, vec_index<char> *b);
-
-static octave_value_list
-mx_sort (ArrayN<char> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<char> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_int8>;
-template class vec_index<octave_int8>;
-template class octave_sort<vec_index<octave_int8> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_int8 a, octave_int8 b);
-
-bool
-ascending_compare (vec_index<octave_int8> *a, vec_index<octave_int8> *b);
-
-bool
-descending_compare (octave_int8 a, octave_int8 b);
-
-bool
-descending_compare (vec_index<octave_int8> *a, vec_index<octave_int8> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_int8> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_int8> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_uint8>;
-template class vec_index<octave_uint8>;
-template class octave_sort<vec_index<octave_uint8> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_uint8 a, octave_uint8 b);
-
-bool
-ascending_compare (vec_index<octave_uint8> *a, vec_index<octave_uint8> *b);
-
-bool
-descending_compare (octave_uint8 a, octave_uint8 b);
-
-bool
-descending_compare (vec_index<octave_uint8> *a, vec_index<octave_uint8> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_uint8> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_uint8> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_int16>;
-template class vec_index<octave_int16>;
-template class octave_sort<vec_index<octave_int16> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_int16 a, octave_int16 b);
-
-bool
-ascending_compare (vec_index<octave_int16> *a, vec_index<octave_int16> *b);
-
-bool
-descending_compare (octave_int16 a, octave_int16 b);
-
-bool
-descending_compare (vec_index<octave_int16> *a, vec_index<octave_int16> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_int16> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_int16> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_uint16>;
-template class vec_index<octave_uint16>;
-template class octave_sort<vec_index<octave_uint16> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_uint16 a, octave_uint16 b);
-
-bool
-ascending_compare (vec_index<octave_uint16> *a, vec_index<octave_uint16> *b);
-
-bool
-descending_compare (octave_uint16 a, octave_uint16 b);
-
-bool
-descending_compare (vec_index<octave_uint16> *a, vec_index<octave_uint16> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_uint16> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_uint16> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_int32>;
-template class vec_index<octave_int32>;
-template class octave_sort<vec_index<octave_int32> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_int32 a, octave_int32 b);
-
-bool
-ascending_compare (vec_index<octave_int32> *a, vec_index<octave_int32> *b);
-
-bool
-descending_compare (octave_int32 a, octave_int32 b);
-
-bool
-descending_compare (vec_index<octave_int32> *a, vec_index<octave_int32> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_int32> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_int32> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_uint32>;
-template class vec_index<octave_uint32>;
-template class octave_sort<vec_index<octave_uint32> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_uint32 a, octave_uint32 b);
-
-bool
-ascending_compare (vec_index<octave_uint32> *a, vec_index<octave_uint32> *b);
-
-bool
-descending_compare (octave_uint32 a, octave_uint32 b);
-
-bool
-descending_compare (vec_index<octave_uint32> *a, vec_index<octave_uint32> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_uint32> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_uint32> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_int64>;
-template class vec_index<octave_int64>;
-template class octave_sort<vec_index<octave_int64> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_int64 a, octave_int64 b);
-
-bool
-ascending_compare (vec_index<octave_int64> *a, vec_index<octave_int64> *b);
-
-bool
-descending_compare (octave_int64 a, octave_int64 b);
-
-bool
-descending_compare (vec_index<octave_int64> *a, vec_index<octave_int64> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_int64> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_int64> &m, int dim, sortmode mode);
-#endif
-
-template class octave_sort<octave_uint64>;
-template class vec_index<octave_uint64>;
-template class octave_sort<vec_index<octave_uint64> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-bool
-ascending_compare (octave_uint64 a, octave_uint64 b);
-
-bool
-ascending_compare (vec_index<octave_uint64> *a, vec_index<octave_uint64> *b);
-
-bool
-descending_compare (octave_uint64 a, octave_uint64 b);
-
-bool
-descending_compare (vec_index<octave_uint64> *a, vec_index<octave_uint64> *b);
-
-static octave_value_list
-mx_sort (ArrayN<octave_uint64> &m, int dim, sortmode mode);
-
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_uint64> &m, int dim, sortmode mode);
-#endif
-
-template <>
-bool
-ascending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b)
-{
-  return (a->vec.string_value () < b->vec.string_value ());
-}
-
-template <>
-bool
-descending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b)
-{
-  return (a->vec.string_value () > b->vec.string_value ());
-}
-
-template class vec_index<octave_value>;
-template class octave_sort<vec_index<octave_value> *>;
-
-#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
-static octave_value_list
-mx_sort_indexed (ArrayN<octave_value> &m, int dim, sortmode mode);
-#endif
-
-DEFUN_DLD (sort, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x})\n\
-@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim})\n\
-@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{mode})\n\
-@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim}, @var{mode})\n\
-Return a copy of @var{x} with the elements arranged in increasing\n\
-order.  For matrices, @code{sort} orders the elements in each column.\n\
-\n\
-For example,\n\
-\n\
-@example\n\
-@group\n\
-sort ([1, 2; 2, 3; 3, 1])\n\
-     @result{}  1  1\n\
-         2  2\n\
-         3  3\n\
-@end group\n\
-@end example\n\
-\n\
-The @code{sort} function may also be used to produce a matrix\n\
-containing the original row indices of the elements in the sorted\n\
-matrix.  For example,\n\
-\n\
-@example\n\
-@group\n\
-[s, i] = sort ([1, 2; 2, 3; 3, 1])\n\
-     @result{} s = 1  1\n\
-            2  2\n\
-            3  3\n\
-     @result{} i = 1  3\n\
-            2  1\n\
-            3  2\n\
-@end group\n\
-@end example\n\
-\n\
-If the optional argument @var{dim} is given, then the matrix is sorted\n\
-along the dimension defined by @var{dim}. The optional argument @code{mode}\n\
-defines the order in which the values will be sorted. Valid values of\n\
-@code{mode} are `ascend' or `descend'.\n\
-\n\
-For equal elements, the indices are such that the equal elements are listed\n\
-in the order that appeared in the original list.\n\
-\n\
-The @code{sort} function may also be used to sort strings and cell arrays\n\
-of strings, in which case the dictionary order of the strings is used.\n\
-\n\
-The algorithm used in @code{sort} is optimized for the sorting of partially\n\
-ordered lists.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-  sortmode smode = ASCENDING;
-
-  if (nargin < 1 || nargin > 3)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  bool return_idx = nargout > 1;
-
-  octave_value arg = args(0);
-
-  int dim = 0;
-  if (nargin > 1)
-    {
-      if (args(1).is_string ())
-	{
-	  std::string mode = args(1).string_value();
-	  if (mode == "ascend")
-	    smode = ASCENDING;
-	  else if (mode == "descend")
-	    smode = DESCENDING;
-	  else
-	    {
-	      error ("sort: mode must be either \"ascend\" or \"descend\"");
-	      return retval;
-	    }
-	}
-      else
-	dim = args(1).nint_value () - 1;
-    }
-
-  if (nargin > 2)
-    {
-      if (args(1).is_string ())
-	{
-	  print_usage ();
-	  return retval;
-	}
-
-      if (! args(2).is_string ())
-	{
-	  error ("sort: mode must be a string");
-	  return retval;
-	}
-      std::string mode = args(2).string_value();
-      if (mode == "ascend")
-	smode = ASCENDING;
-      else if (mode == "descend")
-	smode = DESCENDING;
-      else
-	{
-	  error ("sort: mode must be either \"ascend\" or \"descend\"");
-	  return retval;
-	}
-    }
-
-  dim_vector dv = arg.dims ();
-  if (error_state)
-    {
-      gripe_wrong_type_arg ("sort", arg);
-      return retval;
-    }
-  if (nargin == 1 || args(1).is_string ())
-    {
-      // Find first non singleton dimension
-      for (int i = 0; i < dv.length (); i++)
-	if (dv(i) > 1)
-	  {
-	    dim = i;
-	    break;
-	  }
-    }
-  else
-    {
-      if (dim < 0 || dim > dv.length () - 1)
-	{
-	  error ("sort: dim must be a valid dimension");
-	  return retval;
-	}
-    }
-
-  // FIXME -- Perhaps sort should be made a method of the octave_value
-  // classes and then the mess of if statements below might be
-  // replaced with
-  //
-  //   retval = arg.sort (dim, smode, return_idx);
-
-  if (arg.is_real_type ())
-    {
-      if (arg.is_sparse_type ())
-	{
-	  SparseMatrix m = arg.sparse_matrix_value ();
-
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_sparse_indexed (m, dim, smode);
-	      else
-		retval = mx_sort_sparse (m, dim, smode);
-	    }
-	}
-      else if (arg.is_int8_type ())
-	{
-	  int8NDArray m = arg.int8_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_uint8_type ())
-	{
-	  uint8NDArray m = arg.uint8_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_int16_type ())
-	{
-	  int16NDArray m = arg.int16_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_uint16_type ())
-	{
-	  uint16NDArray m = arg.uint16_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_int32_type ())
-	{
-	  int32NDArray m = arg.int32_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_uint32_type ())
-	{
-	  uint32NDArray m = arg.uint32_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_int64_type ())
-	{
-	  int64NDArray m = arg.int64_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else if (arg.is_uint64_type ())
-	{
-	  uint64NDArray m = arg.uint64_array_value ();
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (m, dim, smode);
-	      else
-		retval = mx_sort (m, dim, smode);
-	    }
-	}
-      else
-	{
-	  NDArray m = arg.array_value ();
-
-	  if (! error_state)
-	    {
-#ifdef HAVE_IEEE754_DATA_FORMAT
-	      // As operator > gives the right result, can special case here
-	      if (! return_idx && smode == ASCENDING)
-		retval = mx_sort (m, dim);
-	      else
-#endif
-		{
-		  if (return_idx)
-		    retval = mx_sort_indexed (m, dim, smode);
-		  else
-		    retval = mx_sort (m, dim, smode);
-		}
-	    }
-	}
-    }
-  else if (arg.is_complex_type ())
-    {
-      if (arg.is_sparse_type ())
-	{
-	  SparseComplexMatrix cm = arg.sparse_complex_matrix_value ();
-
-	  if (! error_state)
-	    {
-	      if (return_idx)
-		retval = mx_sort_sparse_indexed (cm, dim, smode);
-	      else
-		retval = mx_sort_sparse (cm, dim, smode);
-	    }
-	}
-      else
-	{
-	  ComplexNDArray cm = arg.complex_array_value ();
-
-	  if (! error_state)
-	    {
-	      // The indexed version seems to be slightly faster
-	      retval = mx_sort_indexed (cm, dim, smode);
-	    }
-	}
-    }
-  else if (arg.is_string ())
-    {
-      charNDArray chm = arg.char_array_value ();
-
-      if (! error_state)
-	{
-	  // As operator > gives the right result, can special case here
-	  if (! return_idx && smode == ASCENDING)
-	    retval = mx_sort (chm, dim);
-	  else
-	    {
-	      if (return_idx)
-		retval = mx_sort_indexed (chm, dim, smode);
-	      else
-		retval = mx_sort (chm, dim, smode);
-	    }
-
-	  // FIXME It would have been better to call 
-	  // "octave_value(m, true)" but how can that be done 
-	  // within the template
-	  retval(0) = retval(0).convert_to_str (false, true);
-	}
-    }
-  else if (arg.is_cell ())
-    {
-      Cell cellm = arg.cell_value ();
-
-      // Need to check that all elements are strings
-      for (octave_idx_type i = 0; i < cellm.numel (); i++)
-	if (! cellm(i).is_string ())
-	  {
-	    gripe_wrong_type_arg ("sort", arg);
-	    break;
-	  }
-
-      // Don't have unindexed version as ">" operator doesn't return bool
-      if (!error_state)
-	retval = mx_sort_indexed (cellm, dim, smode);
-    }
-  else
-    gripe_wrong_type_arg ("sort", arg);
-
-  return retval;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
diff -r 3c999b2b5de8 -r 402168152bb9 src/Makefile.in
--- a/src/Makefile.in	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/Makefile.in	Thu Jan 31 18:59:11 2008 +0000
@@ -67,7 +67,7 @@
 	gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
 	givens.cc hess.cc inv.cc kron.cc lsode.cc \
 	lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \
-	quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
+	quad.cc qz.cc rand.cc regexp.cc schur.cc sparse.cc \
 	spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \
 	sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc typecast.cc \
 	urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \
diff -r 3c999b2b5de8 -r 402168152bb9 src/TEMPLATE-INST/Array-os.cc
--- a/src/TEMPLATE-INST/Array-os.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/TEMPLATE-INST/Array-os.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -34,10 +34,13 @@
 typedef scanf_format_elt* scanf_format_elt_ptr;
 typedef printf_format_elt* printf_format_elt_ptr;
 
+NO_INSTANTIATE_ARRAY_SORT (scanf_format_elt_ptr);
 INSTANTIATE_ARRAY (scanf_format_elt_ptr, OCTINTERP_API);
 
+NO_INSTANTIATE_ARRAY_SORT (printf_format_elt_ptr);
 INSTANTIATE_ARRAY (printf_format_elt_ptr, OCTINTERP_API);
 
+NO_INSTANTIATE_ARRAY_SORT (octave_stream);
 INSTANTIATE_ARRAY (octave_stream, OCTINTERP_API);
 
 /*
diff -r 3c999b2b5de8 -r 402168152bb9 src/TEMPLATE-INST/Array-tc.cc
--- a/src/TEMPLATE-INST/Array-tc.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/TEMPLATE-INST/Array-tc.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -40,6 +40,40 @@
 
 #include "oct-obj.h"
 
+#define IFLT(a, b) if (compare ? compare ((a), (b)) : true)
+
+#include "oct-sort.cc"
+
+template <>
+bool
+ascending_compare (octave_value a, octave_value b)
+{
+  return (a.string_value () < b.string_value ());
+}
+
+template <>
+bool
+ascending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b)
+{
+  return (a->vec.string_value () < b->vec.string_value ());
+}
+
+template <>
+bool
+descending_compare (octave_value a, octave_value b)
+{
+  return (a.string_value () > b.string_value ());
+}
+
+template <>
+bool
+descending_compare (vec_index<octave_value> *a, vec_index<octave_value> *b)
+{
+  return (a->vec.string_value () > b->vec.string_value ());
+}
+
+INSTANTIATE_ARRAY_SORT (octave_value);
+
 template class OCTINTERP_API Array<octave_value>;
 
 INSTANTIATE_ARRAY_ASSIGN (octave_value, octave_value, OCTINTERP_API);
diff -r 3c999b2b5de8 -r 402168152bb9 src/data.cc
--- a/src/data.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/data.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -3305,6 +3305,300 @@
   return retval;
 }
 
+DEFUN (sort, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x})\n\
+@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim})\n\
+@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{mode})\n\
+@deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim}, @var{mode})\n\
+Return a copy of @var{x} with the elements arranged in increasing\n\
+order.  For matrices, @code{sort} orders the elements in each column.\n\
+\n\
+For example,\n\
+\n\
+@example\n\
+@group\n\
+sort ([1, 2; 2, 3; 3, 1])\n\
+     @result{}  1  1\n\
+         2  2\n\
+         3  3\n\
+@end group\n\
+@end example\n\
+\n\
+The @code{sort} function may also be used to produce a matrix\n\
+containing the original row indices of the elements in the sorted\n\
+matrix.  For example,\n\
+\n\
+@example\n\
+@group\n\
+[s, i] = sort ([1, 2; 2, 3; 3, 1])\n\
+     @result{} s = 1  1\n\
+            2  2\n\
+            3  3\n\
+     @result{} i = 1  3\n\
+            2  1\n\
+            3  2\n\
+@end group\n\
+@end example\n\
+\n\
+If the optional argument @var{dim} is given, then the matrix is sorted\n\
+along the dimension defined by @var{dim}. The optional argument @code{mode}\n\
+defines the order in which the values will be sorted. Valid values of\n\
+@code{mode} are `ascend' or `descend'.\n\
+\n\
+For equal elements, the indices are such that the equal elements are listed\n\
+in the order that appeared in the original list.\n\
+\n\
+The @code{sort} function may also be used to sort strings and cell arrays\n\
+of strings, in which case the dictionary order of the strings is used.\n\
+\n\
+The algorithm used in @code{sort} is optimized for the sorting of partially\n\
+ordered lists.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+  sortmode smode = ASCENDING;
+
+  if (nargin < 1 || nargin > 3)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  bool return_idx = nargout > 1;
+
+  octave_value arg = args(0);
+
+  int dim = 0;
+  if (nargin > 1)
+    {
+      if (args(1).is_string ())
+	{
+	  std::string mode = args(1).string_value();
+	  if (mode == "ascend")
+	    smode = ASCENDING;
+	  else if (mode == "descend")
+	    smode = DESCENDING;
+	  else
+	    {
+	      error ("sort: mode must be either \"ascend\" or \"descend\"");
+	      return retval;
+	    }
+	}
+      else
+	dim = args(1).nint_value () - 1;
+    }
+
+  if (nargin > 2)
+    {
+      if (args(1).is_string ())
+	{
+	  print_usage ();
+	  return retval;
+	}
+
+      if (! args(2).is_string ())
+	{
+	  error ("sort: mode must be a string");
+	  return retval;
+	}
+      std::string mode = args(2).string_value();
+      if (mode == "ascend")
+	smode = ASCENDING;
+      else if (mode == "descend")
+	smode = DESCENDING;
+      else
+	{
+	  error ("sort: mode must be either \"ascend\" or \"descend\"");
+	  return retval;
+	}
+    }
+
+  dim_vector dv = arg.dims ();
+  if (error_state)
+    {
+      gripe_wrong_type_arg ("sort", arg);
+      return retval;
+    }
+  if (nargin == 1 || args(1).is_string ())
+    {
+      // Find first non singleton dimension
+      for (int i = 0; i < dv.length (); i++)
+	if (dv(i) > 1)
+	  {
+	    dim = i;
+	    break;
+	  }
+    }
+  else
+    {
+      if (dim < 0 || dim > dv.length () - 1)
+	{
+	  error ("sort: dim must be a valid dimension");
+	  return retval;
+	}
+    }
+
+  if (return_idx)
+    {
+      Array<octave_idx_type> sidx;
+
+      retval (0) = arg.sort (sidx, dim, smode);
+
+      octave_idx_type *ps = sidx.fortran_vec ();
+      NDArray midx (sidx.dims ());
+      double *pm = midx.fortran_vec ();
+
+      for (octave_idx_type i = 0; i < sidx.numel (); i++)
+	pm [i] = static_cast<double> 
+	  (ps [i] + static_cast<octave_idx_type> (1));
+
+      retval (1) = midx;
+    }
+  else
+    retval(0) = arg.sort (dim, smode);
+
+  return retval;
+}
+
+/*
+
+%% Double
+%!assert (sort ([NaN, 1, -1, 2, Inf]), [-1, 1, 2, Inf, NaN])
+%!assert (sort ([NaN, 1, -1, 2, Inf], 1), [NaN, 1, -1, 2, Inf])
+%!assert (sort ([NaN, 1, -1, 2, Inf], 2), [-1, 1, 2, Inf, NaN])
+%!error (sort ([NaN, 1, -1, 2, Inf], 3))
+%!assert (sort ([NaN, 1, -1, 2, Inf], "ascend"), [-1, 1, 2, Inf, NaN])
+%!assert (sort ([NaN, 1, -1, 2, Inf], 2, "ascend"), [-1, 1, 2, Inf, NaN])
+%!assert (sort ([NaN, 1, -1, 2, Inf], "descend"), [NaN, Inf, 2, 1, -1])
+%!assert (sort ([NaN, 1, -1, 2, Inf], 2, "descend"), [NaN, Inf, 2, 1, -1])
+%!assert (sort ([3, 1, 7, 5; 8, 2, 6, 4]), [3, 1, 6, 4; 8, 2, 7, 5])
+%!assert (sort ([3, 1, 7, 5; 8, 2, 6, 4], 1), [3, 1, 6, 4; 8, 2, 7, 5])
+%!assert (sort ([3, 1, 7, 5; 8, 2, 6, 4], 2), [1, 3, 5, 7; 2, 4, 6, 8])
+%!assert (sort (1), 1)
+
+%!test
+%! [v, i] = sort ([NaN, 1, -1, Inf, 1]);
+%! assert (v, [-1, 1, 1, Inf, NaN])
+%! assert (i, [3, 2, 5, 4, 1])
+
+%% Complex
+%!assert (sort ([NaN, 1i, -1, 2, Inf]), [1i, -1, 2, Inf, NaN])
+%!assert (sort ([NaN, 1i, -1, 2, Inf], 1), [NaN, 1i, -1, 2, Inf])
+%!assert (sort ([NaN, 1i, -1, 2, Inf], 2), [1i, -1, 2, Inf, NaN])
+%!error (sort ([NaN, 1i, -1, 2, Inf], 3))
+%!assert (sort ([NaN, 1i, -1, 2, Inf], "ascend"), [1i, -1, 2, Inf, NaN])
+%!assert (sort ([NaN, 1i, -1, 2, Inf], 2, "ascend"), [1i, -1, 2, Inf, NaN])
+%!assert (sort ([NaN, 1i, -1, 2, Inf], "descend"), [NaN, Inf, 2, -1, 1i])
+%!assert (sort ([NaN, 1i, -1, 2, Inf], 2, "descend"), [NaN, Inf, 2, -1, 1i])
+%!assert (sort ([3, 1i, 7, 5; 8, 2, 6, 4]), [3, 1i, 6, 4; 8, 2, 7, 5])
+%!assert (sort ([3, 1i, 7, 5; 8, 2, 6, 4], 1), [3, 1i, 6, 4; 8, 2, 7, 5])
+%!assert (sort ([3, 1i, 7, 5; 8, 2, 6, 4], 2), [1i, 3, 5, 7; 2, 4, 6, 8])
+%!assert (sort (1i), 1i)
+
+%!test
+%! [v, i] = sort ([NaN, 1i, -1, Inf, 1, 1i]);
+%! assert (v, [1, 1i, 1i, -1, Inf, NaN])
+%! assert (i, [5, 2, 6, 3, 4, 1])
+
+%% Bool
+%!assert (sort ([true, false, true, false]), [false, false, true, true])
+%!assert (sort ([true, false, true, false], 1), [true, false, true, false])
+%!assert (sort ([true, false, true, false], 2), [false, false, true, true])
+%!error (sort ([true, false, true, false], 3))
+%!assert (sort ([true, false, true, false], "ascend"), [false, false, true, true])
+%!assert (sort ([true, false, true, false], 2, "ascend"), [false, false, true, true])
+%!assert (sort ([true, false, true, false], "descend"), [true, true, false, false])
+%!assert (sort ([true, false, true, false], 2, "descend"), [true, true, false, false])
+%!assert (sort (true), true)
+
+%!test
+%! [v, i] = sort ([true, false, true, false]);
+%! assert (v, [false, false, true, true])
+%! assert (i, [2, 4, 1, 3])
+
+%% Sparse Double
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf])), sparse ([-1, 0, 0, 1, 2, Inf, NaN]))
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), 1), sparse ([0, NaN, 1, 0, -1, 2, Inf]))
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), 2), sparse ([-1, 0, 0, 1, 2, Inf, NaN]))
+%!error (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), 3))
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), "ascend"), sparse ([-1, 0, 0, 1, 2, Inf, NaN]))
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), 2, "ascend"), sparse ([-1, 0, 0, 1, 2, Inf, NaN]))
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), "descend"), sparse ([NaN, Inf, 2, 1, 0, 0, -1]))
+%!assert (sort (sparse ([0, NaN, 1, 0, -1, 2, Inf]), 2, "descend"), sparse ([NaN, Inf, 2, 1, 0, 0, -1]))
+
+%!shared a
+%! a = randn (10, 10);
+%! a (a < 0) = 0;
+%!assert (sort (sparse (a)), sparse (sort (a)))
+%!assert (sort (sparse (a), 1), sparse (sort (a, 1)))
+%!assert (sort (sparse (a), 2), sparse (sort (a, 2)))
+%!test
+%! [v, i] = sort (a);
+%! [vs, is] = sort (sparse (a));
+%! assert (vs, sparse (v))
+%! assert (is, i)
+
+%% Sparse Complex
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf])), sparse ([0, 0, 1i, -1, 2, Inf, NaN]))
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), 1), sparse ([0, NaN, 1i, 0, -1, 2, Inf]))
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), 2), sparse ([0, 0, 1i, -1, 2, Inf, NaN]))
+%!error (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), 3))
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), "ascend"), sparse ([0, 0, 1i, -1, 2, Inf, NaN]))
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), 2, "ascend"), sparse ([0, 0, 1i, -1, 2, Inf, NaN]))
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), "descend"), sparse ([NaN, Inf, 2, -1, 1i, 0, 0]))
+%!assert (sort (sparse ([0, NaN, 1i, 0, -1, 2, Inf]), 2, "descend"), sparse ([NaN, Inf, 2, -1, 1i, 0, 0]))
+
+%!shared a
+%! a = randn (10, 10); 
+%! a (a < 0) = 0;
+%! a = 1i * a;
+%!assert (sort (sparse (a)), sparse (sort (a)))
+%!assert (sort (sparse (a), 1), sparse (sort (a, 1)))
+%!assert (sort (sparse (a), 2), sparse (sort (a, 2)))
+%!test
+%! [v, i] = sort (a);
+%! [vs, is] = sort (sparse (a));
+%! assert (vs, sparse (v))
+%! assert (is, i)
+
+%% Sparse Bool
+%!assert (sort (sparse ([true, false, true, false])), sparse ([false, false, true, true]))
+%!assert (sort (sparse([true, false, true, false]), 1), sparse ([true, false, true, false]))
+%!assert (sort (sparse ([true, false, true, false]), 2), sparse ([false, false, true, true]))
+%!error (sort (sparse ([true, false, true, false], 3)))
+%!assert (sort (sparse ([true, false, true, false]), "ascend"), sparse([false, false, true, true]))
+%!assert (sort (sparse ([true, false, true, false]), 2, "ascend"), sparse([false, false, true, true]))
+%!assert (sort (sparse ([true, false, true, false]), "descend"), sparse ([true, true, false, false]))
+%!assert (sort (sparse ([true, false, true, false]), 2, "descend"), sparse([true, true, false, false]))
+
+%!test
+%! [v, i] = sort (sparse([true, false, true, false]));
+%! assert (v, sparse([false, false, true, true]))
+%! assert (i, [2, 4, 1, 3])
+
+%% Cell string array
+%!shared a, b, c
+%! a = {"Alice", "Cecile", "Eric", "Barry", "David"};
+%! b = {"Alice", "Barry", "Cecile", "David", "Eric"};
+%! c = {"Eric", "David", "Cecile", "Barry", "Alice"};
+%!assert (sort (a), b);
+%!assert (sort (a, 1), a)
+%!assert (sort (a, 2), b)
+%!error (sort (a, 3))
+%!assert (sort (a, "ascend"), b)
+%!assert (sort (a, 2, "ascend"), b)
+%!assert (sort (a, "descend"), c)
+%!assert (sort (a, 2, "descend"), c)
+
+%!test
+%! [v, i] = sort (a);
+%! assert (i, [1, 4, 2, 5, 3])
+
+*/
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff -r 3c999b2b5de8 -r 402168152bb9 src/oct-stream.cc
--- a/src/oct-stream.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/oct-stream.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -3220,6 +3220,7 @@
 typedef octave_value (*read_fptr) (octave_stream&, octave_idx_type, octave_idx_type, octave_idx_type, octave_idx_type, bool,
 				   oct_mach_info::float_format ffmt, octave_idx_type&);
 
+NO_INSTANTIATE_ARRAY_SORT (read_fptr);
 INSTANTIATE_ARRAY (read_fptr,);
 template class Array2<read_fptr>;
 
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-base-mat.h
--- a/src/ov-base-mat.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-base-mat.h	Thu Jan 31 18:59:11 2008 +0000
@@ -110,6 +110,12 @@
   MatrixType matrix_type (const MatrixType& _typ) const
     { MatrixType ret = typ; typ = _typ; return ret; }
 
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+    { return octave_value (matrix.sort (dim, mode)); }
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		     sortmode mode = UNDEFINED) const
+    { return octave_value (matrix.sort (sidx, dim, mode)); }
+
   bool is_matrix_type (void) const { return true; }
 
   bool is_numeric_type (void) const { return true; }
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-base-scalar.h
--- a/src/ov-base-scalar.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-base-scalar.h	Thu Jan 31 18:59:11 2008 +0000
@@ -93,6 +93,16 @@
 
   octave_value any (int = 0) const { return (scalar != ST ()); }
 
+  octave_value sort (octave_idx_type, sortmode) const
+    { return octave_value (scalar); }
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type,
+		     sortmode) const
+    { 
+      sidx.resize (dim_vector (1, 1)); 
+      sidx(0) = 0; 
+      return octave_value (scalar); 
+    }
+
   MatrixType matrix_type (void) const { return typ; }
   MatrixType matrix_type (const MatrixType& _typ) const
     { MatrixType ret = typ; typ = _typ; return ret; }
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-base-sparse.h
--- a/src/ov-base-sparse.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-base-sparse.h	Thu Jan 31 18:59:11 2008 +0000
@@ -116,6 +116,12 @@
   octave_value all (int dim = 0) const { return matrix.all (dim); }
   octave_value any (int dim = 0) const { return matrix.any (dim); }
 
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+    { return octave_value (matrix.sort (dim, mode)); }
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		     sortmode mode = UNDEFINED) const
+    { return octave_value (matrix.sort (sidx, dim, mode)); }
+
   MatrixType matrix_type (void) const { return typ; }
   MatrixType matrix_type (const MatrixType& _typ) const
     { MatrixType ret = typ; typ = _typ; return ret; }
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-base.cc
--- a/src/ov-base.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-base.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -878,6 +878,23 @@
   return 0;
 }
 
+octave_value
+octave_base_value::sort (octave_idx_type, sortmode) const
+{
+  gripe_wrong_type_arg ("octave_base_value::sort ()", type_name ());
+
+  return octave_value();
+}
+
+octave_value
+octave_base_value::sort (Array<octave_idx_type> &, 
+			 octave_idx_type, sortmode) const
+{
+  gripe_wrong_type_arg ("octave_base_value::sort ()", type_name ());
+
+  return octave_value();
+}
+
 static void
 gripe_indexed_assignment (const std::string& tn1, const std::string& tn2)
 {
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-base.h
--- a/src/ov-base.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-base.h	Thu Jan 31 18:59:11 2008 +0000
@@ -456,6 +456,12 @@
 
   virtual mxArray *as_mxArray (void) const;
 
+  virtual octave_value sort (octave_idx_type dim = 0, 
+			     sortmode mode = UNDEFINED) const;
+  virtual octave_value sort (Array<octave_idx_type> &sidx, 
+			     octave_idx_type dim = 0,
+			     sortmode mode = UNDEFINED) const;
+
 protected:
 
   // This should only be called for derived types.
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-str-mat.h
--- a/src/ov-str-mat.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-str-mat.h	Thu Jan 31 18:59:11 2008 +0000
@@ -126,6 +126,12 @@
 
   std::string string_value (bool force = false) const;
 
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+  { return octave_value (matrix.sort (dim, mode), true); }
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		     sortmode mode = UNDEFINED) const
+  { return octave_value (matrix.sort (sidx, dim, mode), true); }
+
   bool print_as_scalar (void) const { return (rows () <= 1); }
 
   void print_raw (std::ostream& os, bool pr_as_read_syntax = false) const;
@@ -220,6 +226,13 @@
 			    bool resize_ok = false)
     { return do_index_op_internal (idx, resize_ok, '\''); }
 
+
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+  { return octave_value (matrix.sort (dim, mode), true, '\''); }
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		     sortmode mode = UNDEFINED) const
+  { return octave_value (matrix.sort (sidx, dim, mode), true, '\''); }
+
 private:
 
   DECLARE_OCTAVE_ALLOCATOR
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-streamoff.h
--- a/src/ov-streamoff.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-streamoff.h	Thu Jan 31 18:59:11 2008 +0000
@@ -77,6 +77,18 @@
 
   void print_raw (std::ostream& os, bool pr_as_read_syntax = false) const;
 
+  octave_value sort (octave_idx_type, sortmode) const
+    { 
+      gripe_wrong_type_arg ("octave_base_value::sort ()", type_name ());
+      return octave_value();
+    }
+  octave_value sort (Array<octave_idx_type> &, octave_idx_type,
+		     sortmode) const
+    { 
+      gripe_wrong_type_arg ("octave_base_value::sort ()", type_name ());
+      return octave_value();
+    }
+
 private:
 
   DECLARE_OCTAVE_ALLOCATOR
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov-typeinfo.cc
--- a/src/ov-typeinfo.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov-typeinfo.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -43,26 +43,33 @@
 
 #include <Array.cc>
 
+NO_INSTANTIATE_ARRAY_SORT (octave_value_typeinfo::unary_op_fcn);
 INSTANTIATE_ARRAY (octave_value_typeinfo::unary_op_fcn, );
 template class Array2<octave_value_typeinfo::unary_op_fcn>;
 
+NO_INSTANTIATE_ARRAY_SORT (octave_value_typeinfo::non_const_unary_op_fcn);
 INSTANTIATE_ARRAY (octave_value_typeinfo::non_const_unary_op_fcn, );
 template class Array2<octave_value_typeinfo::non_const_unary_op_fcn>;
 
+NO_INSTANTIATE_ARRAY_SORT (octave_value_typeinfo::binary_op_fcn);
 INSTANTIATE_ARRAY (octave_value_typeinfo::binary_op_fcn, );
 template class Array2<octave_value_typeinfo::binary_op_fcn>;
 template class Array3<octave_value_typeinfo::binary_op_fcn>;
 
+NO_INSTANTIATE_ARRAY_SORT (octave_value_typeinfo::cat_op_fcn);
 INSTANTIATE_ARRAY (octave_value_typeinfo::cat_op_fcn, );
 template class Array2<octave_value_typeinfo::cat_op_fcn>;
 
+NO_INSTANTIATE_ARRAY_SORT (octave_value_typeinfo::assign_op_fcn);
 INSTANTIATE_ARRAY (octave_value_typeinfo::assign_op_fcn, );
 template class Array2<octave_value_typeinfo::assign_op_fcn>;
 template class Array3<octave_value_typeinfo::assign_op_fcn>;
 
+NO_INSTANTIATE_ARRAY_SORT (octave_value_typeinfo::assignany_op_fcn);
 INSTANTIATE_ARRAY (octave_value_typeinfo::assignany_op_fcn, );
 template class Array2<octave_value_typeinfo::assignany_op_fcn>;
 
+NO_INSTANTIATE_ARRAY_SORT (octave_base_value::type_conv_fcn);
 INSTANTIATE_ARRAY (octave_base_value::type_conv_fcn, );
 template class Array2<octave_base_value::type_conv_fcn>;
 
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov.cc
--- a/src/ov.cc	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov.cc	Thu Jan 31 18:59:11 2008 +0000
@@ -588,6 +588,12 @@
   maybe_mutate ();
 }
 
+octave_value::octave_value (const ArrayN<bool>& bnda)
+  : rep (new octave_bool_matrix (bnda))
+{
+  maybe_mutate ();
+}
+
 octave_value::octave_value (char c, char type)
   : rep (type == '"'
 	 ? new octave_char_matrix_dq_str (c)
@@ -680,6 +686,12 @@
   maybe_mutate ();
 }
 
+octave_value::octave_value (const Sparse<bool>& bm, const MatrixType &t)
+  : rep (new octave_sparse_bool_matrix (bm, t))
+{
+  maybe_mutate ();
+}
+
 octave_value::octave_value (const octave_int8& i)
   : rep (new octave_int8_scalar (i))
 {
@@ -851,6 +863,11 @@
 {
 }
 
+octave_value::octave_value (const ArrayN<std::streamoff>& inda)
+  : rep (new octave_streamoff (inda))
+{
+}
+
 octave_value::octave_value (const octave_value_list& l, bool is_csl)
   : rep (is_csl
 	 ? dynamic_cast<octave_base_value *> (new octave_cs_list (l))
diff -r 3c999b2b5de8 -r 402168152bb9 src/ov.h
--- a/src/ov.h	Wed Jan 30 09:11:58 2008 +0000
+++ b/src/ov.h	Thu Jan 31 18:59:11 2008 +0000
@@ -44,6 +44,8 @@
 #include "oct-time.h"
 #include "str-vec.h"
 
+#include "oct-sort.h"
+
 class Cell;
 class streamoff_array;
 class Octave_map;
@@ -183,6 +185,7 @@
   octave_value (bool b);
   octave_value (const boolMatrix& bm, const MatrixType& t = MatrixType());
   octave_value (const boolNDArray& bnda);
+  octave_value (const ArrayN<bool>& bnda);
   octave_value (char c, char type = '"');
   octave_value (const char *s, char type = '"');
   octave_value (const std::string& s, char type = '"');
@@ -200,6 +203,7 @@
   octave_value (const Sparse<Complex>& m, const MatrixType& t = MatrixType ());
   octave_value (const SparseBoolMatrix& bm, 
 		const MatrixType& t = MatrixType ());
+  octave_value (const Sparse<bool>& m, const MatrixType& t = MatrixType ());
   octave_value (const octave_int8& i);
   octave_value (const octave_int16& i);
   octave_value (const octave_int32& i);
@@ -229,6 +233,7 @@
   octave_value (const Octave_map& m);
   octave_value (const Octave_map& m, const std::string& id);
   octave_value (const streamoff_array& off);
+  octave_value (const ArrayN<std::streamoff>& inda);
   octave_value (const octave_value_list& m, bool is_cs_list = false);
   octave_value (octave_value::magic_colon);
 
@@ -860,6 +865,12 @@
 
   mxArray *as_mxArray (void) const { return rep->as_mxArray (); }
 
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = UNDEFINED) const
+    { return rep->sort (dim, mode); } 
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		 sortmode mode = UNDEFINED) const
+    { return rep->sort (sidx, dim, mode); } 
+
 protected:
 
   // The real representation.