[Scipy-svn] r3478 - trunk/scipy/sandbox/multigrid

scipy-svn@scip... scipy-svn@scip...
Wed Oct 31 20:31:40 CDT 2007


Author: wnbell
Date: 2007-10-31 20:31:37 -0500 (Wed, 31 Oct 2007)
New Revision: 3478

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/dec_test.py
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
Log:
added new method for curl-curl problem


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -1,18 +1,18 @@
 import numpy,scipy,scipy.sparse
 from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
-                  asarray, hstack, ascontiguousarray, isinf
+                  asarray, hstack, ascontiguousarray, isinf, dot
 from numpy.random import randn
 from scipy.sparse import csr_matrix,coo_matrix
 
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
 from sa import sa_constant_interpolation,sa_fit_candidates
-from utils import approximate_spectral_radius,hstack_csr,vstack_csr,expand_into_blocks
+from utils import approximate_spectral_radius,hstack_csr,vstack_csr,expand_into_blocks,diag_sparse
 
 
 
 def augment_candidates(AggOp, old_Q, old_R, new_candidate):
-    #TODO update P also
+    #TODO update P and A also
 
     K = old_R.shape[1]
 
@@ -124,7 +124,7 @@
 
 
 class adaptive_sa_solver:
-    def __init__(self, A, blocks=None, max_levels=10, max_coarse=100,\
+    def __init__(self, A, blocks=None, aggregation=None, max_levels=10, max_coarse=100,\
                  max_candidates=1, mu=5, epsilon=0.1):
 
         self.A = A
@@ -138,8 +138,9 @@
         x,AggOps = self.__initialization_stage(A, blocks = blocks, \
                                                max_levels = max_levels, \
                                                max_coarse = max_coarse, \
-                                               mu = mu, epsilon = epsilon)
-
+                                               mu = mu, epsilon = epsilon, \
+                                               aggregation = aggregation ) 
+        
         #create SA using x here
         As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)
 
@@ -147,8 +148,20 @@
             x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
 
             #TODO which is faster?
-            #As,Ps,Ts,Bs = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
-            B = hstack((Bs[0],x))
+            As,Ps,Ts,Bs = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
+
+            #B = hstack((Bs[0],x))
+            #As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+
+        #improve candidates?
+        if True:
+            print "improving candidates"
+            B = Bs[0]
+            for i in range(max_candidates):
+                B = B[:,1:]
+                As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+                x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)    
+                B = hstack((B,x))
             As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
 
         self.Ts = Ts
@@ -156,12 +169,12 @@
         self.AggOps = AggOps
         self.Bs = Bs
 
+    def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon,aggregation):
+        if aggregation is not None:
+            max_coarse = 0
+            max_levels = len(aggregation) + 1
 
-    def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
-        AggOps = []
-        Ps     = []
-
-        # aSA parameters
+        # aSA parameters 
         # mu      - number of test relaxation iterations
         # epsilon - minimum acceptable relaxation convergence factor
 
@@ -177,9 +190,14 @@
         #TODO test convergence rate here
 
         As = [A]
+        AggOps = []
+        Ps = []
 
-        while len(AggOps) + 1 < max_levels  and A_l.shape[0] > max_coarse:
-            W_l   = sa_constant_interpolation(A_l,epsilon=0,blocks=blocks) #step 4b
+        while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
+            if aggregation is None:
+                W_l   = sa_constant_interpolation(A_l,epsilon=0,blocks=blocks) #step 4b
+            else:
+                W_l   = aggregation[len(AggOps)]
             P_l,x = sa_fit_candidates(W_l,x)                             #step 4c
             I_l   = smoothed_prolongator(P_l,A_l)                          #step 4d
             A_l   = I_l.T.tocsr() * A_l * I_l                              #step 4e
@@ -286,8 +304,11 @@
     from multilevel import poisson_problem1D,poisson_problem2D
 
     blocks = None
+    aggregation = None
 
-    A = poisson_problem2D(100)
+    #A = poisson_problem2D(200,1e-2)
+    #aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ] 
+
     #A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
     #A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
     #A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
@@ -301,17 +322,17 @@
     blocks = arange(A.shape[0]/2).repeat(2)
 
     from time import clock; start = clock()
-    asa = adaptive_sa_solver(A,max_candidates=5,mu=6,blocks=blocks)
+    asa = adaptive_sa_solver(A,max_candidates=3,mu=5,blocks=blocks,aggregation=aggregation)
     print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
 
-    #scipy.random.seed(0)  #make tests repeatable
+    scipy.random.seed(0)  #make tests repeatable
     x = randn(A.shape[0])
     b = A*randn(A.shape[0])
     #b = zeros(A.shape[0])
 
 
     print "solving"
-    if True:
+    if False:
         x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
     else:
         residuals = []
@@ -352,11 +373,12 @@
         pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
         show()
 
+  
     for c in asa.Bs[0].T:
+        #plot2d(c)
         plot2d_arrows(c)
         print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
 
 
-
     ##W = asa.AggOps[0]*asa.AggOps[1]
     ##pcolor((W * rand(W.shape[1])).reshape((200,200)))

Modified: trunk/scipy/sandbox/multigrid/dec_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/dec_test.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/dec_test.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -1,10 +1,12 @@
 
 from scipy import *
+from scipy.sparse import *
 from pydec import *
-from pydec.multigrid import *
 from pydec.multigrid.discrete_laplacian import boundary_hierarchy, discrete_laplacian_solver, hodge_solver
 
-from scipy.sandbox.multigrid import smoothed_aggregation_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
 
 
@@ -13,8 +15,9 @@
 #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')
-for i in range(3):
+#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'])
 
@@ -24,8 +27,34 @@
 #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),dims=(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)),dims=(num_edges,num_basis*num_aggs)).tocsr()
+
+
+    
+
+
 def whitney_innerproduct_cache(cmplx,k):
     h = hash(cmplx.vertices.tostring()) ^ hash(cmplx.simplices.tostring()) ^ hash(k)
 
@@ -73,47 +102,83 @@
         Mi = whitney_innerproduct_cache(cmplx,i+1)
     else:
         Mi = regular_cube_innerproduct(cmplx,i+1)
+        
+        
+    dimension = mesh['vertices'].shape[1]
 
-    ##print "constructing solver"
-    ##ss = discrete_laplacian_solver(cochain_complex,len(cochain_complex)-i-1,innerproduct=Mi)
-    ##print ss
-    ##
-    ##print "solving"
-    ##x,res = ss.solve(b=zeros(ss.A.shape[0]),x0=rand(ss.A.shape[0]),return_residuals=True)
-
-    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
-    #A = B.T.tocsr() * Mi * B
-
-    constant_prolongators = [lvl[N - i].I for lvl in bh[:-1]]
-
-    if i == 0:
+    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:
-        #candidates = [ones(A.shape[0])]
 
-        #TODO test
-        candidates = []
-        for coord in range(mesh['vertices'].shape[1]):
-            candidates.append( bh[0][N-i+1].B * mesh['vertices'][:,coord] )
+        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'
 
-        K = len(candidates)
+        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:] ]
 
-        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)
 
-
-    ml = smoothed_aggregation_solver(A,candidates,aggregation=constant_prolongators)
-    #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])
@@ -125,9 +190,11 @@
         def add_resid(x):
             residuals.append(linalg.norm(b - A*x))
         A.psolve = ml.psolve
-        x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
 
+        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

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -24,14 +24,14 @@
     O =  -ones(N)
     return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate explicit zeros
 
-def poisson_problem2D(N):
+def poisson_problem2D(N,epsilon=1.0):
     """
     Return a sparse CSR matrix for the 2d poisson problem
     with standard 5-point finite difference stencil on a
     square N-by-N grid.
     """
-    D = 4*ones(N*N)
-    T =  -ones(N*N)
+    D = (2 + 2*epsilon)*ones(N*N)
+    T =  -epsilon * ones(N*N)
     O =  -ones(N*N)
     T[N-1::N] = 0
     return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate explicit zeros
@@ -208,7 +208,9 @@
             #use direct solver on coarsest level
             #TODO reuse factors for efficiency?
             coarse_x[:] = spsolve(self.As[-1],coarse_b).reshape(coarse_x.shape)
-            #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0]
+            #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0].reshape(coarse_x.shape)
+            #A_inv = asarray(scipy.linalg.pinv2(self.As[-1].todense()))
+            #coarse_x[:] = scipy.dot(A_inv,coarse_b)
             #print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
         else:
             self.__solve(lvl+1,coarse_x,coarse_b)
@@ -230,18 +232,17 @@
     from scipy import *
     candidates = None
     blocks = None
-    #A = poisson_problem2D(100)
+    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')
-    #candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
-    blocks = arange(A.shape[0]/2).repeat(2)
+    #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)
 
-    ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=100,max_levels=2)
+    ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0.08,max_coarse=100,max_levels=10)
     #ml = ruge_stuben_solver(A)
 
     x = rand(A.shape[0])

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-11-01 01:31:37 UTC (rev 3478)
@@ -182,7 +182,7 @@
     P = T - (D_inv_A*T)
 
     #S = (spidentity(A.shape[0]).tocsr() - D_inv_A) #TODO drop this?
-    #P = S * ( S * T)
+    #P = S *(S * ( S * T))
 
     return P
 
@@ -201,8 +201,6 @@
     T,coarse_candidates = sa_fit_candidates(AggOp,candidates)
     #T = AggOp #TODO test
 
-    A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
-
     P = sa_smoothed_prolongator(A,T,epsilon,omega,blocks)
 
     if blocks is not None:



More information about the Scipy-svn mailing list