[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