librsb is a high performance sparse matrix library implementing the Recursive Sparse Blocks format, which is especially well suited for multiplications in iterative methods on huge sparse matrices.
PyRSB is a Cython-based Python interface to librsb.
So far, PyRSB is a working prototype, so prospective users and collaborators feedback is sought; please contact me to feedback and help.
The following functionality is implemented:
- Initialization with
rsb.rsb_matrix()styled asscipy.sparse.csr_matrix(). - Conversion from
scipy.sparse.csr_matrix(). - Multiplication by vector/multivector.
- Rows/columns through
nr=a.shape()[0]/nr=a.shape()[1]. find(),find_block(),tril(),triu(),shape(),nnz().print'able.- PyRSB-Specific:
autotune(),do_print().
- If you have librsb installed:
makewill attempt building and testing. Make sure you havecython,scipy,numpy. installed. - If you want the
Makefileto build librsb (in this directory):make all-localwill attempt downloading librsb-1.2.0-rc7 from the web and building it here before building pyrsb. If the file is in place, it won't download it a second time. After that,make local-librsb-pyrsb(ormake lp) will build pyrsb using local librsb, then run it. make testwill test benchmark code usingtest.py(to compare speed to SciPy)make bwill also produce graphs (requiresgnuplot)
import numpy
import scipy
from scipy.sparse import csr_matrix
from rsb import rsb_matrix
V=[11.,12.,22.]
I=[ 0, 0, 1]
J=[ 0, 1, 1]
c=csr_matrix((V,(I,J)))
print(c)
# several constructor forms, as with csr_matrix:
a=rsb_matrix((V,(I,J)))
a=rsb_matrix((V,(I,J)),[3,3])
a=rsb_matrix((V, I,J))
a=rsb_matrix((V, I,J),sym='S') # symmetric example
print(a)
a=rsb_matrix( (4,4))
a=rsb_matrix(c)
nrhs=1 # set to nrhs>1 to multiply by multiple vectors at once
nr=a.shape()[0]
nc=a.shape()[1]
x=numpy.empty([nc,nrhs],dtype=scipy.double)
y=numpy.empty([nr,nrhs],dtype=scipy.double)
x[:,:]=1.0
y[:,:]=0.0
print(a)
print(x)
print(y)
import rsb # import operators
# a.autotune() # makes only sense for large matrices
# matrix-vector multiply
y=y+a*x; # equivalent to y=y+c*x
print(y)
del a
GPLv3+