[Scipy-svn] r3702 - in trunk/scipy/sparse: . sparsetools tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Dec 24 01:41:41 EST 2007
Author: wnbell
Date: 2007-12-24 00:41:28 -0600 (Mon, 24 Dec 2007)
New Revision: 3702
Modified:
trunk/scipy/sparse/compressed.py
trunk/scipy/sparse/sparsetools/sparsetools.h
trunk/scipy/sparse/tests/test_sparse.py
Log:
updated sparsetools comments
added sort unittest
Modified: trunk/scipy/sparse/compressed.py
===================================================================
--- trunk/scipy/sparse/compressed.py 2007-12-24 05:05:51 UTC (rev 3701)
+++ trunk/scipy/sparse/compressed.py 2007-12-24 06:41:28 UTC (rev 3702)
@@ -192,7 +192,7 @@
elif isspmatrix(other):
if (other.shape != self.shape):
raise ValueError, "inconsistent shapes"
-
+
return self._binopt(other,'_plus_')
elif isdense(other):
# Convert this matrix to a dense matrix and add them
Modified: trunk/scipy/sparse/sparsetools/sparsetools.h
===================================================================
--- trunk/scipy/sparse/sparsetools/sparsetools.h 2007-12-24 05:05:51 UTC (rev 3701)
+++ trunk/scipy/sparse/sparsetools/sparsetools.h 2007-12-24 06:41:28 UTC (rev 3702)
@@ -78,7 +78,7 @@
* Output array Bi must be preallocated
*
* Note:
- * Complexity: Linear.
+ * Complexity: Linear
*
*/
template <class I>
@@ -107,7 +107,7 @@
* I num_blocks - number of blocks
*
* Note:
- * Complexity: Linear.
+ * Complexity: Linear
*
*/
template <class I>
@@ -287,12 +287,13 @@
* I Bj[nnz(B)] - column indices
* T Bx[nnz(B)] - nonzeros
* Output Arguments:
- * vec<I> Cp - row pointer
- * vec<I> Cj - column indices
- * vec<T> Cx - nonzeros
+ * I Cp[n_row+1] - row pointer
+ * I Cj[nnz(C)] - column indices
+ * T Cx[nnz(C)] - nonzeros
*
* Note:
* Output arrays Cp, Cj, and Cx must be preallocated
+ * The value of nnz(C) will be stored in Ap[n_row] after the first pass.
*
* Note:
* Input: A and B column indices *are not* assumed to be in sorted order
@@ -347,7 +348,8 @@
}
/*
- * Pass 2 computes CSR entries for C using the row pointer computed in Pass 1
+ * Pass 2 computes CSR entries for matrix C = A*B using the
+ * row pointer Cp[] computed in Pass 1.
*
*/
template <class I, class T>
@@ -539,80 +541,6 @@
}
(*Cp)[i+1] = Cx->size();
}
-
-
-// //Method that works for unsorted indices
-//
-// Cp->resize(n_brow + 1, 0);
-//
-// const I RC = R*C;
-//
-// std::vector<I> next(n_bcol,-1);
-// std::vector<T> A_row(n_bcol*RC, 0);
-// std::vector<T> B_row(n_bcol*RC, 0);
-//
-// for(I i = 0; i < n_brow; i++){
-// I head = -2;
-// I length = 0;
-//
-// //add a row of A to A_row
-// for(I jj = Ap[i]; jj < Ap[i+1]; jj++){
-// I j = Aj[jj];
-//
-// for(I n = 0; n < RC; n++)
-// A_row[RC*j + n] += Ax[RC*jj + n];
-//
-// if(next[j] == -1){
-// next[j] = head;
-// head = j;
-// length++;
-// }
-// }
-//
-// //add a row of B to B_row
-// for(I jj = Bp[i]; jj < Bp[i+1]; jj++){
-// I j = Bj[jj];
-//
-// for(I n = 0; n < RC; n++)
-// B_row[RC*j + n] += Bx[RC*jj + n];
-//
-// if(next[j] == -1){
-// next[j] = head;
-// head = j;
-// length++;
-// }
-// }
-//
-//
-// for(I jj = 0; jj < length; jj++){
-// bool nonzero_block = false;
-// for(I n = 0; n < RC; n++){
-// T result = op(A_row[RC*head + n],B_row[RC*head + n]);
-// A_row[RC*head + n] = result;
-// if(result != 0)
-// nonzero_block = true;
-// }
-//
-//
-// if(nonzero_block){
-// Cj->push_back(head);
-// for(I n = 0; n < RC; n++){
-// Cx->push_back(A_row[RC*head + n]);
-// }
-// }
-//
-// for(I n = 0; n < RC; n++){
-// A_row[RC*head + n] = 0;
-// B_row[RC*head + n] = 0;
-// }
-//
-// I temp = head;
-// head = next[head];
-// next[temp] = -1;
-// }
-//
-// (*Cp)[i+1] = Cj->size();
-// }
}
/* element-wise binary operations*/
@@ -659,7 +587,7 @@
/*
* Compute C = A (bin_op) B for CSR matrices A,B
*
- * (bin_op) - binary operator to apply elementwise
+ * bin_op(x,y) - binary operator to apply elementwise
*
*
* Input Arguments:
@@ -672,16 +600,18 @@
* I Bj[nnz(B)] - column indices
* T Bx[nnz(B)] - nonzeros
* Output Arguments:
- * vec<I> Cp - row pointer
- * vec<I> Cj - column indices
- * vec<T> Cx - nonzeros
+ * I Cp[n_row+1] - row pointer
+ * I Cj[nnz(C)] - column indices
+ * T Cx[nnz(C)] - nonzeros
*
* Note:
- * Output arrays Cp, Cj, and Cx will be allocated within in the method
+ * Output arrays Cp, Cj, and Cx must be preallocated
+ * If nnz(C) is not known a priori, a conservative bound is:
+ * nnz(C) <= nnz(A) + nnz(B)
*
* Note:
- * Input: A and B column indices *are not* assumed to be in sorted order
- * Output: C column indices *are not* assumed to be in sorted order
+ * Input: A and B column indices are assumed to be in sorted order
+ * Output: C column indices are assumed to be in sorted order
* Cx will not contain any zero entries
*
*/
@@ -699,12 +629,10 @@
T Cx[],
const bin_op& op)
{
- //Method that works for sorted indices
- assert( csr_has_sorted_indices(n_row,Ap,Aj) );
- assert( csr_has_sorted_indices(n_row,Bp,Bj) );
+ //Method that works for sorted indices
+ // assert( csr_has_sorted_indices(n_row,Ap,Aj) );
+ // assert( csr_has_sorted_indices(n_row,Bp,Bj) );
- //Cp->resize(n_row + 1, 0);
- //(*Cp)[0] = 0;
Cp[0] = 0;
I nnz = 0;
@@ -716,7 +644,7 @@
I A_j = Aj[A_pos];
I B_j = Bj[B_pos];
-
+
//while not finished with either row
while(A_pos < A_end && B_pos < B_end){
if(A_j == B_j){
@@ -769,68 +697,6 @@
}
Cp[i+1] = nnz;
}
-
-
-// //Method that works for unsorted indices
-// Cp->resize(n_row + 1, 0);
-//
-// std::vector<I> next(n_col,-1);
-// std::vector<T> A_row(n_col, 0);
-// std::vector<T> B_row(n_col, 0);
-//
-// for(I i = 0; i < n_row; i++){
-// I head = -2;
-// I length = 0;
-//
-// //add a row of A to A_row
-// I i_start = Ap[i];
-// I i_end = Ap[i+1];
-// for(I jj = i_start; jj < i_end; jj++){
-// I j = Aj[jj];
-//
-// A_row[j] += Ax[jj];
-//
-// if(next[j] == -1){
-// next[j] = head;
-// head = j;
-// length++;
-// }
-// }
-//
-// //add a row of B to B_row
-// i_start = Bp[i];
-// i_end = Bp[i+1];
-// for(I jj = i_start; jj < i_end; jj++){
-// I j = Bj[jj];
-//
-// B_row[j] += Bx[jj];
-//
-// if(next[j] == -1){
-// next[j] = head;
-// head = j;
-// length++;
-// }
-// }
-//
-//
-// for(I jj = 0; jj < length; jj++){
-// T result = op(A_row[head],B_row[head]);
-//
-// if(result != 0){
-// Cj->push_back(head);
-// Cx->push_back(result);
-// }
-//
-// I temp = head;
-// head = next[head];
-//
-// next[temp] = -1;
-// A_row[temp] = 0;
-// B_row[temp] = 0;
-// }
-//
-// (*Cp)[i+1] = Cj->size();
-// }
}
/* element-wise binary operations*/
Modified: trunk/scipy/sparse/tests/test_sparse.py
===================================================================
--- trunk/scipy/sparse/tests/test_sparse.py 2007-12-24 05:05:51 UTC (rev 3701)
+++ trunk/scipy/sparse/tests/test_sparse.py 2007-12-24 06:41:28 UTC (rev 3702)
@@ -12,6 +12,13 @@
from scipy.linsolve import splu
restore_path()
+def random_sparse(m,n,nnz_per_row):
+ rows = numpy.arange(m).repeat(nnz_per_row)
+ cols = numpy.random.random_integers(low=0,high=n-1,size=nnz_per_row*m)
+ vals = numpy.random.random_sample(m*nnz_per_row)
+ return coo_matrix((vals,(rows,cols)),(m,n)).tocsr()
+
+
#TODO move this to a matrix gallery and add unittests
def poisson2d(N,dtype='d',format=None):
"""
@@ -39,14 +46,12 @@
class TestSparseTools(NumpyTestCase):
"""Simple benchmarks for sparse matrix module"""
- def test_arithmetic(self,level=4):
+ def test_arithmetic(self,level=5):
matrices = []
#matrices.append( ('A','Identity', spidentity(500**2,format='csr')) )
matrices.append( ('A','Poisson5pt', poisson2d(500,format='csr')) )
matrices.append( ('B','Poisson5pt^2', poisson2d(500,format='csr')**2) )
- #matrices = [ (a,b,c.astype('int8')) for (a,b,c) in matrices ]
-
print
print ' Sparse Matrix Arithmetic'
print '===================================================================='
@@ -87,9 +92,43 @@
operation = (X + '.' + op + '(' + Y + ')').center(17)
print fmt % (format,operation,msec_per_it)
+
- def test_matvec(self,level=5):
+ def bench_sort(self,level=4):
+ """sort CSR column indices"""
matrices = []
+ matrices.append( ('Rand10', 1e4, 10) )
+ matrices.append( ('Rand25', 1e4, 25) )
+ matrices.append( ('Rand50', 1e4, 50) )
+ matrices.append( ('Rand100', 1e4, 100) )
+ matrices.append( ('Rand200', 1e4, 200) )
+
+ print
+ print ' Sparse Matrix Index Sorting'
+ print '====================================================================='
+ print ' type | name | shape | nnz | time (msec) '
+ print '---------------------------------------------------------------------'
+ fmt = ' %3s | %12s | %20s | %8d | %6.2f '
+
+ for name,N,K in matrices:
+ N = int(N)
+ A = random_sparse(N,N,K)
+
+ start = time.clock()
+ iter = 0
+ while iter < 5 and time.clock() - start < 1:
+ A.sort_indices(check_first =False)
+ iter += 1
+ end = time.clock()
+
+ name = name.center(12)
+ shape = ("%s" % (A.shape,)).center(20)
+
+ print fmt % (A.format,name,shape,A.nnz,1e3*(end-start)/float(iter) )
+
+
+ def bench_matvec(self,level=5):
+ matrices = []
matrices.append(('Identity', spidentity(10**4,format='dia')))
matrices.append(('Identity', spidentity(10**4,format='csr')))
matrices.append(('Poisson5pt', poisson2d(300,format='dia')))
More information about the Scipy-svn
mailing list