Mercurial > octave-antonio
changeset 7433:402168152bb9
[project @ 2008-01-31 18:59:09 by dbateman]
| author | dbateman |
|---|---|
| date | Thu, 31 Jan 2008 18:59:11 +0000 |
| parents | 3c999b2b5de8 |
| children | 600a3af7963e |
| files | liboctave/Array-C.cc liboctave/Array-b.cc liboctave/Array-ch.cc liboctave/Array-d.cc liboctave/Array-i.cc liboctave/Array-idx-vec.cc liboctave/Array-s.cc liboctave/Array-str.cc liboctave/Array.cc liboctave/Array.h liboctave/Array2.h liboctave/Array3.h liboctave/ArrayN.h liboctave/CSparse.cc liboctave/ChangeLog liboctave/Sparse-C.cc liboctave/Sparse-b.cc liboctave/Sparse-d.cc liboctave/Sparse.cc liboctave/Sparse.h liboctave/dSparse.cc liboctave/oct-sort.cc liboctave/oct-sort.h src/ChangeLog src/DLD-FUNCTIONS/sort.cc src/Makefile.in src/TEMPLATE-INST/Array-os.cc src/TEMPLATE-INST/Array-tc.cc src/data.cc src/oct-stream.cc src/ov-base-mat.h src/ov-base-scalar.h src/ov-base-sparse.h src/ov-base.cc src/ov-base.h src/ov-str-mat.h src/ov-streamoff.h src/ov-typeinfo.cc src/ov.cc src/ov.h |
| diffstat | 40 files changed, 1616 insertions(+), 1521 deletions(-) [+] |
line wrap: on
line diff
--- 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);
--- 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);
--- 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);
--- 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);
--- 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);
--- 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); /*
--- 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);
--- 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); /*
--- 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>
--- 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 /*
--- 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
--- 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
--- 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>
--- 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
--- 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,
--- 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);
--- 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);
--- 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);
--- 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.
--- 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 /*
--- 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
--- 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; }
--- 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 /*
--- 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.
--- 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: *** -*/
--- 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 \
--- 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); /*
--- 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);
--- 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++ ***
--- 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>;
--- 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; }
--- 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; }
--- 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; }
--- 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) {
--- 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.
--- 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
--- 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
--- 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>;
--- 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))
--- 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.