# [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