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

scipy-svn@scip... scipy-svn@scip...
Tue Jan 15 10:15:17 CST 2008


Author: wnbell
Date: 2008-01-15 10:15:12 -0600 (Tue, 15 Jan 2008)
New Revision: 3838

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
Log:
updated adaptive SA to use BSR


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-15 16:15:12 UTC (rev 3838)
@@ -1,3 +1,6 @@
+__all__ = ['adaptive_sa_solver']
+
+
 import numpy,scipy,scipy.sparse
 
 from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
@@ -2,9 +5,10 @@
                   asarray, hstack, ascontiguousarray, isinf, dot
-from numpy.random import randn
+from numpy.random import randn, rand
 from scipy.sparse import csr_matrix,coo_matrix
 
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
-from sa import sa_constant_interpolation,sa_fit_candidates
 from utils import approximate_spectral_radius,hstack_csr,vstack_csr,diag_sparse
+from sa import sa_constant_interpolation, sa_smoothed_prolongator, \
+        sa_fit_candidates
 
@@ -19,13 +23,14 @@
     Construct multilevel hierarchy using Smoothed Aggregation
         Inputs:
           A  - matrix
-          Ps - list of constant prolongators
-          B  - "candidate" basis function to be approximated
+          B  - fine-level near nullspace candidates be approximated
+
         Ouputs:
-          (As,Ps,Ts) - tuple of lists
-                  - As - [A, Ts[0].T*A*Ts[0], Ts[1].T*A*Ts[1], ... ]
+          (As,Ps,Ts,Bs) - tuple of lists
+                  - As -  
                   - Ps - smoothed prolongators
                   - Ts - tentative prolongators
+                  - Bs - near nullspace candidates
     """
     As = [A]
     Ps = []
@@ -34,7 +39,7 @@
 
     for AggOp in AggOps:
         P,B = sa_fit_candidates(AggOp,B)
-        I   = sa_smoothed_prolongator(P,A)
+        I   = sa_smoothed_prolongator(A,P)
         A   = I.T.tocsr() * A * I
         As.append(A)
         Ts.append(P)
@@ -44,19 +49,17 @@
 
 
 class adaptive_sa_solver:
-    def __init__(self, A, blocks=None, aggregation=None, max_levels=10, max_coarse=100,\
+    def __init__(self, A, aggregation=None, max_levels=10, max_coarse=100,\
                  max_candidates=1, mu=5, epsilon=0.1):
-
-        self.A = A
-
+        
+        self.A  = A
         self.Rs = []
 
         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, \
+        x,AggOps = self.__initialization_stage(A, max_levels = max_levels, \
                                                max_coarse = max_coarse, \
                                                mu = mu, epsilon = epsilon, \
                                                aggregation = aggregation )
@@ -67,9 +70,6 @@
         for i in range(max_candidates - 1):
             x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
 
-            #TODO which is faster?
-            #As,Ps,Ts,Bs = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
-
             B = hstack((Bs[0],x))
             As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
 
@@ -89,7 +89,7 @@
         self.AggOps = AggOps
         self.Bs = Bs
 
-    def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon,aggregation):
+    def __initialization_stage(self,A,max_levels,max_coarse,mu,epsilon,aggregation):
         if aggregation is not None:
             max_coarse = 0
             max_levels = len(aggregation) + 1
@@ -100,31 +100,28 @@
 
         #step 1
         A_l = A
-        x   = randn(A_l.shape[0],1)
+        x   = rand(A_l.shape[0],1) # TODO see why randn() fails here
         skip_f_to_i = False
 
         #step 2
-        gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')
+        gauss_seidel(A_l, x, zeros_like(x), iterations=mu, sweep='symmetric')
 
         #step 3
         #TODO test convergence rate here
 
-        As = [A]
+        As     = [A]
+        Ps     = []
         AggOps = []
-        Ps = []
 
         while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
             if aggregation is None:
-                W_l   = sa_constant_interpolation(A_l,epsilon=0,blocks=blocks) #step 4b
+                W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b
             else:
-                W_l   = aggregation[len(AggOps)]
-            P_l,x = sa_fit_candidates(W_l,x)                             #step 4c
-            I_l   = sa_smoothed_prolongator(P_l,A_l)                          #step 4d
-            A_l   = I_l.T.tocsr() * A_l * I_l                              #step 4e
+                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
 
-            blocks = None #not needed on subsequent levels
-
-            print 'A_l.shape',A_l.shape
             AggOps.append(W_l)
             Ps.append(I_l)
             As.append(A_l)
@@ -132,7 +129,6 @@
             if A_l.shape <= max_coarse:  break
 
             if not skip_f_to_i:
-                print "."
                 x_hat = x.copy()                                                   #step 4g
                 gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')  #step 4h
                 x_A_x = dot(x.T,A_l*x)
@@ -176,7 +172,7 @@
             B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
             T,R = sa_fit_candidates(AggOp,B)
 
-            P = sa_smoothed_prolongator(T,A)
+            P = sa_smoothed_prolongator(A,T)
             A = P.T.tocsr() * A * P
 
             temp_Ps.append(P)
@@ -206,7 +202,7 @@
 #        for i in range(len(As) - 1):
 #            T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
 #
-#            P = sa_smoothed_prolongator(T,A)
+#            P = sa_smoothed_prolongator(A,T)
 #            A = P.T.tocsr() * A * P
 #
 #            new_As.append(A)
@@ -219,9 +215,6 @@
 #        return new_As,new_Ps,new_Ts,new_Bs
 
 
-
-
-
 #def augment_candidates(AggOp, old_Q, old_R, new_candidate):
 #    #TODO update P and A also
 #

Modified: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py	2008-01-15 16:15:12 UTC (rev 3838)
@@ -2,36 +2,32 @@
 from scipy import *
 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
 
+#A,B = linear_elasticity( (100,100) )
 
-#A = poisson( (200,200), spacing=(1,1e-2)
-#aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ]
+A = poisson( (200,200), format='csr' )
 
-#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 = 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
 
-#A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
 
-A,B = linear_elasticity( (100,100) )
 
 from time import clock; start = clock()
-asa = adaptive_sa_solver(A,max_candidates=3,mu=5,blocks=blocks,aggregation=aggregation)
+
+asa = adaptive_sa_solver(A, max_levels=2, max_candidates=1, mu=10)
+
 print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
 
-scipy.random.seed(0)  #make tests repeatable
+random.seed(0)  #make tests repeatable
 x = randn(A.shape[0])
-b = A*randn(A.shape[0])
-#b = zeros(A.shape[0])
+#b = A*randn(A.shape[0])
+b = zeros(A.shape[0])
 
 
 print "solving"
-if False:
+if True:
     x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
 else:
     residuals = []
@@ -74,10 +70,11 @@
 
 
 for c in asa.Bs[0].T:
-    #plot2d(c)
-    plot2d_arrows(c)
+    plot2d(c)
+    #plot2d_arrows(c)
     print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
 
+#plot2d(x_sol)
 
 ##W = asa.AggOps[0]*asa.AggOps[1]
 ##pcolor((W * rand(W.shape[1])).reshape((200,200)))

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/sa.py	2008-01-15 16:15:12 UTC (rev 3838)
@@ -97,6 +97,9 @@
         sa_constant_interpolation(csr_matrix(A),epsilon)
 
 def sa_fit_candidates(AggOp,candidates,tol=1e-10):
+    if not isspmatrix_csr(AggOp):
+        raise TypeError,'expected csr_matrix for argument AggOp'
+
     if candidates.dtype != 'float32':
         candidates = asarray(candidates,dtype='float64')
 
@@ -147,7 +150,7 @@
 
     return Q,R
 
-def sa_smoothed_prolongator(A,T,epsilon,omega):
+def sa_smoothed_prolongator(A,T,epsilon=0.0,omega=4.0/3.0):
     """For a given matrix A and tentative prolongator T return the
     smoothed prolongator P
 
@@ -161,7 +164,6 @@
         A_filtered - sa_filtered_matrix(A,epsilon)
     """
 
-
     A_filtered = sa_filtered_matrix(A,epsilon) #use filtered matrix for anisotropic problems
 
     # TODO use scale_rows()

Modified: trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_adaptive.py	2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/tests/test_adaptive.py	2008-01-15 16:15:12 UTC (rev 3838)
@@ -5,7 +5,7 @@
 
 
 from scipy.sandbox.multigrid.sa import sa_fit_candidates
-from scipy.sandbox.multigrid.adaptive import augment_candidates
+#from scipy.sandbox.multigrid.adaptive import augment_candidates
 
 
 



More information about the Scipy-svn mailing list