[Scipy-svn] r3826 - in trunk/scipy: sandbox/multigrid sparse sparse/benchmarks sparse/tests
scipy-svn@scip...
scipy-svn@scip...
Sat Jan 12 09:41:35 CST 2008
Author: wnbell
Date: 2008-01-12 09:41:27 -0600 (Sat, 12 Jan 2008)
New Revision: 3826
Added:
trunk/scipy/sparse/benchmarks/
trunk/scipy/sparse/benchmarks/test_sparse.py
Removed:
trunk/scipy/sparse/tests/test_sparse.py
Modified:
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/simple_example.py
Log:
updated example
fixed SA BSR usage
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2008-01-12 11:23:18 UTC (rev 3825)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2008-01-12 15:41:27 UTC (rev 3826)
@@ -155,7 +155,7 @@
while len(As) < max_levels and A.shape[0] > max_coarse:
P,B = sa_interpolation(A,B,epsilon*0.5**(len(As)-1),omega=omega)
- A = (P.T.tocsr() * A) * P #galerkin operator
+ A = (P.T.asformat(P.format) * A) * P #galerkin operator
As.append(A)
Ps.append(P)
@@ -276,45 +276,3 @@
gauss_seidel(A,x,b,iterations=1,sweep="symmetric")
-
-if __name__ == '__main__':
- from scipy import *
- from scipy.sandbox.multigrid import *
- candidates = None
- A = poisson_problem2D(40,1e-2)
- #A = io.mmread("rocker_arm_surface.mtx").tocsr()
- #A = io.mmread("9pt-100x100.mtx").tocsr()
- #A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
- #A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
-
- #A = io.mmread('tests/sample_data/elas30_A.mtx').tocsr()
- #candidates = io.mmread('tests/sample_data/elas30_nullspace.mtx')
- #blocks = arange(A.shape[0]/2).repeat(2)
-
- mats = io.loadmat('/home/nathan/Work/ogroup/matrices/elasticity/simple_grid_2d/elasticity_50x50.mat')
- A = mats['A']
- candidates = mats['B']
-
- ml = smoothed_aggregation_solver(A,candidates,epsilon=0.08,max_coarse=100,max_levels=10)
- #ml = ruge_stuben_solver(A)
-
- x = rand(A.shape[0])
- b = zeros_like(x)
- #b = A*rand(A.shape[0])
-
- if True:
- x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-8,return_residuals=True)
- else:
- residuals = []
- def add_resid(x):
- residuals.append(linalg.norm(b - A*x))
- A.psolve = ml.psolve
- x_sol = linalg.cg(A,b,x0=x,maxiter=25,tol=1e-12,callback=add_resid)[0]
-
-
- residuals = array(residuals)/residuals[0]
- avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
- print "average convergence ratio",avg_convergence_ratio
- print "last convergence ratio",residuals[-1]/residuals[-2]
-
- print residuals
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2008-01-12 11:23:18 UTC (rev 3825)
+++ trunk/scipy/sandbox/multigrid/sa.py 2008-01-12 15:41:27 UTC (rev 3826)
@@ -170,6 +170,7 @@
A_filtered = sa_filtered_matrix(A,epsilon) #use filtered matrix for anisotropic problems
+ # TODO use scale_rows()
D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
D_inv_A = D_inv * A_filtered
D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
Modified: trunk/scipy/sandbox/multigrid/simple_example.py
===================================================================
--- trunk/scipy/sandbox/multigrid/simple_example.py 2008-01-12 11:23:18 UTC (rev 3825)
+++ trunk/scipy/sandbox/multigrid/simple_example.py 2008-01-12 15:41:27 UTC (rev 3826)
@@ -1,22 +1,28 @@
-from scipy import rand
+from scipy import *
from scipy.sandbox.multigrid.sa import *
from scipy.sandbox.multigrid import *
from scipy.sandbox.multigrid.utils import *
from time import clock
-mats = io.loadmat('/home/nathan/Work/ogroup/matrices/elasticity/simple_grid_2d/elasticity_50x50.mat')
+mats = io.loadmat('/home/nathan/Work/ogroup/matrices/elasticity/simple_grid_2d/elasticity_500x500.mat')
A = mats['A'].tobsr(blocksize=(2,2))
B = mats['B']
-#A = poisson_problem2D(500)
+#A = poisson_problem2D(50)
#B = None
start = clock()
sa = smoothed_aggregation_solver(A,B=B)
print "constructed solver in %s seconds" % (clock() - start)
-b = rand(A.shape[0])
+x0 = rand(A.shape[0])
+b = zeros_like(x0)
start = clock()
-x,residuals = sa.solve(b,return_residuals=True)
+x,residuals = sa.solve(b,x0=x0,return_residuals=True)
print "solved in %s seconds" % (clock() - start)
+
+residuals = array(residuals)/residuals[0]
+avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
+print "average convergence ratio",avg_convergence_ratio
+print "last convergence ratio",residuals[-1]/residuals[-2]
print 'residuals',residuals
Copied: trunk/scipy/sparse/benchmarks/test_sparse.py (from rev 3822, trunk/scipy/sparse/tests/test_sparse.py)
Deleted: trunk/scipy/sparse/tests/test_sparse.py
===================================================================
--- trunk/scipy/sparse/tests/test_sparse.py 2008-01-12 11:23:18 UTC (rev 3825)
+++ trunk/scipy/sparse/tests/test_sparse.py 2008-01-12 15:41:27 UTC (rev 3826)
@@ -1,286 +0,0 @@
-"""general tests and simple benchmarks for the sparse module"""
-
-import numpy
-from numpy import ones, array, asarray, empty
-
-import random
-from scipy.testing import *
-
-from scipy.sparse import csc_matrix, csr_matrix, dok_matrix, \
- coo_matrix, lil_matrix, dia_matrix, spidentity, spdiags, \
- spkron
-from scipy.linsolve import splu
-
-
-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):
- """
- Return a sparse matrix for the 2d poisson problem
- with standard 5-point finite difference stencil on a
- square N-by-N grid.
- """
- if N == 1:
- diags = asarray( [[4]],dtype=dtype)
- return dia_matrix((diags,[0]), shape=(1,1)).asformat(format)
-
- offsets = array([0,-N,N,-1,1])
-
- diags = empty((5,N**2),dtype=dtype)
-
- diags[0] = 4 #main diagonal
- diags[1:] = -1 #all offdiagonals
-
- diags[3,N-1::N] = 0 #first lower diagonal
- diags[4,N::N] = 0 #first upper diagonal
-
- return dia_matrix((diags,offsets),shape=(N**2,N**2)).asformat(format)
-
-import time
-class TestSparseTools(TestCase):
- """Simple benchmarks for sparse matrix module"""
-
- @dec.bench
- def test_arithmetic(self):
- 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) )
-
- print
- print ' Sparse Matrix Arithmetic'
- print '===================================================================='
- print ' var | name | shape | dtype | nnz '
- print '--------------------------------------------------------------------'
- fmt = ' %1s | %14s | %20s | %9s | %8d '
-
- for var,name,mat in matrices:
- name = name.center(14)
- shape = ("%s" % (mat.shape,)).center(20)
- dtype = mat.dtype.name.center(9)
- print fmt % (var,name,shape,dtype,mat.nnz)
-
- space = ' ' * 10
- print
- print space+' Timings'
- print space+'=========================================='
- print space+' format | operation | time (msec) '
- print space+'------------------------------------------'
- fmt = space+' %3s | %17s | %7.1f '
-
- for format in ['csr']:
- vars = dict( [(var,mat.asformat(format)) for (var,name,mat) in matrices ] )
- for X,Y in [ ('A','A'),('A','B'),('B','A'),('B','B') ]:
- x,y = vars[X],vars[Y]
- for op in ['__add__','__sub__','multiply','__div__','__mul__']:
- fn = getattr(x,op)
- fn(y) #warmup
-
- start = time.clock()
- iter = 0
- while iter < 5 or time.clock() < start + 1:
- fn(y)
- iter += 1
- end = time.clock()
-
- msec_per_it = 1000*(end - start)/float(iter)
- operation = (X + '.' + op + '(' + Y + ')').center(17)
- print fmt % (format,operation,msec_per_it)
-
-
- @dec.bench
- def test_sort(self):
- """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._has_sorted_indices = False
- A.sort_indices()
- 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) )
-
- @dec.bench
- def test_matvec(self):
- matrices = []
- matrices.append(('Identity', spidentity(10**4,format='dia')))
- matrices.append(('Identity', spidentity(10**4,format='csr')))
- matrices.append(('Poisson5pt', poisson2d(300,format='dia')))
- matrices.append(('Poisson5pt', poisson2d(300,format='csr')))
- matrices.append(('Poisson5pt', poisson2d(300,format='bsr')))
-
- A = spkron(poisson2d(150),ones((2,2))).tobsr(blocksize=(2,2))
- matrices.append( ('Block2x2', A.tocsr()) )
- matrices.append( ('Block2x2', A) )
-
- A = spkron(poisson2d(100),ones((3,3))).tobsr(blocksize=(3,3))
- matrices.append( ('Block3x3', A.tocsr()) )
- matrices.append( ('Block3x3', A) )
-
- print
- print ' Sparse Matrix Vector Product'
- print '=================================================================='
- print ' type | name | shape | nnz | MFLOPs '
- print '------------------------------------------------------------------'
- fmt = ' %3s | %12s | %20s | %8d | %6.1f '
-
- for name,A in matrices:
- x = ones(A.shape[1],dtype=A.dtype)
-
- y = A*x #warmup
-
- start = time.clock()
- iter = 0
- while iter < 5 or time.clock() < start + 1:
- try:
- #avoid creating y if possible
- A.matvec(x,y)
- except:
- y = A*x
- iter += 1
- end = time.clock()
-
- del y
-
- name = name.center(12)
- shape = ("%s" % (A.shape,)).center(20)
- MFLOPs = (2*A.nnz*iter/(end-start))/float(1e6)
-
- print fmt % (A.format,name,shape,A.nnz,MFLOPs)
-
- @dec.bench
- def test_construction(self):
- """build matrices by inserting single values"""
- matrices = []
- matrices.append( ('Empty',csr_matrix((10000,10000))) )
- matrices.append( ('Identity',spidentity(10000)) )
- matrices.append( ('Poisson5pt', poisson2d(100)) )
-
- print
- print ' Sparse Matrix Construction'
- print '===================================================================='
- print ' type | name | shape | nnz | time (sec) '
- print '--------------------------------------------------------------------'
- fmt = ' %3s | %12s | %20s | %8d | %6.4f '
-
- for name,A in matrices:
- A = A.tocoo()
-
- for format in ['lil','dok']:
-
- start = time.clock()
-
- iter = 0
- while time.clock() < start + 0.5:
- T = eval(format + '_matrix')(A.shape)
- for i,j,v in zip(A.row,A.col,A.data):
- T[i,j] = v
- iter += 1
- end = time.clock()
-
- del T
- name = name.center(12)
- shape = ("%s" % (A.shape,)).center(20)
-
- print fmt % (format,name,shape,A.nnz,(end-start)/float(iter))
-
- @dec.bench
- def test_conversion(self):
- A = poisson2d(100)
-
- formats = ['csr','csc','coo','lil','dok']
-
- print
- print ' Sparse Matrix Conversion'
- print '=========================================================='
- print ' format | tocsr() | tocsc() | tocoo() | tolil() | todok() '
- print '----------------------------------------------------------'
-
- for fromfmt in formats:
- base = getattr(A,'to' + fromfmt)()
-
- times = []
-
- for tofmt in formats:
- try:
- fn = getattr(base,'to' + tofmt)
- except:
- times.append(None)
- else:
- x = fn() #warmup
- start = time.clock()
- iter = 0
- while time.clock() < start + 0.2:
- x = fn()
- iter += 1
- end = time.clock()
- del x
- times.append( (end - start)/float(iter))
-
- output = " %3s " % fromfmt
- for t in times:
- if t is None:
- output += '| n/a '
- else:
- output += '| %5.1fms ' % (1000*t)
- print output
-
-
-class TestLarge(TestCase):
- def test_large(self):
- # Create a 100x100 matrix with 100 non-zero elements
- # and play around with it
- #TODO move this out of Common since it doesn't use spmatrix
- random.seed(0)
- A = dok_matrix((100,100))
- for k in range(100):
- i = random.randrange(100)
- j = random.randrange(100)
- A[i,j] = 1.
- csr = A.tocsr()
- csc = A.tocsc()
- csc2 = csr.tocsc()
- coo = A.tocoo()
- csr2 = coo.tocsr()
- assert_array_equal(A.transpose().todense(), csr.transpose().todense())
- assert_array_equal(csc.todense(), csr.todense())
- assert_array_equal(csr.todense(), csr2.todense())
- assert_array_equal(csr2.todense().transpose(), coo.transpose().todense())
- assert_array_equal(csr2.todense(), csc2.todense())
- csr_plus_csc = csr + csc
- csc_plus_csr = csc + csr
- assert_array_equal(csr_plus_csc.todense(), (2*A).todense())
- assert_array_equal(csr_plus_csc.todense(), csc_plus_csr.todense())
-
-
-if __name__ == "__main__":
- unittest.main()
-
More information about the Scipy-svn
mailing list