[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