From octave-maintainers-request at bevo dot che dot wisc dot edu Wed Jan 14 17:07:40 2004
Subject: Re: benchmarks - sort
From: David Bateman <David dot Bateman at motorola dot com>
To: Schloegl Alois <alois dot schloegl at tugraz dot at>
Cc: Paul Kienzle <pkienzle at users dot sourceforge dot net>,    octave-maintainers@bevo.che.wisc.edu
Date: Wed, 14 Jan 2004 12:12:51 -0600


--i7F3eY7HS/tUJxUd
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline


Well, for the hell of it I programmed a simple sort oct-file to test 
what would be faster (IEEE754 integer compare or floating point compare
with checking for Inf and NaN), you can check the attached oct-file.
One complication of IEEE754 is the symmetry of negative and postive numbers.
This means that positive and negative tests need to be treated seperately, as
opposed to a 2's complement respresentation 

The code at the moment is only for a row vector and is attached the example
I used was

a=randn(1,100000); a(1) = NaN; a(2) = -Inf; a(3) = Inf; tic; mysort(a); t = toc; tic; mysort(a,1); t1= toc; t, t1

And I got the result

t = 0.27054
t1 = 0.23212

showing that IEEE754 integer compares are about 20% faster....  Well
IEEE754 is faster, but it is certainly not a factor of 10 faster...
Seems like we need to look elsewhere for this order of magnitude
increase...

Regards
D.


According to Schloegl Alois <alois dot schloegl at tugraz dot at> (on 01/06/04):
> Zitat von Paul Kienzle <pkienzle at users dot sourceforge dot net>:
> 
> >
> > Paul Thomas pointed me to these benchmarks.  Anybody want to do
> > something about them?
> >
> > 	http://www.sciviews.org/other/benchmark.htm
> >
> 
> ...
> 
> >
> > I.C Matlab 0.89  - Octave 7.77
> >      b = sort(a);
> >
> > Octave's sort is surprisingly slow.  3x worse than any other package
> > mentioned.  Anyone know a fast stable sort algorithm?
> >
> 
> 
> One reason, for the bad performance of the sort algorithm might be the exception
> handling of NaNs. See also
> http://www.octave.org/octave-lists/archive/bug-octave.2001/msg00047.html
> 
> Sorting on the binary level might be helpful, because the bit patterns of the
> IEEE754 numbers provide the correct sorting order. This is
> -inf < -1 < 0 < 1 < inf < NaN
> 
> Just my few thougths.
> 
> 
> Alois
> 

-- 
David Bateman                                David dot Bateman at motorola dot com
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

--i7F3eY7HS/tUJxUd
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="mysort.cc"

#include <octave/oct.h>
#include <octave/lo-mappers.h>
#include <octave/config.h>
#include <octave/quit.h>

static Array<int>
create_index_array (int n)
{
  Array<int> l (n+2);

  l (0) = 1;

  for (int i = 1; i < n - 1; i++)
    l (i) = -(i+2);

  l (n-1) = 0;
  l (n) = 0;
  l (n+1) = 2;

  return l;
}

#define SORT_INIT_PHASE(n) \
  int s = 0; \
  int t = n + 1; \
  int p = l (s); \
  int q = l (t); \
  if (q == 0) \
     break

#define SORT_COMMON_CODE \
  p = -p; \
  q = -q; \
  if (q == 0) \
    { \
      l (s) = (l (s) < 0) \
	? ((p < 0) ? p : -p) \
	  : ((p >= 0) ? p : -p); \
      l (t) = 0; \
      break; \
    } \

#define SORT_REORDER_PHASE_ONE \
  l (s) = (l (s) < 0) \
    ? ((q < 0) ? q : -q) \
      : ((q >= 0) ? q : -q); \
  s = q; \
  q = l (q); \
  if (q <= 0) \
    { \
      l (s) = p; \
      s = t; \
      do \
	{ \
	  t = p; \
	  p = l (p); \
	} \
      while (p > 0); \
      SORT_COMMON_CODE; \
    } \

#define SORT_REORDER_PHASE_TWO \
  l (s) = (l (s) < 0) \
    ? ((p < 0) ? p : -p) \
      : ((p >= 0) ? p : -p); \
  s = p; \
  p = l (p); \
  if (p <= 0) \
    { \
      l (s) = q; \
      s = t; \
      do \
	{ \
	  t = q; \
	  q = l (q); \
	} \
      while (q > 0); \
      SORT_COMMON_CODE; \
    }

#define DO_SORT(n, condition) \
  while (1) \
    { \
      SORT_INIT_PHASE(n); \
      while (1) \
	{ \
          OCTAVE_QUIT; \
	  if (condition) \
	    { \
	      SORT_REORDER_PHASE_ONE; \
	    } \
	  else \
	    { \
	      SORT_REORDER_PHASE_TWO; \
	    } \
	} \
    }

#define DO_SORT_IEEE754(n) \
  while (1) \
    { \
      SORT_INIT_PHASE(n); \
      while (1) \
	{ \
          OCTAVE_QUIT; \
	  unsigned long long v1 = vec_int[p-1]; \
	  unsigned long long v2 = vec_int[q-1]; \
	  bool s1 = (mask & v1) != 0; \
	  bool s2 = (mask & v2) != 0; \
	  if ((!s1 & s2) || ( s1 & s2 & (v1 < v2)) || \
		( !s1 & !s2 & (v1 > v2))) \
	    { \
	      SORT_REORDER_PHASE_ONE; \
	    } \
	  else \
	    { \
	      SORT_REORDER_PHASE_TWO; \
	    } \
	} \
    }

#define VECTOR_CREATE_RETURN_VALUES(vs, v) \
  int k = l (0); \
  idx (0) = k; \
  vs (0) = v (k-1); \
  for (int i = 1; i < n; i++) \
    { \
      OCTAVE_QUIT; \
      k = l (static_cast<int> (idx (i-1))); \
      idx (i) = k; \
      vs (i) = v (k-1); \
    }

DEFUN_DLD (mysort, args, , "") {
  octave_value_list retval;

  RowVector v = args(0).row_vector_value();
  int n = v.capacity();
  RowVector vs (n);
  RowVector idx (n);
  
  // Simple bubble sort for now
  if (n == 1)
    {
      retval(0) = v;
      retval(1) = RowVector (n, 1.0);
      return retval;
    }
  else if (args.length() > 1)
    {
      // Use integer compare al la IEEE754
      double * vec_dbl = v.fortran_vec();
      unsigned long long *vec_int = (unsigned long long *)vec_dbl;
      unsigned long long mask = (unsigned long long)1 << ((sizeof(unsigned long long)<<3) -1);

      Array<int> l = create_index_array (n);

      DO_SORT_IEEE754 (n);

      VECTOR_CREATE_RETURN_VALUES (vs, v);

      retval(1) = idx;
      retval(0) = vs;
    }
  else
    {
      Array<int> l = create_index_array (n);

      DO_SORT (n, (xisnan (v (p-1)) || v (p-1) > v (q-1)));

      VECTOR_CREATE_RETURN_VALUES (vs, v);

      retval(1) = idx;
      retval(0) = vs;
    }

  return retval;
}

--i7F3eY7HS/tUJxUd--