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

scipy-svn@scip... scipy-svn@scip...
Tue Oct 9 16:30:34 CDT 2007


Author: wnbell
Date: 2007-10-09 16:30:30 -0500 (Tue, 09 Oct 2007)
New Revision: 3427

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
fixed bug in approximate_spectral_radius where A was assumed to be symmetric
adaptive SA works for symmetric diagonal scaled case now


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-09 21:30:30 UTC (rev 3427)
@@ -4,7 +4,7 @@
 
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
-from sa import sa_constant_interpolation
+from sa import sa_constant_interpolation,sa_fit_candidates
 #from utils import infinity_norm
 from utils import approximate_spectral_radius
 
@@ -144,7 +144,7 @@
     Ps = []
 
     for W in Ws:
-        P,x = fit_candidates(W,x)
+        P,x = sa_fit_candidates(W,x)
         I   = smoothed_prolongator(P,A)  
         A   = I.T.tocsr() * A * I
         As.append(A)
@@ -301,7 +301,8 @@
         while len(AggOps) + 1 < max_levels  and A_l.shape[0] > max_coarse:
             #W_l   = sa_constant_interpolation(A_l,epsilon=0.08*0.5**(len(AggOps)-1))           #step 4b  #TEST
             W_l   = sa_constant_interpolation(A_l,epsilon=0)              #step 4b
-            P_l,x = fit_candidate(W_l,x)                                  #step 4c  
+            P_l,x = sa_fit_candidates(W_l,[x])                            #step 4c  
+            x = x[0]  #TODO make sa_fit_candidates accept a single x
             I_l   = smoothed_prolongator(P_l,A_l)                         #step 4d
             A_l   = I_l.T.tocsr() * A_l * I_l                             #step 4e
             
@@ -337,15 +338,15 @@
 from scipy import *
 from utils import diag_sparse
 from multilevel import poisson_problem1D,poisson_problem2D
-A = poisson_problem2D(50)
+A = poisson_problem2D(200)
 #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()
 #A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
 
 #A = A*A
-#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-#A = D * A * D
+D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+A = D * A * D
 #A = io.mmread("nos2.mtx").tocsr()
 asa = adaptive_sa_solver(A,max_candidates=1)
 #x = arange(A.shape[0]).astype('d') + 1

Modified: trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h	2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h	2007-10-09 21:30:30 UTC (rev 3427)
@@ -42,18 +42,26 @@
 
     T eps_Aii = epsilon*epsilon*diags[i];
 
+    T weak_sum = 0.0;
+
     for(int jj = row_start; jj < row_end; jj++){
       const int   j = Aj[jj];
       const T   Aij = Ax[jj];
 
-      if(i == j){continue;}
+      if(i == j){continue;} //skip diagonal until end of row
 
       //  |A(i,j)| < epsilon * sqrt(|A(i,i)|*|A(j,j)|) 
       if(Aij*Aij >= std::abs(eps_Aii * diags[j])){    
           Sj->push_back(j);
           Sx->push_back(Aij);
+      } else {
+          weak_sum += Aij;
       }
     }
+    //Add modified diagonal entry
+    Sj->push_back(i);
+    Sx->push_back(diags[i] + weak_sum); //filtered matrix
+
     Sp->push_back(Sj->size());
   }
 }

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-09 21:30:30 UTC (rev 3427)
@@ -4,8 +4,9 @@
 
 import scipy
 import numpy
-from numpy import zeros,zeros_like,array
+from numpy import ones,zeros,zeros_like,array
 from numpy.linalg import norm
+from scipy.linsolve import spsolve
 
 from sa import sa_interpolation
 from rs import rs_interpolation
@@ -19,8 +20,8 @@
     with standard 3-point finite difference stencil on a
     grid with N points.
     """
-    D = 2*numpy.ones(N)
-    O =  -numpy.ones(N)
+    D = 2*ones(N)
+    O =  -ones(N)
     return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate zeros
 
 def poisson_problem2D(N):
@@ -29,9 +30,9 @@
     with standard 5-point finite difference stencil on a
     square N-by-N grid.
     """
-    D = 4*numpy.ones(N*N)
-    T =  -numpy.ones(N*N)
-    O =  -numpy.ones(N*N)
+    D = 4*ones(N*N)
+    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
     
@@ -130,12 +131,12 @@
 
         #TODO change use of tol (relative tolerance) to agree with other iterative solvers
         A = self.As[0]
-        residuals = [scipy.linalg.norm(b-A*x)]
+        residuals = [ norm(b-A*x) ]
 
         while len(residuals) <= maxiter and residuals[-1]/residuals[0] > tol:
             self.__solve(0,x,b)
 
-            residuals.append(scipy.linalg.norm(b-A*x))
+            residuals.append( norm(b-A*x) )
 
             if callback is not None:
                 callback(x)
@@ -150,7 +151,7 @@
         A = self.As[lvl]
         
         if len(self.As) == 1:
-            x[:] = scipy.linsolve.spsolve(A,b)
+            x[:] = spsolve(A,b)
             return 
 
         self.presmoother(A,x,b)
@@ -162,7 +163,7 @@
         
         if lvl == len(self.As) - 2:
             #use direct solver on coarsest level
-            coarse_x[:] = scipy.linsolve.spsolve(self.As[-1],coarse_b)
+            coarse_x[:] = spsolve(self.As[-1],coarse_b)
             #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:   

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-09 21:30:30 UTC (rev 3427)
@@ -1,43 +1,60 @@
 from numpy.testing import *
 
-from numpy import sqrt,empty,ones,arange,array_split,eye,array,zeros,diag
+from numpy import sqrt,empty,ones,arange,array_split,eye,array, \
+                  zeros,diag,zeros_like
+from numpy.linalg import norm                  
 from scipy import rand
-from scipy.sparse import spdiags,csr_matrix,lil_matrix
+from scipy.sparse import spdiags,csr_matrix,lil_matrix, \
+                         isspmatrix_csr,isspmatrix_csc,isspmatrix_coo, \
+                         isspmatrix_lil
 import numpy
 
 set_package_path()
 import scipy.sandbox.multigrid
 from scipy.sandbox.multigrid.sa import sa_strong_connections, sa_constant_interpolation, \
                                         sa_interpolation, sa_fit_candidates
-from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D
+from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D, \
+                                        smoothed_aggregation_solver
+from scipy.sandbox.multigrid.utils import diag_sparse
 restore_path()
 
 
+#def sparsity(A):
+#    A = A.copy()
+#
+#    if isspmatrix_csr(A) or isspmatrix_csc(A) or isspmatrix_coo(A):
+#        A.data[:] = 1
+#    elif isspmatrix_lil:
+#        for row in A.data:
+#            row[:] = [1]*len(row)
+#    else:
+#        raise ValueError,'expected csr,csc,coo, or lil'
+#
+#    return A
+
 def reference_sa_strong_connections(A,epsilon):
     A_coo = A.tocoo()
     S = lil_matrix(A.shape)
 
     for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
-        if i == j: continue #skip diagonal
+        if i == j or abs(v) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
+            S[i,j] += v
+        else:
+            S[i,i] += v
 
-        if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
-            S[i,j] = v
+    return S
 
-    return S.tocsr()
-
 class TestSAStrongConnections(NumpyTestCase):
-    def check_simple(self):
-        N = 4
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
-        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
-        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
-        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
-       
-        N = 100
-        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
-        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
-        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
-        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+#    def check_simple(self):
+#        N = 4
+#        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+#        assert_array_equal(sa_strong_connections(A,0.50).todense(),A.todense())   #all connections are strong
+#        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*A.todense()) #no connections are strong
+#       
+#        N = 100
+#        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+#        assert_array_equal(sa_strong_connections(A,0.50).todense(),A.todense())   #all connections are strong
+#        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*A.todense()) #no connections are strong
 
     def check_random(self):
         numpy.random.seed(0)
@@ -47,7 +64,8 @@
             for epsilon in [0.0,0.1,0.5,1.0,10.0]:
                 S_result = sa_strong_connections(A,epsilon)
                 S_expected = reference_sa_strong_connections(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
+                assert_almost_equal(S_result.todense(),S_expected.todense())
+                #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
 
     def check_poisson1D(self):
         for N in [2,3,5,7,10,11,19]:
@@ -56,6 +74,7 @@
                 S_result   = sa_strong_connections(A,epsilon)
                 S_expected = reference_sa_strong_connections(A,epsilon)
                 assert_array_equal(S_result.todense(),S_expected.todense())
+                #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
 
     def check_poisson2D(self):
         for N in [2,3,5,7,10,11]:
@@ -64,6 +83,7 @@
                 S_result   = sa_strong_connections(A,epsilon)
                 S_expected = reference_sa_strong_connections(A,epsilon)
                 assert_array_equal(S_result.todense(),S_expected.todense())
+                #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
 
 
 
@@ -209,7 +229,89 @@
 
 
 
+class TestSASolver(NumpyTestCase):
+    def setUp(self):
+        self.cases = []
 
+        #self.cases.append((poisson_problem1D(10),None))
+
+        self.cases.append((poisson_problem1D(500),None))
+        self.cases.append((poisson_problem2D(50),None))
+        
+    def check_basic(self):
+        """check that method converges at a reasonable rate"""
+
+        for A,candidates in self.cases:
+            ml = smoothed_aggregation_solver(A,candidates,max_coarse=10,max_levels=10)
+
+            numpy.random.seed(0) #make tests repeatable
+            
+            x = rand(A.shape[0])
+            b = A*rand(A.shape[0]) #zeros_like(x)
+            
+            x_sol,residuals = ml.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
+
+            avg_convergence_ratio = (residuals[-1]/residuals[0])**(1.0/len(residuals))
+            
+            assert(avg_convergence_ratio < 0.5)
+
+    def check_DAD(self):
+
+        for A,candidates in self.cases:
+
+            x = rand(A.shape[0])
+            b = A*rand(A.shape[0]) #zeros_like(x)
+            
+            D = diag_sparse(rand(A.shape[0]))
+            D_inv = diag_sparse(1.0/D.data)
+            DAD = D*A*D
+           
+            if candidates is None:
+                candidates = [ ones(A.shape[0]) ]
+           
+            DAD_candidates = [ (D_inv * c) for c in candidates ]
+            
+            #ml = smoothed_aggregation_solver(A,candidates,max_coarse=1,max_levels=2)
+
+            ml = smoothed_aggregation_solver(DAD,DAD_candidates,max_coarse=100,max_levels=2)
+
+            #print (D_inv*ml.Ps[0]).todense()
+            
+            x_sol,residuals = ml.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
+
+            avg_convergence_ratio = (residuals[-1]/residuals[0])**(1.0/len(residuals))
+            print avg_convergence_ratio
+
+            assert(avg_convergence_ratio < 0.5)
+        
+##    def check_DAD(self):
+##        """check that method is invariant to symmetric diagonal scaling (A -> DAD)"""
+##
+##        for A,A_candidates in self.cases:
+##            numpy.random.seed(0) #make tests repeatable
+##            
+##            x = rand(A.shape[0])
+##            b = zeros_like(x)
+##
+##            D = diag_sparse(rand(A.shape[0]))
+##            D_inv = diag_sparse(1.0/D.data)
+##            DAD = D*A*D
+##            
+##
+##            if A_candidates is None:
+##                A_candidates = [ ones(A.shape[0]) ]
+##           
+##            DAD_candidates = [ (D_inv * c) for c in A_candidates ]
+##            
+##            ml_A    = smoothed_aggregation_solver(A,     A_candidates, max_coarse=10, max_levels=10, epsilon=0.0)
+##            x_sol_A = ml_A.solve(b, x0=x, maxiter=1, tol=1e-12)
+##
+##            ml_DAD    = smoothed_aggregation_solver(DAD, DAD_candidates, max_coarse=10, max_levels=10, epsilon=0.0)
+##            x_sol_DAD = ml_DAD.solve(b, x0=D*x, maxiter=1, tol=1e-12)
+##            
+##            assert_almost_equal(x_sol_A, D_inv * x_sol_DAD)
+
+
 if __name__ == '__main__':
     NumpyTest().run()
 

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-10-09 21:30:30 UTC (rev 3427)
@@ -10,8 +10,8 @@
     """
     Approximate the spectral radius of a symmetric matrix using ARPACK
     """
-    from scipy.sandbox.arpack import eigen_symmetric
-    return eigen_symmetric(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0]
+    from scipy.sandbox.arpack import eigen
+    return eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0]
 
 
 



More information about the Scipy-svn mailing list