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

scipy-svn@scip... scipy-svn@scip...
Sun Nov 11 04:49:46 CST 2007


Author: wnbell
Date: 2007-11-11 04:48:59 -0600 (Sun, 11 Nov 2007)
New Revision: 3512

Modified:
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
   trunk/scipy/sandbox/multigrid/tests/test_utils.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
added support for diagonal scaling in SA solver


Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-11-09 20:53:50 UTC (rev 3511)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-11-11 10:48:59 UTC (rev 3512)
@@ -11,7 +11,7 @@
 from sa import sa_interpolation
 from rs import rs_interpolation
 from relaxation import gauss_seidel,jacobi,sor
-from utils import infinity_norm
+from utils import symmetric_rescaling, diag_sparse
 
 
 def poisson_problem1D(N):
@@ -43,8 +43,9 @@
 
         References:
             "Multigrid"
-                Trottenberg, U., C. W. Oosterlee, and Anton Schuller. San Diego: Academic Press, 2001.
-                See Appendix A
+                Trottenberg, U., C. W. Oosterlee, and Anton Schuller. 
+                San Diego: Academic Press, 2001.
+                Appendix A
 
     """
     As = [A]
@@ -59,8 +60,14 @@
         Ps.append(P)
 
     return multilevel_solver(As,Ps)
+        
 
-def smoothed_aggregation_solver(A,candidates=None,blocks=None,aggregation=None,max_levels=10,max_coarse=500,epsilon=0.08,omega=4.0/3.0):
+
+def smoothed_aggregation_solver(A, B=None, blocks=None, \
+                                aggregation=None, max_levels=10, \
+                                max_coarse=500, epsilon=0.0, \
+                                omega=4.0/3.0, symmetric=True, \
+                                rescale = True):
     """Create a multilevel solver using Smoothed Aggregation (SA)
 
     *Parameters*:
@@ -83,15 +90,20 @@
             List of csr_matrix objects that describe a user-defined
             multilevel aggregation of the variables.
             TODO ELABORATE
-        max_levels: {integer} : optional
+        max_levels: {integer} : default 10
             Maximum number of levels to be used in the multilevel solver.
-        max_coarse: {integer} : optional
+        max_coarse: {integer} : default 500
             Maximum number of variables permitted on the coarse grid.
-        epsilon: {float} : optional
+        epsilon: {float} : default 0.0
             Strength of connection parameter used in aggregation.
-        omega: {float} : optional
+        omega: {float} : default 4.0/3.0
             Damping parameter used in prolongator smoothing (0 < omega < 2)
-
+        symmetric: {boolean} : default True 
+            True if A is symmetric, False otherwise
+        rescale: {boolean} : default True
+            If True, symmetrically rescale A by the diagonal
+            i.e. A -> D * A * D,  where D is diag(A)^-0.5
+            
     *Example*:
         TODO
 
@@ -101,17 +113,30 @@
             http://citeseer.ist.psu.edu/vanek96algebraic.html
 
     """
+
+    if B is None:
+        B = ones((A.shape[0],1),dtype=A.dtype) # use constant vector
+    else:
+        B = asarray(B)
+
+    pre,post = None,None   #preprocess/postprocess
+    
+    if rescale:
+        D_sqrt,D_sqrt_inv,A = symmetric_rescaling(A)
+        D_sqrt,D_sqrt_inv = diag_sparse(D_sqrt),diag_sparse(D_sqrt_inv)
+        
+        B = D_sqrt * B  #scale candidates
+        def pre(x,b):
+            return D_sqrt*x,D_sqrt_inv*b
+        def post(x):
+            return D_sqrt_inv*x
+        
     As = [A]
     Ps = []
 
-    if candidates is None:
-        candidates = ones((A.shape[0],1),dtype=A.dtype) # use constant vector
-    else:
-        candiates = asarray(candidates)
-
     if aggregation is None:
         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)
+            P,B,blocks = sa_interpolation(A,B,epsilon*0.5**(len(As)-1),omega=omega,blocks=blocks)
 
             A = (P.T.tocsr() * A) * P     #galerkin operator
 
@@ -120,20 +145,22 @@
     else:
         #use user-defined aggregation
         for AggOp in aggregation:
-            P,candidates,blocks = sa_interpolation(A,candidates,omega=omega,AggOp=AggOp)
+            P,B,blocks = sa_interpolation(A,B,omega=omega,AggOp=AggOp)
 
             A = (P.T.tocsr() * A) * P     #galerkin operator
 
             As.append(A)
             Ps.append(P)
 
-    return multilevel_solver(As,Ps)
+    return multilevel_solver(As,Ps,preprocess=pre,postprocess=post)
 
 
 class multilevel_solver:
-    def __init__(self,As,Ps):
+    def __init__(self,As,Ps,preprocess=None,postprocess=None):
         self.As = As
         self.Ps = Ps
+        self.preprocess = preprocess
+        self.postprocess = postprocess
 
     def __repr__(self):
         output = 'multilevel_solver\n'
@@ -170,6 +197,8 @@
         else:
             x = array(x0) #copy
 
+        if self.preprocess is not None:
+            x,b = self.preprocess(x,b)
 
         #TODO change use of tol (relative tolerance) to agree with other iterative solvers
         A = self.As[0]
@@ -183,6 +212,9 @@
             if callback is not None:
                 callback(x)
 
+        if self.postprocess is not None:
+            x = self.postprocess(x)
+
         if return_residuals:
             return x,residuals
         else:

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-11-09 20:53:50 UTC (rev 3511)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-11-11 10:48:59 UTC (rev 3512)
@@ -4,7 +4,8 @@
                   ascontiguousarray
 from scipy.sparse import csr_matrix,isspmatrix_csr,spidentity
 
-from utils import diag_sparse,approximate_spectral_radius,expand_into_blocks
+from utils import diag_sparse, approximate_spectral_radius, \
+                  symmetric_rescaling, expand_into_blocks
 import multigridtools
 
 __all__ = ['sa_filtered_matrix','sa_strong_connections','sa_constant_interpolation',
@@ -103,6 +104,7 @@
 
 def sa_fit_candidates(AggOp,candidates,tol=1e-10):
     #TODO handle non-floating point candidates better
+    candidates = candidates.astype('float64')
 
     K = candidates.shape[1] # num candidates
 
@@ -181,9 +183,6 @@
     # smooth tentative prolongator T
     P = T - (D_inv_A*T)
 
-    #S = (spidentity(A.shape[0]).tocsr() - D_inv_A) #TODO drop this?
-    #P = S *(S * ( S * T))
-
     return P
 
 def sa_interpolation(A,candidates,epsilon=0.0,omega=4.0/3.0,blocks=None,AggOp=None):

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-11-09 20:53:50 UTC (rev 3511)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-11-11 10:48:59 UTC (rev 3512)
@@ -226,7 +226,7 @@
 
             #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)
+            ml = smoothed_aggregation_solver(DAD,DAD_candidates,max_coarse=100,max_levels=2,rescale=False)
 
             #print (D_inv*ml.Ps[0]).todense()
 

Modified: trunk/scipy/sandbox/multigrid/tests/test_utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_utils.py	2007-11-09 20:53:50 UTC (rev 3511)
+++ trunk/scipy/sandbox/multigrid/tests/test_utils.py	2007-11-11 10:48:59 UTC (rev 3512)
@@ -2,12 +2,13 @@
 
 import numpy
 import scipy
-from scipy import matrix,array,diag,zeros
+from numpy import matrix,array,diag,zeros,sqrt
 from scipy.sparse import csr_matrix
 
 
 set_package_path()
 from scipy.sandbox.multigrid.utils import infinity_norm, diag_sparse, \
+                                          symmetric_rescaling, \
                                           expand_into_blocks
 restore_path()
 
@@ -53,6 +54,33 @@
         A = matrix([[1.3,0,0],[0,5.5,0],[0,0,-2]])
         assert_equal(diag_sparse(array([1.3,5.5,-2])).todense(),csr_matrix(A).todense())
 
+
+    def check_symmetric_rescaling(self):
+        cases = []
+        cases.append( diag_sparse(array([1,2,3,4])) )
+        cases.append( diag_sparse(array([1,0,3,4])) )
+        
+        A = array([ [ 5.5,  3.5,  4.8],
+                    [ 2. ,  9.9,  0.5],
+                    [ 6.5,  2.6,  5.7]])
+        A = csr_matrix( A )
+        cases.append(A)
+        P = diag_sparse([1,0,1])
+        cases.append( P*A*P )
+        P = diag_sparse([0,1,0])
+        cases.append( P*A*P )
+        P = diag_sparse([1,-1,1])
+        cases.append( P*A*P )
+
+        for A in cases:
+            D_sqrt,D_sqrt_inv,DAD = symmetric_rescaling(A)
+
+            assert_almost_equal( diag_sparse(A) > 0, diag_sparse(DAD) )
+            assert_almost_equal( diag_sparse(DAD), D_sqrt*D_sqrt_inv )
+
+            D_sqrt,D_sqrt_inv = diag_sparse(D_sqrt),diag_sparse(D_sqrt_inv)
+            assert_almost_equal((D_sqrt_inv*A*D_sqrt_inv).todense(), DAD.todense())
+
     def check_expand_into_blocks(self):
         cases = []
         cases.append( ( matrix([[1]]), (1,2) ) )

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-11-09 20:53:50 UTC (rev 3511)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-11-11 10:48:59 UTC (rev 3512)
@@ -3,7 +3,7 @@
 
 import numpy
 import scipy
-from numpy import ravel,arange,concatenate,tile,asarray
+from numpy import ravel,arange,concatenate,tile,asarray,sqrt,diff
 from scipy.linalg import norm
 from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
                         csr_matrix,csc_matrix,extract_diagonal, \
@@ -26,7 +26,7 @@
 
     if isspmatrix_csr(A) or isspmatrix_csc(A):
         #avoid copying index and ptr arrays
-        abs_A = A.__class__((abs(A.data),A.indices,A.indptr),dims=A.shape,check=False)
+        abs_A = A.__class__((abs(A.data),A.indices,A.indptr),dims=A.shape)
         return (abs_A * numpy.ones(A.shape[1],dtype=A.dtype)).max()
     else:
         return (abs(A) * numpy.ones(A.shape[1],dtype=A.dtype)).max()
@@ -40,12 +40,37 @@
        - return a csr_matrix with A on the diagonal
     """
 
+    #TODO integrate into SciPy?
     if isspmatrix(A):
         return extract_diagonal(A)
     else:
         return csr_matrix((asarray(A),arange(len(A)),arange(len(A)+1)),(len(A),len(A)))
 
 
+def symmetric_rescaling(A):
+    if not (isspmatrix_csr(A) or isspmatrix_csc(A)):
+        raise TypeError,'expected csr_matrix or csc_matrix'
+
+    if A.shape[0] != A.shape[1]:
+        raise ValueError,'expected square matrix'
+
+    D = diag_sparse(A)
+    mask = D <= 0
+
+    D[mask] = 0
+    D_sqrt = sqrt(D)
+    D_sqrt_inv = 1.0/D_sqrt
+    D_sqrt_inv[mask] = 0
+    
+    #TODO time this against simple implementation
+    data = A.data * D_sqrt_inv[A.indices]
+    data *= D_sqrt_inv[arange(A.shape[0]).repeat(diff(A.indptr))]
+    
+    DAD = A.__class__((data,A.indices,A.indptr),dims=A.shape)
+
+    return D_sqrt,D_sqrt_inv,DAD
+
+
 def hstack_csr(A,B):
     if not isspmatrix(A) or not isspmatrix(B):
         raise TypeError,'expected sparse matrix'



More information about the Scipy-svn mailing list