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

scipy-svn@scip... scipy-svn@scip...
Mon Oct 15 15:33:26 CDT 2007


Author: wnbell
Date: 2007-10-15 15:33:20 -0500 (Mon, 15 Oct 2007)
New Revision: 3437

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
Log:
make SA work properly w/ blocks



Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-15 19:23:07 UTC (rev 3436)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-15 20:33:20 UTC (rev 3437)
@@ -16,14 +16,14 @@
      
     """
 
-    #candidate prolongator (assumes every value from x is used)
+    #candidate prolongator (assumes every value from x is used)  #TODO permit gaps
     X = csr_matrix((x_l,W_l.indices,W_l.indptr),dims=W_l.shape,check=False)  
     
     R = (P_l.T.tocsr() * X)  # R has at most 1 nz per row
     X = X - P_l*R            # othogonalize X against P_l
     
-    #DROP REDUNDANT COLUMNS FROM P (AND R?) HERE (NULL OUT R ACCORDINGLY?)    
-    #REMOVE CORRESPONDING COLUMNS FROM W_l AND ROWS FROM A_m ALSO
+    #TODO DROP REDUNDANT COLUMNS FROM P (AND R?) HERE (NULL OUT R ACCORDINGLY?)    
+    #TODO REMOVE CORRESPONDING COLUMNS FROM W_l AND ROWS FROM A_m ALSO
     W_l_new = W_l
     W_m_new = W_m
 
@@ -59,7 +59,7 @@
     D = diag_sparse(A)
     D_inv_A = diag_sparse(1.0/D)*A
     omega = 4.0/(3.0*approximate_spectral_radius(D_inv_A))
-    print "spectral radius",approximate_spectral_radius(D_inv_A)
+    print "spectral radius",approximate_spectral_radius(D_inv_A) #TODO remove this
     D_inv_A *= omega
 
     return P - D_inv_A*P
@@ -124,8 +124,9 @@
             x = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps,mu=mu)    
             
             self.candidates.append(x)
-           
-            #As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)             
+            
+            #TODO which is faster?
+            #As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)  
             As,Is,Ps = sa_hierarchy(A,AggOps,self.candidates)
 
         self.Ps = Ps 

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-15 19:23:07 UTC (rev 3436)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-15 20:33:20 UTC (rev 3437)
@@ -22,7 +22,7 @@
     """
     D = 2*ones(N)
     O =  -ones(N)
-    return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate zeros
+    return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate explicit zeros
 
 def poisson_problem2D(N):
     """
@@ -34,7 +34,7 @@
     T =  -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 zeros
+    return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate explicit zeros
     
 
 def ruge_stuben_solver(A,max_levels=10,max_coarse=500):
@@ -60,7 +60,7 @@
         
     return multilevel_solver(As,Ps)
 
-def smoothed_aggregation_solver(A,candidates=None,blocks=None,max_levels=10,max_coarse=500,epsilon=0.08,omega=4.0/3.0):
+def smoothed_aggregation_solver(A,candidates=None,blocks=None,aggregation=None,max_levels=10,max_coarse=500,epsilon=0.08,omega=4.0/3.0):
     """
     Create a multilevel solver using Smoothed Aggregation (SA)
 
@@ -75,15 +75,25 @@
     
     if candidates is None:
         candidates = [ ones(A.shape[0]) ] # use constant vector
-        
-    while len(As) < max_levels  and A.shape[0] > max_coarse:
-        P,candidates = sa_interpolation(A,candidates,epsilon*0.5**(len(As)-1),omega=omega,blocks=blocks)
-        #P = sa_interpolation(A,epsilon=0.0)
 
-        A = (P.T.tocsr() * A) * P     #galerkin operator
+    if aggregation is None:
+        while len(As) < max_levels  and A.shape[0] > max_coarse:
+            P,candidates = sa_interpolation(A,candidates,epsilon*0.5**(len(As)-1),omega=omega,blocks=blocks)
+            blocks = None #only used for 1st level
 
-        As.append(A)
-        Ps.append(P)
+            A = (P.T.tocsr() * A) * P     #galerkin operator
+
+            As.append(A)
+            Ps.append(P)
+    else:
+        #use user-defined aggregation
+        for AggOp in aggregation:
+            P,candidates = sa_interpolation(A,candidates,omega=omega,AggOp=AggOp)
+
+            A = (P.T.tocsr() * A) * P     #galerkin operator
+
+            As.append(A)
+            Ps.append(P)
         
     return multilevel_solver(As,Ps)
 
@@ -163,7 +173,7 @@
         
         if lvl == len(self.As) - 2:
             #use direct solver on coarsest level
-            coarse_x[:] = spsolve(self.As[-1],coarse_b)
+            coarse_x[:] = spsolve(self.As[-1],coarse_b)  #TODO reuse factors for efficiency?
             #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0]
             #print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
         else:   
@@ -185,17 +195,19 @@
 if __name__ == '__main__':
     from scipy import *
     candidates = None
-    A = poisson_problem2D(100)
+    blocks = None
+    #A = poisson_problem2D(100)
     #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]) ]
-    
-    ml = smoothed_aggregation_solver(A,candidates,max_coarse=10,max_levels=10)
+    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,max_coarse=10,max_levels=10)
     #ml = ruge_stuben_solver(A)
 
     x = rand(A.shape[0])
@@ -210,8 +222,12 @@
             residuals.append(linalg.norm(b - A*x))
         A.psolve = ml.psolve
         x_sol = linalg.cg(A,b,x0=x,maxiter=12,tol=1e-100,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	2007-10-15 19:23:07 UTC (rev 3436)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-10-15 20:33:20 UTC (rev 3437)
@@ -1,6 +1,7 @@
 import scipy
 import numpy
-from numpy import array,arange,ones,zeros,sqrt,isinf,asarray,empty,diff
+from numpy import array,arange,ones,zeros,sqrt,isinf,asarray,empty,diff,\
+                  ascontiguousarray
 from scipy.sparse import csr_matrix,isspmatrix_csr
 
 from utils import diag_sparse,approximate_spectral_radius
@@ -23,7 +24,7 @@
 
         else:
             num_dofs   = A.shape[0]
-            num_blocks = blocks.max()
+            num_blocks = blocks.max() + 1
             
             if num_dofs != len(blocks):
                 raise ValueError,'improper block specification'
@@ -34,11 +35,10 @@
             B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
             Bt = B.T.tocsr()
             
-            #Frobenius norms of blocks entries of A
-            #TODO change to 1-norm ?
-            Block_Frob = Bt * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B 
+            #1-norms of blocks entries of A
+            Block_A = Bt * csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape) * B 
     
-            S = sa_strong_connections(Block_Frob,epsilon)       
+            S = sa_strong_connections(Block_A,epsilon)       
             S.data[:] = 1
 
             Mask = B * S * Bt
@@ -61,20 +61,20 @@
     
     if blocks is not None:
         num_dofs   = A.shape[0]
-        num_blocks = blocks.max()
+        num_blocks = blocks.max() + 1
         
         if num_dofs != len(blocks):
             raise ValueError,'improper block specification'
-        
+       
+        print "SA has blocks"
         # for non-scalar problems, use pre-defined blocks in aggregation
         # the strength of connection matrix is based on the Frobenius norms of the blocks
        
-        #TODO change this to matrix 1 norm?
         B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
-        #Frobenius norms of blocks entries of A
-        Block_Frob = B.T.tocsr() * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B 
+        #1-norms of blocks entries of A
+        Block_A = B.T.tocsr() * csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape) * B 
 
-        S = sa_strong_connections(Block_Frob,epsilon)
+        S = sa_strong_connections(Block_A,epsilon)
     
         Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
         Pj = Pj[blocks] #expand block aggregates into constituent dofs
@@ -108,7 +108,6 @@
         c = c[diff(AggOp.indptr) == 1]  #eliminate DOFs that aggregation misses
         X = csr_matrix((c,AggOp.indices,AggOp.indptr),dims=AggOp.shape)
        
-
         #orthogonalize X against previous
         for j,A in enumerate(candidate_matrices):
             D_AtX = csr_matrix((A.data*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X            
@@ -134,15 +133,22 @@
         Q_data[i::K] = X.data
     Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
 
-    coarse_candidates = [array(R[:,i]) for i in range(K)]
+    coarse_candidates = [ascontiguousarray(R[:,i]) for i in range(K)]
 
     return Q,coarse_candidates
     
     
-def sa_interpolation(A,candidates,epsilon,omega=4.0/3.0,blocks=None):
+def sa_interpolation(A,candidates,epsilon=0.0,omega=4.0/3.0,blocks=None,AggOp=None):
     if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
 
-    AggOp  = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
+    if AggOp is None:
+        AggOp = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
+    else:
+        if not isspmatrix_csr(AggOp):
+            raise TypeError,'aggregation is specified by a list of csr_matrix objects'
+        if A.shape[1] != AggOp.shape[0]:
+            raise ValueError,'incompatible aggregation operator'
+
     T,coarse_candidates = sa_fit_candidates(AggOp,candidates)  
 
     A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems



More information about the Scipy-svn mailing list