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

scipy-svn@scip... scipy-svn@scip...
Tue Oct 16 01:19:15 CDT 2007


Author: wnbell
Date: 2007-10-16 01:19:12 -0500 (Tue, 16 Oct 2007)
New Revision: 3440

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
add block option to adaptive SA


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-16 06:19:12 UTC (rev 3440)
@@ -105,13 +105,13 @@
 
         self.Rs = [] 
         
-        if self.A.shape[0] <= self.opts['coarse: max size']:
+        if self.A.shape[0] <= max_coarse:
             raise ValueError,'small matrices not handled yet'
        
         #first candidate 
-        x,AggOps = self.__initialization_stage(A,blocks=blocks,\
-                            max_levels=max_levels,max_coarse=max_coarse,\
-                            mu=mu,epsilon=epsilon) 
+        x,AggOps = self.__initialization_stage(A, blocks=blocks,\
+                                               max_levels=max_levels, max_coarse=max_coarse,\
+                                               mu=mu, epsilon=epsilon) 
         
         Ws = AggOps
 
@@ -134,7 +134,7 @@
         self.AggOps = AggOps
                 
 
-    def __initialization_stage(self,A,max_levels,max_coarse,mu,epsilon):
+    def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
         AggOps = []
         Is     = []
 
@@ -270,6 +270,9 @@
 from scipy import *
 from utils import diag_sparse
 from multilevel import poisson_problem1D,poisson_problem2D
+
+blocks = None
+
 A = poisson_problem2D(200)
 #A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
 #A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
@@ -280,16 +283,17 @@
 #D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
 #A = D * A * D
 
-#A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
+A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
+blocks = arange(A.shape[0]/2).repeat(2)
 
-asa = adaptive_sa_solver(A,max_candidates=1,mu=5)
+asa = adaptive_sa_solver(A,max_candidates=4,mu=12)
 scipy.random.seed(0)  #make tests repeatable
 x = rand(A.shape[0])
 b = A*rand(A.shape[0])
 
 
 print "solving"
-if False:
+if True:
     x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
 else:
     residuals = []

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-10-16 06:19:12 UTC (rev 3440)
@@ -21,34 +21,37 @@
         if blocks is None:
             Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
             A_filtered = csr_matrix((Sx,Sj,Sp),A.shape)
-
         else:
-            num_dofs   = A.shape[0]
-            num_blocks = blocks.max() + 1
-            
-            if num_dofs != len(blocks):
-                raise ValueError,'improper block specification'
-            
-            # for non-scalar problems, use pre-defined blocks in aggregation
-            # the strength of connection matrix is based on the Frobenius norms of the blocks
-            
-            B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
-            Bt = B.T.tocsr()
-            
-            #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_A,epsilon)       
-            S.data[:] = 1
+            A_filtered = A  #TODO subtract weak blocks from diagonal blocks?
 
-            Mask = B * S * Bt
+##            num_dofs   = A.shape[0]
+##            num_blocks = blocks.max() + 1
+##            
+##            if num_dofs != len(blocks):
+##                raise ValueError,'improper block specification'
+##            
+##            # for non-scalar problems, use pre-defined blocks in aggregation
+##            # the strength of connection matrix is based on the 1-norms of the blocks
+##            
+##            B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
+##            Bt = B.T.tocsr()
+##            
+##            #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_A,epsilon)       
+##            S.data[:] = 1
+##
+##            Mask = B * S * Bt
+##
+##            A_strong = A ** Mask
+##            #A_weak   = A - A_strong
+##            A_filtered = A_strong
 
-            A_filtered = A ** Mask
-
     return A_filtered          
 
 
-def sa_strong_connections(A,epsilon,blocks=None):
+def sa_strong_connections(A,epsilon):
     if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
     
     Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
@@ -66,7 +69,7 @@
         if num_dofs != len(blocks):
             raise ValueError,'improper block specification'
        
-        print "SA has blocks"
+        #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
        

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-16 06:19:12 UTC (rev 3440)
@@ -1,7 +1,7 @@
 from numpy.testing import *
 
 from numpy import sqrt,empty,ones,arange,array_split,eye,array, \
-                  zeros,diag,zeros_like,diff
+                  zeros,diag,zeros_like,diff,matrix
 from numpy.linalg import norm                  
 from scipy import rand
 from scipy.sparse import spdiags,csr_matrix,lil_matrix, \
@@ -51,19 +51,42 @@
     def check_sa_strong_connections(self):
         for A in self.cases:
             for epsilon in [0.0,0.1,0.5,1.0,10.0]:
+                S_expected = reference_sa_strong_connections(A,epsilon)
                 S_result = sa_strong_connections(A,epsilon)
-                S_expected = reference_sa_strong_connections(A,epsilon)
                 assert_almost_equal(S_result.todense(),S_expected.todense())
                 #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
 
     def check_sa_constant_interpolation(self):
         for A in self.cases:
             for epsilon in [0.0,0.1,0.5,1.0]:
-                S_result   = sa_constant_interpolation(A,epsilon)
                 S_expected = reference_sa_constant_interpolation(A,epsilon)
+                
+                S_result   = sa_constant_interpolation(A,epsilon,blocks=None)
                 assert_array_equal(S_result.todense(),S_expected.todense())
+               
+                #blocks=1...N should be the same as blocks=None
+                S_result   = sa_constant_interpolation(A,epsilon,blocks=arange(A.shape[0]))
+                assert_array_equal(S_result.todense(),S_expected.todense())
 
+        #check simple block examples
+        A = csr_matrix(arange(16).reshape(4,4))
+        A = A + A.T
+        blocks = array([0,0,1,1])
 
+        S_result   = sa_constant_interpolation(A,epsilon=0.0,blocks=blocks)
+        S_expected = matrix([1,1,1,1]).T
+        assert_array_equal(S_result.todense(),S_expected)
+
+        S_result   = sa_constant_interpolation(A,epsilon=0.5,blocks=blocks)
+        S_expected = matrix([1,1,1,1]).T
+        assert_array_equal(S_result.todense(),S_expected)
+
+        S_result   = sa_constant_interpolation(A,epsilon=2.0,blocks=blocks)
+        S_expected = matrix([[1,0],[1,0],[0,1],[0,1]])
+        assert_array_equal(S_result.todense(),S_expected)
+
+                  
+
 class TestFitCandidates(NumpyTestCase):
     def setUp(self):
         self.cases = []

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-10-16 06:19:12 UTC (rev 3440)
@@ -3,10 +3,11 @@
 
 import numpy
 import scipy
-from numpy import ravel,arange
+from numpy import ravel,arange,concatenate
 from scipy.linalg import norm
 from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
-                        csr_matrix,csc_matrix,extract_diagonal
+                        csr_matrix,csc_matrix,extract_diagonal, \
+                        coo_matrix
 
 
 def approximate_spectral_radius(A,tol=0.1,maxiter=20):



More information about the Scipy-svn mailing list