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

scipy-svn@scip... scipy-svn@scip...
Sat Jan 19 13:43:07 CST 2008


Author: wnbell
Date: 2008-01-19 13:43:04 -0600 (Sat, 19 Jan 2008)
New Revision: 3852

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/examples/adaptive.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
updated aSA


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-19 19:43:04 UTC (rev 3852)
@@ -6,7 +6,8 @@
 from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
                   asarray, hstack, ascontiguousarray, isinf, dot
 from numpy.random import randn, rand
-from scipy.sparse import csr_matrix,coo_matrix
+from scipy.linalg import norm
+from scipy.sparse import csr_matrix, coo_matrix, bsr_matrix
 
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
@@ -40,7 +41,7 @@
     for AggOp in AggOps:
         P,B = sa_fit_candidates(AggOp,B)
         I   = sa_smoothed_prolongator(A,P)
-        A   = I.T.tocsr() * A * I
+        A   = I.T.asformat(I.format) * A * I
         As.append(A)
         Ts.append(P)
         Ps.append(I)
@@ -99,18 +100,22 @@
     """
     
     if A.shape[0] <= max_coarse:
-        raise ValueError,'small matrices not handled yet'
+        return multilevel_solver( [A], [] )
 
     #first candidate
     x,AggOps = asa_initial_setup_stage(A, max_levels = max_levels, \
             max_coarse = max_coarse, mu = mu, epsilon = epsilon, \
             aggregation = aggregation )
 
+    #TODO make sa_fit_candidates work for small Bs
+    x /= norm(x)
+    
     #create SA using x here
     As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)
 
     for i in range(max_candidates - 1):
         x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu=mu)
+        x /= norm(x)
 
         B = hstack((Bs[0],x))
         As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
@@ -122,7 +127,7 @@
         for i in range(max_candidates):
             B = B[:,1:]
             As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
-            x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu=mu)
+            x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu)
             B = hstack((B,x))
         As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
 
@@ -159,7 +164,7 @@
             W_l = aggregation[len(AggOps)]
         P_l,x = sa_fit_candidates(W_l,x)                   #step 4c
         I_l   = sa_smoothed_prolongator(A_l,P_l)           #step 4d
-        A_l   = I_l.T.tocsr() * A_l * I_l                  #step 4e
+        A_l   = I_l.T.asformat(I_l.format) * A_l * I_l                  #step 4e
 
         AggOps.append(W_l)
         Ps.append(I_l)
@@ -189,7 +194,7 @@
 def asa_general_setup_stage(As, Ps, Ts, Bs, AggOps, mu):
     A = As[0]
 
-    x = randn(A.shape[0],1)
+    x = rand(A.shape[0],1)
     b = zeros_like(x)
 
     x = multilevel_solver(As,Ps).solve(b, x0=x, tol=1e-10, maxiter=mu)
@@ -204,15 +209,20 @@
         K    = P.blocksize[0]
         bnnz = P.indptr[-1]
         data = zeros( (bnnz, K+1, K), dtype=P.dtype )
-        data[:,:-1,:-1] = P.data
+        data[:,:-1,:] = P.data
         return bsr_matrix( (data, P.indices, P.indptr), shape=( (K+1)*(M/K), N) )
 
     for i in range(len(As) - 2):
-        B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
-        T,R = sa_fit_candidates(AggOp,B)
+        B_old = Bs[i]
+        B = zeros( (x.shape[0], B_old.shape[1] + 1), dtype=x.dtype)
 
+        B[:B_old.shape[0],:B_old.shape[1]] = B_old
+        B[:,-1] = x.reshape(-1)
+
+        T,R = sa_fit_candidates(AggOps[i],B)
+
         P = sa_smoothed_prolongator(A,T)
-        A = P.T.tocsr() * A * P
+        A = P.T.asformat(P.format) * A * P
 
         temp_Ps.append(P)
         temp_As.append(A)
@@ -242,7 +252,7 @@
 #            T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
 #
 #            P = sa_smoothed_prolongator(A,T)
-#            A = P.T.tocsr() * A * P
+#            A = P.T.asformat(P.format) * A * P
 #
 #            new_As.append(A)
 #            new_Ps.append(P)

Modified: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-19 19:43:04 UTC (rev 3852)
@@ -3,21 +3,25 @@
 from scipy.sandbox.multigrid.utils import diag_sparse
 from scipy.sandbox.multigrid.gallery import poisson, linear_elasticity
 from scipy.sandbox.multigrid.adaptive import adaptive_sa_solver
+from scipy.sandbox.multigrid import smoothed_aggregation_solver
 
-#A,B = linear_elasticity( (100,100) )
+A,B = linear_elasticity( (20,20) )
 
-A = poisson( (200,200), format='csr' )
+#A = poisson( (200,200), format='csr' )
 
 #A = poisson( (200,200), spacing=(1,1e-2) )  #anisotropic 
-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
 
 
 
 from time import clock; start = clock()
 
-asa = adaptive_sa_solver(A, max_levels=2, max_candidates=1, mu=10)
+#asa = smoothed_aggregation_solver(A,B,max_levels=2)
+#Bs = [B]
 
+asa,Bs = adaptive_sa_solver(A, max_levels=2, max_coarse=10, max_candidates=3, mu=15)
+
 print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
 
 random.seed(0)  #make tests repeatable
@@ -27,14 +31,14 @@
 
 
 print "solving"
-if True:
-    x_sol,residuals = asa.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
+if False:
+    x_sol,residuals = asa.solve(b,x0=x,maxiter=20,tol=1e-10,return_residuals=True)
 else:
     residuals = []
     def add_resid(x):
         residuals.append(linalg.norm(b - A*x))
     A.psolve = asa.psolve
-    x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
+    x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-10,callback=add_resid)[0]
 
 residuals = array(residuals)/residuals[0]
 
@@ -69,10 +73,10 @@
     show()
 
 
-#for c in asa.Bs[0].T:
-#    plot2d(c)
-#    #plot2d_arrows(c)
-#    print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+for c in Bs[0].T:
+    #plot2d(c)
+    plot2d_arrows(c)
+    print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
 
 #plot2d(x_sol)
 

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/sa.py	2008-01-19 19:43:04 UTC (rev 3852)
@@ -103,6 +103,13 @@
     if candidates.dtype != 'float32':
         candidates = asarray(candidates,dtype='float64')
 
+    if len(candidates.shape) != 2:
+        raise ValueError,'expected rank 2 array for argument B'
+
+    if candidates.shape[0] % AggOp.shape[0] != 0:
+        raise ValueError,'dimensions of AggOp %s and B %s are incompatible' % (AggOp.shape, B.shape)
+    
+
     K = candidates.shape[1] # number of near-nullspace candidates
     blocksize = candidates.shape[0] / AggOp.shape[0]
 
@@ -167,8 +174,11 @@
     A_filtered = sa_filtered_matrix(A,epsilon) #use filtered matrix for anisotropic problems
 
     # TODO use scale_rows()
-    D_inv    = diag_sparse(1.0/diag_sparse(A_filtered))
-    D_inv_A  = D_inv * A_filtered
+    D = diag_sparse(A_filtered)
+    D_inv = 1.0 / D
+    D_inv[D == 0] = 0
+
+    D_inv_A  = diag_sparse(D_inv) * A_filtered
     D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
 
     # smooth tentative prolongator T

Modified: trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_adaptive.py	2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/tests/test_adaptive.py	2008-01-19 19:43:04 UTC (rev 3852)
@@ -1,18 +1,56 @@
 from scipy.testing import *
 
+from numpy import arange, ones, zeros, array, eye, vstack, diff
+from scipy import rand
 from scipy.sparse import csr_matrix
-from numpy import arange,ones,zeros,array,eye,vstack,diff
 
 
 from scipy.sandbox.multigrid.sa import sa_fit_candidates
+from scipy.sandbox.multigrid import smoothed_aggregation_solver
 #from scipy.sandbox.multigrid.adaptive import augment_candidates
 
+from scipy.sandbox.multigrid.gallery import *
+from scipy.sandbox.multigrid.adaptive import *
 
 
 class TestAdaptiveSA(TestCase):
     def setUp(self):
-        pass
+        from numpy.random import seed
+        seed(0)
 
+    def test_poisson(self):
+        A = poisson( (100,100), format='csr' )
+
+        asa = adaptive_sa_solver(A, max_candidates = 1)
+        sa  = smoothed_aggregation_solver(A, B = ones((A.shape[0],1)) )
+
+        b = rand(A.shape[0])
+
+        sol0,residuals0 = asa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+        sol1,residuals1 =  sa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+       
+        conv_asa = (residuals0[-1]/residuals0[0])**(1.0/len(residuals0))
+        conv_sa  = (residuals1[-1]/residuals1[0])**(1.0/len(residuals1))
+        
+        assert( conv_asa < 1.1 * conv_sa ) #aSA shouldn't be any worse than SA
+
+#    def test_elasticity(self):
+#        A,B = linear_elasticity( (100,100), format='bsr' )
+#
+#        asa = adaptive_sa_solver(A, max_candidates = 3)
+#        sa  = smoothed_aggregation_solver(A, B=B )
+#
+#        b = rand(A.shape[0])
+#
+#        sol0,residuals0 = asa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+#        sol1,residuals1 =  sa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+#       
+#        conv_asa = (residuals0[-1]/residuals0[0])**(1.0/len(residuals0))
+#        conv_sa  = (residuals1[-1]/residuals1[0])**(1.0/len(residuals1))
+#       
+#        print "ASA convergence",conv_asa
+#        assert( conv_asa < 1.1 * conv_sa ) #aSA shouldn't be any worse than SA
+        
 #class TestAugmentCandidates(TestCase):
 #    def setUp(self):
 #        self.cases = []

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/utils.py	2008-01-19 19:43:04 UTC (rev 3852)
@@ -46,6 +46,8 @@
     if A.shape[0] != A.shape[1]:
         raise ValueError,'expected square matrix'
 
+    maxiter = min(A.shape[0],maxiter)
+
     #TODO make method adaptive
 
     numpy.random.seed(0)  #make results deterministic



More information about the Scipy-svn mailing list