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

scipy-svn@scip... scipy-svn@scip...
Wed Oct 17 15:29:17 CDT 2007


Author: wnbell
Date: 2007-10-17 15:29:14 -0500 (Wed, 17 Oct 2007)
New Revision: 3444

Modified:
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
Log:
added tests cases for user-defined aggregation



Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-17 18:28:32 UTC (rev 3443)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-17 20:29:14 UTC (rev 3444)
@@ -77,9 +77,8 @@
         candidates = [ ones(A.shape[0]) ] # use constant vector
 
     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
+        while len(As) < max_levels and A.shape[0] > max_coarse:
+            P,candidates,blocks = sa_interpolation(A,candidates,epsilon*0.5**(len(As)-1),omega=omega,blocks=blocks)
 
             A = (P.T.tocsr() * A) * P     #galerkin operator
 
@@ -207,12 +206,12 @@
     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 = 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 = 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)
@@ -221,7 +220,7 @@
         def add_resid(x):
             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]
+        x_sol = linalg.cg(A,b,x0=x,maxiter=25,tol=1e-12,callback=add_resid)[0]
             
 
     residuals = array(residuals)/residuals[0]
@@ -233,3 +232,5 @@
 
 
 
+
+

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-10-17 18:28:32 UTC (rev 3443)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-10-17 20:29:14 UTC (rev 3444)
@@ -8,7 +8,7 @@
 import multigridtools
 
 __all__ = ['sa_filtered_matrix','sa_strong_connections','sa_constant_interpolation',
-           'sa_interpolation','sa_fit_candidates']
+           'sa_interpolation','sa_smoothed_prolongator','sa_fit_candidates']
 
 
 def sa_filtered_matrix(A,epsilon,blocks=None):
@@ -83,6 +83,7 @@
         Pj = Pj[blocks] #expand block aggregates into constituent dofs
         Pp = B.indptr
         Px = B.data
+        
     else:
         S = sa_strong_connections(A,epsilon)
         
@@ -140,7 +141,18 @@
 
     return Q,coarse_candidates
     
-    
+def sa_smoothed_prolongator(A,T,epsilon,omega,blocks=None):
+    A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
+
+    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)
+
+    # smooth tentative prolongator T
+    P = T - (D_inv_A*T)
+
+    return P
+
 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')
 
@@ -148,7 +160,7 @@
         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'
+            raise TypeError,'expected csr_matrix for argument AggOp'
         if A.shape[1] != AggOp.shape[0]:
             raise ValueError,'incompatible aggregation operator'
 
@@ -156,14 +168,19 @@
 
     A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
 
-    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)
+    P = sa_smoothed_prolongator(A,T,epsilon,omega,blocks)
 
-    # smooth tentative prolongator T
-    P = T - (D_inv_A*T)
-           
-    return P,coarse_candidates
+##    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)
+##
+##    # smooth tentative prolongator T
+##    P = T - (D_inv_A*T)
 
+    if blocks is not None:
+        blocks = arange(AggOp.shape[1]).repeat(len(candidates))
+          
+    return P,coarse_candidates,blocks
 
 
+

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-17 18:28:32 UTC (rev 3443)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-17 20:29:14 UTC (rev 3444)
@@ -12,7 +12,8 @@
 set_package_path()
 import scipy.sandbox.multigrid
 from scipy.sandbox.multigrid.sa import sa_strong_connections, sa_constant_interpolation, \
-                                        sa_interpolation, sa_fit_candidates
+                                        sa_interpolation, sa_fit_candidates, \
+                                        sa_smoothed_prolongator
 from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D, \
                                         smoothed_aggregation_solver
 from scipy.sandbox.multigrid.utils import diag_sparse
@@ -68,6 +69,15 @@
                 S_result   = sa_constant_interpolation(A,epsilon,blocks=arange(A.shape[0]))
                 assert_array_equal(S_result.todense(),S_expected.todense())
 
+        # two aggregates in 1D
+        A = poisson_problem1D(6)
+        AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),dims=(6,2))
+        candidates = [ones(6)]
+
+        T_result,coarse_candidates_result = sa_fit_candidates(AggOp,candidates)
+        T_expected = csr_matrix((sqrt(1.0/3.0)*ones(6),array([0,0,0,1,1,1]),arange(7)),dims=(6,2))
+        assert_almost_equal(T_result.todense(),T_expected.todense())
+
         #check simple block examples
         A = csr_matrix(arange(16).reshape(4,4))
         A = A + A.T
@@ -85,6 +95,32 @@
         S_expected = matrix([[1,0],[1,0],[0,1],[0,1]])
         assert_array_equal(S_result.todense(),S_expected)
 
+
+    def check_user_aggregation(self):
+        """check that the sa_interpolation accepts user-defined aggregates"""
+
+        user_cases = []
+
+        #simple 1d example w/ two aggregates
+        A = poisson_problem1D(6)
+        AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),dims=(6,2))
+        candidates = [ones(6)]
+        user_cases.append((A,AggOp,candidates))
+
+        #simple 1d example w/ two aggregates (not all nodes are aggregated)
+        A = poisson_problem1D(6)
+        AggOp = csr_matrix((ones(4),array([0,0,1,1]),array([0,1,1,2,3,3,4])),dims=(6,2))
+        candidates = [ones(6)]
+        user_cases.append((A,AggOp,candidates))
+
+        for A,AggOp,candidates in user_cases:
+            T,coarse_candidates_result = sa_fit_candidates(AggOp,candidates)
+
+            P_result = sa_interpolation(A,candidates,omega=4.0/3.0,AggOp=AggOp)[0]
+            P_expected = sa_smoothed_prolongator(A,T,epsilon=0.0,omega=4.0/3.0)
+
+            assert_almost_equal(P_result.todense(),P_expected.todense())
+
                   
 
 class TestFitCandidates(NumpyTestCase):
@@ -149,9 +185,7 @@
 
 
 
-
-
-class TestSASolver(NumpyTestCase):
+class TestSASolverPerformance(NumpyTestCase):
     def setUp(self):
         self.cases = []
 



More information about the Scipy-svn mailing list