[Scipy-svn] r3262 - trunk/scipy/sandbox/multigrid

[scipy-svn at scipy.org]

Author: wnbell
Date: 2007-08-25 18:15:34 -0500 (Sat, 25 Aug 2007)
New Revision: 3262

Modified:
   trunk/scipy/sandbox/multigrid/coarsen.py
   trunk/scipy/sandbox/multigrid/multilevel.py
Log:
change epsilon per level in SA
minor changes



Modified: trunk/scipy/sandbox/multigrid/coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/coarsen.py	2007-08-24 15:55:04 UTC (rev 3261)
+++ trunk/scipy/sandbox/multigrid/coarsen.py	2007-08-25 23:15:34 UTC (rev 3262)
@@ -35,10 +35,13 @@
     Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
     return scipy.sparse.csr_matrix((Sx,Sj,Sp),A.shape)
 
-def sa_constant_interpolation(A,epsilon=0.08):
+def sa_constant_interpolation(A,epsilon=None):
     if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
     
-    S = sa_strong_connections(A,epsilon)
+    if epsilon is not None:
+        S = sa_strong_connections(A,epsilon)
+    else:
+        S = A
     
     #tentative (non-smooth) interpolation operator I
     Ij = multigridtools.sa_get_aggregates(A.shape[0],S.indptr,S.indices)
@@ -48,7 +51,7 @@
     return scipy.sparse.csr_matrix((Ix,Ij,Ip))
 
 
-def sa_interpolation(A,epsilon=0.08,omega=4.0/3.0):
+def sa_interpolation(A,epsilon,omega=4.0/3.0):
     if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
     
     I = sa_constant_interpolation(A,epsilon)

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-08-24 15:55:04 UTC (rev 3261)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-08-25 23:15:34 UTC (rev 3262)
@@ -36,6 +36,15 @@
     return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocsr()
 
 def ruge_stuben_solver(A,max_levels=10,max_coarse=500):
+    """
+    Create a multilevel solver using Ruge-Stuben coarsening (Classical AMG)
+    
+        References:
+            "Multigrid"
+            Trottenberg, U., C. W. Oosterlee, and Anton Schuller. San Diego: Academic Press, 2001.
+            See Appendix A
+    
+    """
     As = [A]
     Ps = []
     
@@ -50,11 +59,20 @@
     return multilevel_solver(As,Ps)
 
 def smoothed_aggregation_solver(A,max_levels=10,max_coarse=500):
+    """
+    Create a multilevel solver using Smoothed Aggregation (SA)
+
+        References:
+            "Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems",
+                Petr Vanek and Jan Mandel and Marian Brezina
+                http://citeseer.ist.psu.edu/vanek96algebraic.html
+    
+    """
     As = [A]
     Ps = []
     
     while len(As) < max_levels  and A.shape[0] > max_coarse:
-        P = sa_interpolation(A)
+        P = sa_interpolation(A,epsilon=0.08*0.5**(len(As)-1))
         
         A = (P.T.tocsr() * A) * P     #galerkin operator
 
@@ -135,6 +153,7 @@
         coarse_b = self.Ps[lvl].T * residual
         
         if lvl == len(self.As) - 2:
+            #direct solver on coarsest level
             coarse_x[:] = scipy.linalg.solve(self.As[-1].todense(),coarse_b)
         else:   
             self.__solve(lvl+1,coarse_x,coarse_b)
@@ -155,15 +174,15 @@
 if __name__ == '__main__':
     from scipy import *
     A = poisson_problem2D(200)
-    asa = smoothed_aggregation_solver(A)
-    #asa = ruge_stuben_solver(A)
+    ml = smoothed_aggregation_solver(A)
+    #ml = ruge_stuben_solver(A)
     x = rand(A.shape[0])
     b = zeros_like(x)
     
     resid = []
     
     for n in range(10):
-        x = asa.solve(b,x,maxiter=1)
+        x = ml.solve(b,x,maxiter=1)
         resid.append(linalg.norm(A*x))