[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: type of sort(sparse(...))
|
From: |
David Bateman |
|
Subject: |
Re: type of sort(sparse(...)) |
|
Date: |
Mon, 03 Sep 2007 00:39:11 +0200 |
|
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
Søren Hauberg wrote:
> David Bateman skrev:
>> Søren Hauberg wrote:
>>> David Bateman skrev:
>>>> Could someone who has a recent copy of matlab tell me what
>>>>
>>>> a = sort (sprand(10,10,0.2));
>>>> whos a
>>>>
>>>> gives? In Octave this converts the sparse matrix to a full one, whereas
>>>> I suspect in matlab it might not.
>>> Matlab Version 7.4.0.129 (R2007a):
>>>
>>>>> a = sort (sprand(10,10,0.2));
>>>>> whos a
>>> Name Size Bytes Class Attributes
>>>
>>> a 10x10 248 double sparse
>>>
>>
>> Now that I think about it, what does
>>
>> [a, i] = sort (sprand (10,10,0.2))
>> whos i
>>
>> give.. For good measure note I ask what the values of a and i are.. My
>> worry here is that as a is sparse above, I don't see how i can be and so
>> it makes no sense to do an indexed sort of a sparse matrix..
>>> [a, i] = sort (sprand (10,10,0.2));
>>> whos
> Name Size Bytes Class Attributes
>
> a 10x10 284 double sparse
> i 10x10 800 double
>
> Søren
>
Ok, then that is an ugly mess and yes indexed sparse sort is not a good
idea for memory reasons as the index matrix occupies so much space.
Taking advantage of the sparsity does however speed up the sort even for
indexed cases. See the attached patch.
D.
*** ./src/DLD-FUNCTIONS/sort.cc.orig17 2007-09-02 21:10:39.712088645 +0200
--- ./src/DLD-FUNCTIONS/sort.cc 2007-09-03 00:09:08.828494641 +0200
***************
*** 236,241 ****
--- 236,423 ----
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.
***************
*** 602,607 ****
--- 784,797 ----
#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)
***************
*** 611,616 ****
--- 801,836 ----
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)
***************
*** 629,640 ****
--- 849,870 ----
&& (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>;
***************
*** 821,851 ****
if (arg.is_real_type ())
{
! 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 ())
{
! ComplexNDArray cm = arg.complex_array_value ();
! // Don't have unindexed version as no ">" operator
! if (! error_state)
! retval = mx_sort_indexed (cm, dim, smode);
}
else if (arg.is_string ())
{
--- 1051,1113 ----
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
! {
! 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 ())
{
*** ./src/ov.cc.orig17 2007-05-21 22:16:48.000000000 +0200
--- ./src/ov.cc 2007-09-02 21:09:01.468093851 +0200
***************
*** 534,545 ****
--- 534,557 ----
maybe_mutate ();
}
+ octave_value::octave_value (const Sparse<double>& m, const MatrixType &t)
+ : rep (new octave_sparse_matrix (m, t))
+ {
+ maybe_mutate ();
+ }
+
octave_value::octave_value (const SparseComplexMatrix& m, const MatrixType &t)
: rep (new octave_sparse_complex_matrix (m, t))
{
maybe_mutate ();
}
+ octave_value::octave_value (const Sparse<Complex>& m, const MatrixType &t)
+ : rep (new octave_sparse_complex_matrix (m, t))
+ {
+ maybe_mutate ();
+ }
+
octave_value::octave_value (const SparseBoolMatrix& bm, const MatrixType &t)
: rep (new octave_sparse_bool_matrix (bm, t))
{
*** ./src/ov-cx-sparse.h.orig17 2007-09-02 21:28:29.009611479 +0200
--- ./src/ov-cx-sparse.h 2007-09-02 21:51:34.476026554 +0200
***************
*** 73,78 ****
--- 73,89 ----
octave_sparse_complex_matrix (const MSparse<Complex>& m)
: octave_base_sparse<SparseComplexMatrix> (m) { }
+ octave_sparse_complex_matrix (const MSparse<Complex>& m,
+ const MatrixType &t)
+ : octave_base_sparse<SparseComplexMatrix> (m, t) { }
+
+ octave_sparse_complex_matrix (const Sparse<Complex>& m,
+ const MatrixType &t)
+ : octave_base_sparse<SparseComplexMatrix> (SparseComplexMatrix (m), t) { }
+
+ octave_sparse_complex_matrix (const Sparse<Complex>& m)
+ : octave_base_sparse<SparseComplexMatrix> (SparseComplexMatrix (m)) { }
+
octave_sparse_complex_matrix (const octave_sparse_complex_matrix& cm)
: octave_base_sparse<SparseComplexMatrix> (cm) { }
*** ./src/ov-re-sparse.h.orig17 2007-09-02 21:28:06.208773106 +0200
--- ./src/ov-re-sparse.h 2007-09-02 21:50:39.809811617 +0200
***************
*** 73,78 ****
--- 73,87 ----
octave_sparse_matrix (const MSparse<double>& m)
: octave_base_sparse<SparseMatrix> (m) { }
+ octave_sparse_matrix (const MSparse<double>& m, const MatrixType& t)
+ : octave_base_sparse<SparseMatrix> (m, t) { }
+
+ octave_sparse_matrix (const Sparse<double>& m)
+ : octave_base_sparse<SparseMatrix> (SparseMatrix (m)) { }
+
+ octave_sparse_matrix (const Sparse<double>& m, const MatrixType& t)
+ : octave_base_sparse<SparseMatrix> (SparseMatrix (m), t) { }
+
octave_sparse_matrix (const octave_sparse_matrix& m)
: octave_base_sparse<SparseMatrix> (m) { }
*** ./src/ov.h.orig17 2007-06-13 07:54:14.000000000 +0200
--- ./src/ov.h 2007-09-02 21:09:46.540797546 +0200
***************
*** 192,199 ****
--- 192,201 ----
octave_value (const ArrayN<char>& chnda, bool is_string = false,
char type = '"');
octave_value (const SparseMatrix& m, const MatrixType& t = MatrixType ());
+ octave_value (const Sparse<double>& m, const MatrixType& t = MatrixType ());
octave_value (const SparseComplexMatrix& m,
const MatrixType& t = MatrixType ());
+ octave_value (const Sparse<Complex>& m, const MatrixType& t = MatrixType
());
octave_value (const SparseBoolMatrix& bm,
const MatrixType& t = MatrixType ());
octave_value (const octave_int8& i);
2007-09-01 David Bateman <address@hidden>
* DLD-FUNCTIONS/sort.cc (mx_sort_sparse, mx_sort_sparse_indexed):
New template classes for sparse sort functions.
(Fsort): Use them.
* ov.h (octave_value (const Sparse<double>&, const MatrixType&),
octave_value (const Sparse<Complex>&, const MatrixType&)): New
constructors.
* ov.cc (octave_value::octave_value (const Sparse<double>&,
const MatrixType&), octave_value::octave_value (const
Sparse<Complex>&, const MatrixType&)): Define them.
* ov-re-sparse.h (octave_sparse_matrix (const MSparse<double>&,
const MatrixType&), octave_sparse_matrix (const Sparse<double>&),
octave_sparse_matrix (const Sparse<double>&, const MatrixType&)):
New constructors.
* ov-cx-sparse.h (octave_sparse_complex_matrix (const MSparse<double>&,
const MatrixType&), octave_sparse_complex_matrix (const
Sparse<double>&), octave_sparse_complex_matrix (const
Sparse<double>&, const MatrixType&)): ditto.