[Scipy-svn] r3695 - in trunk/scipy/sparse: . sparsetools tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Dec 20 21:04:29 EST 2007
Author: wnbell
Date: 2007-12-20 20:04:21 -0600 (Thu, 20 Dec 2007)
New Revision: 3695
Modified:
trunk/scipy/sparse/__init__.py
trunk/scipy/sparse/base.py
trunk/scipy/sparse/block.py
trunk/scipy/sparse/bsr.py
trunk/scipy/sparse/coo.py
trunk/scipy/sparse/csc.py
trunk/scipy/sparse/csr.py
trunk/scipy/sparse/sparsetools/fixed_size.h
trunk/scipy/sparse/sparsetools/sparsetools.h
trunk/scipy/sparse/sparsetools/sparsetools.py
trunk/scipy/sparse/sparsetools/sparsetools_wrap.cxx
trunk/scipy/sparse/tests/test_base.py
Log:
enabled bsr_matrix
improved bsr_matrix.matvec()
Modified: trunk/scipy/sparse/__init__.py
===================================================================
--- trunk/scipy/sparse/__init__.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/__init__.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -9,8 +9,7 @@
from dok import *
from coo import *
from dia import *
-#from bsr import *
-#from bsc import *
+from bsr import *
from construct import *
from spfuncs import *
Modified: trunk/scipy/sparse/base.py
===================================================================
--- trunk/scipy/sparse/base.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/base.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -407,12 +407,9 @@
def todia(self):
return self.tocoo().todia()
-# def tobsr(self,blocksize=None):
-# return self.tocoo().tobsr(blocksize=blocksize)
-#
-# def tobsc(self,blocksize=None):
-# return self.tocoo().tobsc(blocksize=blocksize)
-
+ def tobsr(self,blocksize=None):
+ return self.tocsr().tobsr(blocksize=blocksize)
+
def copy(self):
return self.__class__(self,copy=True)
Modified: trunk/scipy/sparse/block.py
===================================================================
--- trunk/scipy/sparse/block.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/block.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -61,6 +61,21 @@
arg1 = self.__class__( coo_matrix(arg1), blocksize=blocksize )
self._set_self( arg1 )
+ if shape is not None:
+ self.shape = shape # spmatrix will check for errors
+ else:
+ if self.shape is None:
+ # shape not already set, try to infer dimensions
+ try:
+ major_dim = len(self.indptr) - 1
+ minor_dim = self.indices.max() + 1
+ except:
+ raise ValueError,'unable to infer matrix dimensions'
+ else:
+ M,N = self._swap((major_dim,minor_dim))
+ R,C = self.blocksize
+ self.shape = (M*R,N*C)
+
if self.shape is None:
if shape is None:
#infer shape here
@@ -141,8 +156,8 @@
blocksize = property(fget=_get_blocksize)
def getnnz(self):
- X,Y = self.blocksize
- return self.indptr[-1] * X * Y
+ R,C = self.blocksize
+ return self.indptr[-1] * R * C
nnz = property(fget=getnnz)
def __repr__(self):
Modified: trunk/scipy/sparse/bsr.py
===================================================================
--- trunk/scipy/sparse/bsr.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/bsr.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -39,7 +39,6 @@
if other.shape != (N,) and other.shape != (N,1):
raise ValueError, "dimension mismatch"
-
#output array
if output is None:
@@ -76,19 +75,15 @@
-
-
-
-
- def transpose(self,copy=False):
- M,N = self.shape
-
- data = self.data.swapaxes(1,2)
- indices = self.indices
- indptr = self.indptr
-
- from bsc import bsc_matrix
- return bsc_matrix( (data,indices,indptr), shape=(N,M), copy=copy)
+# def transpose(self,copy=False):
+# M,N = self.shape
+#
+# data = self.data.swapaxes(1,2)
+# indices = self.indices
+# indptr = self.indptr
+#
+# from bsc import bsc_matrix
+# return bsc_matrix( (data,indices,indptr), shape=(N,M), copy=copy)
def tocoo(self,copy=True):
"""Convert this matrix to COOrdinate format.
@@ -118,17 +113,11 @@
return coo_matrix( (data,(row,col)), shape=self.shape )
- def tobsc(self,blocksize=None):
- if blocksize is None:
- blocksize = self.blocksize
- elif blocksize != self.blocksize:
- return self.tocoo(copy=False).tobsc(blocksize=blocksize)
-
- #maintian blocksize
+ def transpose(self):
X,Y = self.blocksize
M,N = self.shape
- #use CSR->CSC to determine a permutation for BSR<->BSC
+ #use CSR.T to determine a permutation for BSR.T
from csr import csr_matrix
data = arange(len(self.indices), dtype=self.indices.dtype)
proxy = csr_matrix((data,self.indices,self.indptr),shape=(M/X,N/Y))
@@ -138,8 +127,7 @@
indices = proxy.indices
indptr = proxy.indptr
- from bsc import bsc_matrix
- return bsc_matrix( (data,indices,indptr), shape=self.shape )
+ return bsr_matrix( (data,indices,indptr), shape=(N,M) )
def tobsr(self,blocksize=None,copy=False):
Modified: trunk/scipy/sparse/coo.py
===================================================================
--- trunk/scipy/sparse/coo.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/coo.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -293,58 +293,61 @@
return dok
-
# def tobsc(self,blocksize=None):
# if blocksize in [None, (1,1)]:
# return self.tocsc().tobsc(blocksize)
# else:
# return self.transpose().tobsr().transpose()
-#
-# def tobsr(self,blocksize=None):
-# if blocksize in [None, (1,1)]:
-# return self.tocsr().tobsr(blocksize)
-#
-# M,N = self.shape
-# X,Y = blocksize
-#
-# if (M % X) != 0 or (N % Y) != 0:
-# raise ValueError, 'shape must be multiple of blocksize'
-#
-# i_block,i_sub = divmod(self.row, X)
-# j_block,j_sub = divmod(self.col, Y)
-#
-# perm = lexsort( keys=[j_block,i_block] )
-#
-# i_block = i_block[perm]
-# j_block = j_block[perm]
-#
-# mask = (i_block[1:] != i_block[:-1]) + (j_block[1:] != j_block[:-1])
-# mask = concatenate((array([True]),mask))
-#
-# #map self.data[n] -> data[map[n],i_sub[n],j_sub[n]]
-# map = cumsum(mask)
-# num_blocks = map[-1]
-# map -= 1
-#
-# iperm = empty_like(perm) #inverse permutation
-# iperm[perm] = arange(len(perm))
-#
-# data = zeros( (num_blocks,X,Y), dtype=self.dtype )
-# data[map[iperm],i_sub,j_sub] = self.data
-#
-# row = i_block[mask]
-# col = j_block[mask]
-#
-# # now row,col,data form BOO format
-#
-# temp = cumsum(bincount(row))
-# indptr = zeros( M/X + 1, dtype=intc )
-# indptr[1:len(temp)+1] = temp
-# indptr[len(temp)+1:] = temp[-1]
-#
-# from bsr import bsr_matrix
-# return bsr_matrix((data,col,indptr),shape=self.shape)
+ def tobsr(self,blocksize=None):
+ from bsr import bsr_matrix
+
+ if self.nnz == 0:
+ return bsr_matrix(self.shape,blocksize=blocksize,dtype=self.dtype)
+
+ if blocksize in [None, (1,1)]:
+ return self.tocsr().tobsr(blocksize)
+
+ M,N = self.shape
+ X,Y = blocksize
+
+ if (M % X) != 0 or (N % Y) != 0:
+ raise ValueError, 'shape must be multiple of blocksize'
+
+ i_block,i_sub = divmod(self.row, X)
+ j_block,j_sub = divmod(self.col, Y)
+
+ perm = lexsort( keys=[j_block,i_block] )
+
+ i_block = i_block[perm]
+ j_block = j_block[perm]
+
+ mask = (i_block[1:] != i_block[:-1]) + (j_block[1:] != j_block[:-1])
+ mask = concatenate((array([True]),mask))
+
+ #map self.data[n] -> data[map[n],i_sub[n],j_sub[n]]
+ map = cumsum(mask)
+ num_blocks = map[-1]
+ map -= 1
+
+ iperm = empty_like(perm) #inverse permutation
+ iperm[perm] = arange(len(perm))
+
+ data = zeros( (num_blocks,X,Y), dtype=self.dtype )
+ data[map[iperm],i_sub,j_sub] = self.data
+
+ row = i_block[mask]
+ col = j_block[mask]
+
+ # now row,col,data form BOO format
+
+ temp = cumsum(bincount(row))
+ indptr = zeros( M/X + 1, dtype=intc )
+ indptr[1:len(temp)+1] = temp
+ indptr[len(temp)+1:] = temp[-1]
+
+ return bsr_matrix((data,col,indptr),shape=self.shape)
+
# needed by _data_matrix
def _with_data(self,data,copy=True):
"""Returns a matrix with the same sparsity structure as self,
Modified: trunk/scipy/sparse/csc.py
===================================================================
--- trunk/scipy/sparse/csc.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/csc.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -225,17 +225,17 @@
# return bsc_matrix( arg1, shape=self.shape, copy=copy )
# else:
# #TODO make this more efficient
-# return self.tocoo(copy=False).tobsr(blocksize=blocksize)
+# return self.tocoo(copy=False).tobsc(blocksize=blocksize)
#
-# def tobsr(self, blocksize=None):
-# if blocksize in [None, (1,1)]:
-# from bsr import bsr_matrix
-# csr = self.tocsr()
-# arg1 = (csr.data.reshape(-1,1,1),csr.indices,csr.indptr)
-# return bsr_matrix( arg1, shape=self.shape )
-# else:
-# #TODO make this more efficient
-# return self.tocoo(copy=False).tobsr(blocksize=blocksize)
+ def tobsr(self, blocksize=None):
+ if blocksize in [None, (1,1)]:
+ from bsr import bsr_matrix
+ csr = self.tocsr()
+ arg1 = (csr.data.reshape(-1,1,1),csr.indices,csr.indptr)
+ return bsr_matrix( arg1, shape=self.shape )
+ else:
+ #TODO make this more efficient
+ return self.tocoo(copy=False).tobsr(blocksize=blocksize)
def get_submatrix( self, slice0, slice1 ):
"""Return a submatrix of this matrix (new matrix is created).
Modified: trunk/scipy/sparse/csr.py
===================================================================
--- trunk/scipy/sparse/csr.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/csr.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -224,15 +224,15 @@
from csc import csc_matrix
return csc_matrix((data, indices, indptr), self.shape)
-# def tobsr(self,blocksize=None,copy=True):
-# if blocksize in [None, (1,1)]:
-# from bsr import bsr_matrix
-# arg1 = (self.data.reshape(-1,1,1),self.indices,self.indptr)
-# return bsr_matrix( arg1, shape=self.shape, copy=copy )
-# else:
-# #TODO make this more efficient
-# return self.tocoo(copy=False).tobsr(blocksize=blocksize)
-#
+ def tobsr(self,blocksize=None,copy=True):
+ if blocksize in [None, (1,1)]:
+ from bsr import bsr_matrix
+ arg1 = (self.data.reshape(-1,1,1),self.indices,self.indptr)
+ return bsr_matrix( arg1, shape=self.shape, copy=copy )
+ else:
+ #TODO make this more efficient
+ return self.tocoo(copy=False).tobsr(blocksize=blocksize)
+
# def tobsc(self,blocksize=None):
# if blocksize in [None, (1,1)]:
# from bsc import bsc_matrix
Modified: trunk/scipy/sparse/sparsetools/fixed_size.h
===================================================================
--- trunk/scipy/sparse/sparsetools/fixed_size.h 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/sparsetools/fixed_size.h 2007-12-21 02:04:21 UTC (rev 3695)
@@ -18,16 +18,16 @@
};
template<class T>
-class Dot<0,T>
+class Dot<1,T>
{
public:
inline T operator()(const T * lhs, const T * rhs)
{
- return 0;
+ return *lhs * *rhs;
}
};
- template<int N, class T>
+template<int N, class T>
inline T dot(const T * lhs, const T * rhs)
{
Dot<N,T> d;
@@ -35,5 +35,4 @@
}
-
#endif
Modified: trunk/scipy/sparse/sparsetools/sparsetools.h
===================================================================
--- trunk/scipy/sparse/sparsetools/sparsetools.h 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/sparsetools/sparsetools.h 2007-12-21 02:04:21 UTC (rev 3695)
@@ -727,49 +727,88 @@
+//template <class I, class T>
+//void bsr_tocsr(const I n_brow,
+// const I n_bcol,
+// const I R,
+// const I C,
+// const I Ap[],
+// const I Aj[],
+// const T Ax[],
+// I Bp[],
+// I Bj[]
+// T Bx[])
+//{
+// const I RC = R*C;
+//
+// for(I brow = 0; brow < n_brow; brow++){
+// I row_size = C * (Ap[brow + 1] - Ap[brow]);
+// for(I r = 0; r < R; r++){
+// Bp[R*brow + r] = RC * Ap[brow] + r * row_size
+// }
+// }
+//}
template <class I, class T, int R, int C>
-void bsr_matvec_fixed(const I n_row,
- const I n_col,
+void bsr_matvec_fixed(const I n_brow,
+ const I n_bcol,
const I Ap[],
const I Aj[],
const T Ax[],
const T Xx[],
T Yx[])
{
- for(I i = 0; i < n_row; i++) {
+ for(I i = 0; i < n_brow; i++) {
T r0 = 0;
T r1 = 0;
T r2 = 0;
T r3 = 0;
T r4 = 0;
T r5 = 0;
+ T r6 = 0;
+ T r7 = 0;
for(I jj = Ap[i]; jj < Ap[i+1]; jj++) {
I j = Aj[jj];
const T * base = Ax + jj*(R*C);
- if (R >= 1) r0 += dot<C,T>(base + 0*C, Xx + j*C);
- if (R >= 2) r1 += dot<C,T>(base + 1*C, Xx + j*C);
- if (R >= 3) r2 += dot<C,T>(base + 2*C, Xx + j*C);
- if (R >= 4) r3 += dot<C,T>(base + 3*C, Xx + j*C);
- if (R >= 5) r4 += dot<C,T>(base + 4*C, Xx + j*C);
- if (R >= 6) r5 += dot<C,T>(base + 5*C, Xx + j*C);
+ if (R > 0) r0 += dot<C>(base + 0*C, Xx + j*C);
+ if (R > 1) r1 += dot<C>(base + 1*C, Xx + j*C);
+ if (R > 2) r2 += dot<C>(base + 2*C, Xx + j*C);
+ if (R > 3) r3 += dot<C>(base + 3*C, Xx + j*C);
+ if (R > 4) r4 += dot<C>(base + 4*C, Xx + j*C);
+ if (R > 5) r5 += dot<C>(base + 5*C, Xx + j*C);
+ if (R > 6) r6 += dot<C>(base + 6*C, Xx + j*C);
+ if (R > 7) r7 += dot<C>(base + 7*C, Xx + j*C);
}
- if (R >= 1) Yx[R*i+0] = r0;
- if (R >= 2) Yx[R*i+1] = r1;
- if (R >= 3) Yx[R*i+2] = r2;
- if (R >= 4) Yx[R*i+3] = r3;
- if (R >= 5) Yx[R*i+4] = r4;
- if (R >= 6) Yx[R*i+5] = r5;
+ if (R > 0) Yx[R*i+0] = r0;
+ if (R > 1) Yx[R*i+1] = r1;
+ if (R > 2) Yx[R*i+2] = r2;
+ if (R > 3) Yx[R*i+3] = r3;
+ if (R > 4) Yx[R*i+4] = r4;
+ if (R > 5) Yx[R*i+5] = r5;
+ if (R > 6) Yx[R*i+6] = r6;
+ if (R > 7) Yx[R*i+7] = r7;
}
}
+#define F(X,Y) bsr_matvec_fixed<I,T,X,Y>
+/*
+ * Generate the table below with:
+ * out = ''
+ * N = 8
+ * for i in range(N):
+ * out += '{'
+ * for j in range(N-1):
+ * out += ' F(%d,%d),' % (i+1,j+1)
+ * out += ' F(%d,%d) },\n' % (i+1,j+2)
+ * out = out[:-2]
+ *
+ */
-
template <class I, class T>
-void bsr_matvec(const I n_row,
- const I n_col,
+void bsr_matvec(const I n_brow,
+ const I n_bcol,
const I R,
const I C,
const I Ap[],
@@ -778,45 +817,38 @@
const T Xx[],
T Yx[])
{
- /*
- //use fixed version for R,C <= 4,4
- if (R == 1){
- if (C == 1) { bsr_matvec_fixed<I,T,1,1>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 2) { bsr_matvec_fixed<I,T,1,2>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 3) { bsr_matvec_fixed<I,T,1,3>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 4) { bsr_matvec_fixed<I,T,1,4>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
+ assert(R > 0 && C > 0);
+
+ void (*dispatch[8][8])(I,I,const I*,const I*,const T*,const T*,T*) = \
+ {
+ { F(1,1), F(1,2), F(1,3), F(1,4), F(1,5), F(1,6), F(1,7), F(1,8) },
+ { F(2,1), F(2,2), F(2,3), F(2,4), F(2,5), F(2,6), F(2,7), F(2,8) },
+ { F(3,1), F(3,2), F(3,3), F(3,4), F(3,5), F(3,6), F(3,7), F(3,8) },
+ { F(4,1), F(4,2), F(4,3), F(4,4), F(4,5), F(4,6), F(4,7), F(4,8) },
+ { F(5,1), F(5,2), F(5,3), F(5,4), F(5,5), F(5,6), F(5,7), F(5,8) },
+ { F(6,1), F(6,2), F(6,3), F(6,4), F(6,5), F(6,6), F(6,7), F(6,8) },
+ { F(7,1), F(7,2), F(7,3), F(7,4), F(7,5), F(7,6), F(7,7), F(7,8) },
+ { F(8,1), F(8,2), F(8,3), F(8,4), F(8,5), F(8,6), F(8,7), F(8,8) }
+ };
+
+ if (R <= 8 && C <= 8){
+ dispatch[R-1][C-1](n_brow,n_bcol,Ap,Aj,Ax,Xx,Yx);
+ return;
}
- if (R == 2){
- if (C == 1) { bsr_matvec_fixed<I,T,2,1>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 2) { bsr_matvec_fixed<I,T,2,2>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 3) { bsr_matvec_fixed<I,T,2,3>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 4) { bsr_matvec_fixed<I,T,2,4>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- }
- if (R == 3){
- if (C == 1) { bsr_matvec_fixed<I,T,3,1>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 2) { bsr_matvec_fixed<I,T,3,2>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 3) { bsr_matvec_fixed<I,T,3,3>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 4) { bsr_matvec_fixed<I,T,3,4>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- }
- if (R == 4){
- if (C == 1) { bsr_matvec_fixed<I,T,4,1>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 2) { bsr_matvec_fixed<I,T,4,2>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 3) { bsr_matvec_fixed<I,T,4,3>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- if (C == 4) { bsr_matvec_fixed<I,T,4,4>(n_row,n_col,Ap,Aj,Ax,Xx,Yx); return; }
- } */
//otherwise use general method
- for(I i = 0; i < n_row; i++){
+ for(I i = 0; i < n_brow; i++){
Yx[i] = 0;
}
- for(I i = 0; i < n_row; i++){
+ for(I i = 0; i < n_brow; i++){
const T * A = Ax + R * C * Ap[i];
T * y = Yx + R * i;
for(I jj = Ap[i]; jj < Ap[i+1]; jj++){
const T * x = Xx + C*Aj[jj];
+ //TODO replace this with a proper matvec
for( I r = 0; r < R; r++ ){
T sum = 0;
for( I c = 0; c < C; c++ ){
@@ -829,6 +861,7 @@
}
}
}
+#undef F
Modified: trunk/scipy/sparse/sparsetools/sparsetools.py
===================================================================
--- trunk/scipy/sparse/sparsetools/sparsetools.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/sparsetools/sparsetools.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -302,24 +302,24 @@
def bsr_matvec(*args):
"""
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
signed char Ax, signed char Xx, signed char Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
unsigned char Ax, unsigned char Xx, unsigned char Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
short Ax, short Xx, short Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
int Ax, int Xx, int Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
long long Ax, long long Xx, long long Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
float Ax, float Xx, float Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
double Ax, double Xx, double Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
npy_cfloat_wrapper Ax, npy_cfloat_wrapper Xx,
npy_cfloat_wrapper Yx)
- bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj,
+ bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj,
npy_cdouble_wrapper Ax, npy_cdouble_wrapper Xx,
npy_cdouble_wrapper Yx)
"""
Modified: trunk/scipy/sparse/sparsetools/sparsetools_wrap.cxx
===================================================================
--- trunk/scipy/sparse/sparsetools/sparsetools_wrap.cxx 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/sparsetools/sparsetools_wrap.cxx 2007-12-21 02:04:21 UTC (rev 3695)
@@ -43397,24 +43397,24 @@
" npy_cdouble_wrapper Xx, npy_cdouble_wrapper Yx)\n"
""},
{ (char *)"bsr_matvec", _wrap_bsr_matvec, METH_VARARGS, (char *)"\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" signed char Ax, signed char Xx, signed char Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" unsigned char Ax, unsigned char Xx, unsigned char Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" short Ax, short Xx, short Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" int Ax, int Xx, int Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" long long Ax, long long Xx, long long Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" float Ax, float Xx, float Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" double Ax, double Xx, double Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" npy_cfloat_wrapper Ax, npy_cfloat_wrapper Xx, \n"
" npy_cfloat_wrapper Yx)\n"
- "bsr_matvec(int n_row, int n_col, int R, int C, int Ap, int Aj, \n"
+ "bsr_matvec(int n_brow, int n_bcol, int R, int C, int Ap, int Aj, \n"
" npy_cdouble_wrapper Ax, npy_cdouble_wrapper Xx, \n"
" npy_cdouble_wrapper Yx)\n"
""},
Modified: trunk/scipy/sparse/tests/test_base.py
===================================================================
--- trunk/scipy/sparse/tests/test_base.py 2007-12-20 21:42:31 UTC (rev 3694)
+++ trunk/scipy/sparse/tests/test_base.py 2007-12-21 02:04:21 UTC (rev 3695)
@@ -21,7 +21,8 @@
from numpy.testing import *
set_package_path()
from scipy.sparse import csc_matrix, csr_matrix, dok_matrix, \
- coo_matrix, lil_matrix, dia_matrix, extract_diagonal, speye
+ coo_matrix, lil_matrix, dia_matrix, bsr_matrix, \
+ extract_diagonal, speye
from scipy.linsolve import splu
restore_path()
@@ -302,9 +303,7 @@
def check_conversions(self):
- #TODO add bsr/bsc
- #for format in ['bsc','bsr','coo','csc','csr','dia','dok','lil']:
- for format in ['coo','csc','csr','dia','dok','lil']:
+ for format in ['bsr','coo','csc','csr','dia','dok','lil']:
a = self.datsp.asformat(format)
assert_equal(a.format,format)
assert_array_almost_equal(a.todense(), self.dat)
@@ -318,18 +317,17 @@
#TODO, add and test .todia(maxdiags)
pass
-# def check_tocompressedblock(self):
-# #TODO more extensively test .tobsc() and .tobsr() w/ blocksizes
-# x = array([[1,0,2,0],[0,0,0,0],[0,0,4,5]])
-# y = array([[0,1,2],[3,0,5]])
-# A = kron(x,y)
-# Asp = self.spmatrix(A)
-# for format in ['bsc','bsr']:
-# fn = getattr(Asp, 'to' + format )
-#
-# for X in [ 1, 2, 3, 6 ]:
-# for Y in [ 1, 2, 3, 4, 6, 12]:
-# assert_equal( fn(blocksize=(X,Y)).todense(), A)
+ def check_tocompressedblock(self):
+ x = array([[1,0,2,0],[0,0,0,0],[0,0,4,5]])
+ y = array([[0,1,2],[3,0,5]])
+ A = kron(x,y)
+ Asp = self.spmatrix(A)
+ for format in ['bsr']:
+ fn = getattr(Asp, 'to' + format )
+
+ for X in [ 1, 2, 3, 6 ]:
+ for Y in [ 1, 2, 3, 4, 6, 12]:
+ assert_equal( fn(blocksize=(X,Y)).todense(), A)
def check_transpose(self):
@@ -1222,33 +1220,60 @@
pass
#TODO add test
-#class TestBSR(_TestCommon, _TestArithmetic, NumpyTestCase):
-# spmatrix = bsr_matrix
-#
-# def check_constructor1(self):
-# indptr = array([0,2,2,4])
-# indices = array([0,2,2,3])
-# data = zeros((4,2,3))
-#
-# data[0] = array([[ 0, 1, 2],
-# [ 3, 0, 5]])
-# data[1] = array([[ 0, 2, 4],
-# [ 6, 0, 10]])
-# data[2] = array([[ 0, 4, 8],
-# [12, 0, 20]])
-# data[3] = array([[ 0, 5, 10],
-# [15, 0, 25]])
-#
-# A = kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
-#
-# Asp = bsr_matrix((data,indices,indptr),shape=(6,12))
-#
-# assert_equal(Asp.todense(),A)
-# #TODO add infer from shape example
-#
-#
-#
-#
+class TestBSR(_TestCommon, _TestArithmetic, NumpyTestCase):
+ spmatrix = bsr_matrix
+
+ def check_constructor1(self):
+ """check native BSR format constructor"""
+ indptr = array([0,2,2,4])
+ indices = array([0,2,2,3])
+ data = zeros((4,2,3))
+
+ data[0] = array([[ 0, 1, 2],
+ [ 3, 0, 5]])
+ data[1] = array([[ 0, 2, 4],
+ [ 6, 0, 10]])
+ data[2] = array([[ 0, 4, 8],
+ [12, 0, 20]])
+ data[3] = array([[ 0, 5, 10],
+ [15, 0, 25]])
+
+ A = kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
+ Asp = bsr_matrix((data,indices,indptr),shape=(6,12))
+ assert_equal(Asp.todense(),A)
+
+ #infer shape from arrays
+ Asp = bsr_matrix((data,indices,indptr))
+ assert_equal(Asp.todense(),A)
+
+ def check_constructor2(self):
+ """construct from dense"""
+
+ #test zero mats
+ for shape in [ (1,1), (5,1), (1,10), (10,4), (3,7), (2,1)]:
+ A = zeros(shape)
+ assert_equal(bsr_matrix(A).todense(),A)
+ A = zeros((4,6))
+ assert_equal(bsr_matrix(A,blocksize=(2,2)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)
+
+ A = kron( [[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]] )
+ assert_equal(bsr_matrix(A).todense(),A)
+ assert_equal(bsr_matrix(A,shape=(6,12)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(1,1)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(2,6)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(2,12)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(3,12)).todense(),A)
+ assert_equal(bsr_matrix(A,blocksize=(6,12)).todense(),A)
+
+ A = kron( [[1,0,2,0],[0,1,0,0],[0,0,0,0]], [[0,1,2],[3,0,5]] )
+ assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)
+
+
+
+
+
#class TestBSC(_TestCommon, _TestArithmetic, NumpyTestCase):
# spmatrix = bsc_matrix
#
More information about the Scipy-svn
mailing list