[Scipy-svn] r3825 - in trunk/scipy/sandbox/multigrid: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Jan 12 05:23:23 CST 2008


Author: wnbell
Date: 2008-01-12 05:23:18 -0600 (Sat, 12 Jan 2008)
New Revision: 3825

Added:
   trunk/scipy/sandbox/multigrid/dec_example.py
   trunk/scipy/sandbox/multigrid/simple_example.py
Removed:
   trunk/scipy/sandbox/multigrid/dec_test.py
   trunk/scipy/sandbox/multigrid/simple_test.py
Modified:
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
Log:
change names of example scripts


Copied: trunk/scipy/sandbox/multigrid/dec_example.py (from rev 3822, trunk/scipy/sandbox/multigrid/dec_test.py)

Deleted: trunk/scipy/sandbox/multigrid/dec_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/dec_test.py	2008-01-12 11:20:31 UTC (rev 3824)
+++ trunk/scipy/sandbox/multigrid/dec_test.py	2008-01-12 11:23:18 UTC (rev 3825)
@@ -1,241 +0,0 @@
-
-from scipy import *
-from scipy.sparse import *
-from pydec import *
-from pydec.multigrid.discrete_laplacian import boundary_hierarchy, discrete_laplacian_solver, hodge_solver
-
-from scipy.sandbox.multigrid import smoothed_aggregation_solver,multigridtools,multilevel_solver
-from scipy.sandbox.multigrid.adaptive import adaptive_sa_solver
-from scipy.sandbox.multigrid.sa import sa_smoothed_prolongator
-from scipy.sandbox.multigrid.utils import expand_into_blocks
-
-
-## Load mesh from file
-mesh_path = '../../../../../hirani_group/wnbell/meshes/'
-#mesh = read_mesh(mesh_path + 'rocket/rocket.xml')
-#mesh = read_mesh(mesh_path + 'genus3/genus3_168k.xml')
-#mesh = read_mesh(mesh_path + 'genus3/genus3_455k.xml')
-#mesh = read_mesh(mesh_path + '/torus/torus.xml')
-mesh = read_mesh(mesh_path + '/sq14tri/sq14tri.xml')
-for i in range(5):
-    mesh['vertices'],mesh['elements'] = loop_subdivision(mesh['vertices'],mesh['elements'])
-cmplx = simplicial_complex(mesh['vertices'],mesh['elements'])
-
-## Construct mesh manually
-#bitmap = ones((60,60),dtype='bool')
-#bitmap[1::10,1::10] = False
-#bitmap[100:150,100:400] = False
-#cmplx = regular_cube_complex(regular_cube_mesh(bitmap))
-
-def curl_curl_prolongator(D_nodal,vertices):
-    if not isspmatrix_csr(D_nodal):
-        raise TypeError('expected csr_matrix')
-
-    A = D_nodal.T.tocsr() * D_nodal
-    aggs = multigridtools.sa_get_aggregates(A.shape[0],A.indptr,A.indices)
-
-    num_edges = D_nodal.shape[0]
-    num_basis = vertices.shape[1]
-    num_aggs  = aggs.max() + 1
-
-    # replace with CSR + eliminate duplicates
-    #indptr  = (2*num_basis) * arange(num_edges+1)
-    ## same same
-    #csr_matrix((data,indices,indptr),shape=(num_edges,num_aggs))
-
-    row  = arange(num_edges).repeat(2*num_basis)
-    col  = (num_basis*aggs[D_nodal.indices]).repeat(num_basis)
-    col = col.reshape(-1,num_basis) + arange(num_basis)
-    col = col.reshape(-1)
-    data = tile(0.5 * (D_nodal*vertices),(1,2)).reshape(-1)
-
-    return coo_matrix((data,(row,col)),shape=(num_edges,num_basis*num_aggs)).tocsr()
-
-
-
-
-
-def whitney_innerproduct_cache(cmplx,k):
-    h = hash(cmplx.vertices.tostring()) ^ hash(cmplx.simplices.tostring()) ^ hash(k)
-
-    filename = "/home/nathan/.pydec/cache/whitney_" + str(h) + ".mtx"
-
-    try:
-        import pydec
-        M = pydec.io.read_array(filename)
-    except:
-        import pydec
-        M = whitney_innerproduct(cmplx,k)
-        pydec.io.write_array(filename,M)
-
-    return M
-
-
-
-def cube_innerproduct_cache(cmplx,k):
-    h = hash(cmplx.mesh.bitmap.tostring()) ^ hash(cmplx.mesh.bitmap.shape) ^ hash(k)
-
-    filename = "/home/nathan/.pydec/cache/cube_" + str(h) + ".mtx"
-
-    try:
-        import pydec
-        M = pydec.io.read_array(filename)
-    except:
-        import pydec
-        M = regular_cube_innerproduct(cmplx,k)
-        pydec.io.write_array(filename,M)
-
-    return M
-
-
-
-#solve d_k d_k problem for all reasonable k
-#from pylab import semilogy,show,xlabel,ylabel,legend,ylim,xlim
-#from matplotlib.font_manager import fontManager, FontProperties
-
-cochain_complex = cmplx.cochain_complex()
-
-for i in [1]: #range(len(cochain_complex)-1):
-    print "computing mass matrix"
-
-    if isinstance(cmplx,simplicial_complex):
-        Mi = whitney_innerproduct_cache(cmplx,i+1)
-    else:
-        Mi = regular_cube_innerproduct(cmplx,i+1)
-
-
-    dimension = mesh['vertices'].shape[1]
-
-    if True:
-
-        d0 = cmplx[0].d
-        d1 = cmplx[1].d
-
-        #A = (d1.T.tocsr() * d1 + d0 * d0.T.tocsr()).astype('d')
-        A = (d1.T.tocsr()  * d1).astype('d')
-
-        P = curl_curl_prolongator(d0,mesh['vertices'])
-
-        num_blocks = P.shape[1]/dimension
-        blocks = arange(num_blocks).repeat(dimension)
-
-        P = sa_smoothed_prolongator(A,P,epsilon=0,omega=4.0/3.0)
-
-        PAP = P.T.tocsr() * A * P
-
-        candidates = None
-        candidates = zeros((num_blocks,dimension,dimension))
-        for n in range(dimension):
-            candidates[:,n,n] = 1.0
-        candidates = candidates.reshape(-1,dimension)
-
-        ml = smoothed_aggregation_solver(PAP,epsilon=0.0,candidates=candidates,blocks=blocks)
-        #A = PAP
-        ml = multilevel_solver([A] + ml.As, [P] + ml.Ps)
-    else:
-
-        bh = boundary_hierarchy(cochain_complex)
-        while len(bh) < 3:
-            bh.coarsen()
-        print repr(bh)
-
-        N = len(cochain_complex) - 1
-
-        B =  bh[0][N - i].B
-
-        A = (B.T.tocsr() * B).astype('d')
-        #A = B.T.tocsr() * Mi * B
-
-        constant_prolongators = [lvl[N - i].I for lvl in bh[:-1]]
-
-        method = 'aSA'
-
-        if method == 'RS':
-            As = [A]
-            Ps = []
-            for T in constant_prolongators:
-                Ps.append( sa_smoothed_prolongator(As[-1],T,epsilon=0.0,omega=4.0/3.0) )
-                As.append(Ps[-1].T.tocsr() * As[-1] * Ps[-1])
-            ml = multilevel_solver(As,Ps)
-
-        else:
-            if method == 'BSA':
-                if i == 0:
-                    candidates = None
-                else:
-                    candidates = cmplx[0].d * mesh['vertices']
-                    K = candidates.shape[1]
-
-                    constant_prolongators = [constant_prolongators[0]] + \
-                            [expand_into_blocks(T,K,1).tocsr() for T in constant_prolongators[1:] ]
-
-                    ml = smoothed_aggregation_solver(A,candidates,aggregation=constant_prolongators)
-            elif method == 'aSA':
-                asa = adaptive_sa_solver(A,aggregation=constant_prolongators,max_candidates=dimension,epsilon=0.0)
-                ml = asa.solver
-            else:
-                raise ValuerError,'unknown method'
-
-        #ml = smoothed_aggregation_solver(A,candidates)
-
-    #x = d0 * mesh['vertices'][:,0]
-    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=50,tol=1e-12,return_residuals=True)
-    else:
-        residuals = []
-        def add_resid(x):
-            residuals.append(linalg.norm(b - A*x))
-        A.psolve = ml.psolve
-
-        from pydec import cg
-        x_sol = cg(A,b,x0=x,maxiter=40,tol=1e-8,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
-
-
-
-
-##candidates = None
-##blocks = None
-##
-##
-##
-##A = io.mmread('tests/sample_data/elas30_A.mtx').tocsr()
-##candidates = io.mmread('tests/sample_data/elas30_nullspace.mtx')
-##candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
-##blocks = arange(A.shape[0]/2).repeat(2)
-##
-##ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=10,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-12,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
-##

Copied: trunk/scipy/sandbox/multigrid/simple_example.py (from rev 3822, trunk/scipy/sandbox/multigrid/simple_test.py)

Deleted: trunk/scipy/sandbox/multigrid/simple_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/simple_test.py	2008-01-12 11:20:31 UTC (rev 3824)
+++ trunk/scipy/sandbox/multigrid/simple_test.py	2008-01-12 11:23:18 UTC (rev 3825)
@@ -1,22 +0,0 @@
-from scipy import rand
-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')
-A = mats['A'].tobsr(blocksize=(2,2))
-B = mats['B']
-
-#A = poisson_problem2D(500)
-#B = None
-
-start = clock()
-sa = smoothed_aggregation_solver(A,B=B)
-print "constructed solver in %s seconds" % (clock() - start)
-
-b = rand(A.shape[0])
-start = clock()
-x,residuals = sa.solve(b,return_residuals=True)
-print "solved in %s seconds" % (clock() - start)
-print 'residuals',residuals

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2008-01-12 11:20:31 UTC (rev 3824)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2008-01-12 11:23:18 UTC (rev 3825)
@@ -179,21 +179,7 @@
             assert_almost_equal(fine_candidates,Q*coarse_candidates)
             assert_almost_equal(Q*(Q.T*fine_candidates),fine_candidates)
 
-##            #aggregate one more level (to a single aggregate)
-##            K = coarse_candidates.shape[1]
-##            N = K*AggOp.shape[1]
-##            AggOp = csr_matrix((ones(N),zeros(N),arange(N + 1)),shape=(N,1)) #aggregate to a single point
-##            fine_candidates = coarse_candidates
-##
-##            #mask_candidate(AggOp,fine_candidates)  #not needed
-##
-##            #now check the coarser problem
-##            Q,coarse_candidates = sa_fit_candidates(AggOp,fine_candidates)
-##
-##            assert_almost_equal(fine_candidates,Q*coarse_candidates)
-##            assert_almost_equal(Q*(Q.T*fine_candidates),fine_candidates)
 
-
 class TestSASolverPerformance(TestCase):
     def setUp(self):
         self.cases = []



More information about the Scipy-svn mailing list