[Scipy-svn] r3444 - in trunk/scipy/sandbox/multigrid: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Oct 17 16:29:17 EDT 2007
Author: wnbell
Date: 2007-10-17 15:29:14 -0500 (Wed, 17 Oct 2007)
New Revision: 3444
Modified:
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/tests/test_sa.py
Log:
added tests cases for user-defined aggregation
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-17 18:28:32 UTC (rev 3443)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-17 20:29:14 UTC (rev 3444)
@@ -77,9 +77,8 @@
candidates = [ ones(A.shape[0]) ] # use constant vector
if aggregation is None:
- while len(As) < max_levels and A.shape[0] > max_coarse:
- P,candidates = sa_interpolation(A,candidates,epsilon*0.5**(len(As)-1),omega=omega,blocks=blocks)
- blocks = None #only used for 1st level
+ 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)
A = (P.T.tocsr() * A) * P #galerkin operator
@@ -207,12 +206,12 @@
candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
blocks = arange(A.shape[0]/2).repeat(2)
- ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,max_coarse=10,max_levels=10)
+ ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=10,max_levels=10)
#ml = ruge_stuben_solver(A)
x = rand(A.shape[0])
- b = zeros_like(x)
- #b = rand(A.shape[0])
+ #b = zeros_like(x)
+ b = A*rand(A.shape[0])
if True:
x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-12,return_residuals=True)
@@ -221,7 +220,7 @@
def add_resid(x):
residuals.append(linalg.norm(b - A*x))
A.psolve = ml.psolve
- x_sol = linalg.cg(A,b,x0=x,maxiter=12,tol=1e-100,callback=add_resid)[0]
+ x_sol = linalg.cg(A,b,x0=x,maxiter=25,tol=1e-12,callback=add_resid)[0]
residuals = array(residuals)/residuals[0]
@@ -233,3 +232,5 @@
+
+
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2007-10-17 18:28:32 UTC (rev 3443)
+++ trunk/scipy/sandbox/multigrid/sa.py 2007-10-17 20:29:14 UTC (rev 3444)
@@ -8,7 +8,7 @@
import multigridtools
__all__ = ['sa_filtered_matrix','sa_strong_connections','sa_constant_interpolation',
- 'sa_interpolation','sa_fit_candidates']
+ 'sa_interpolation','sa_smoothed_prolongator','sa_fit_candidates']
def sa_filtered_matrix(A,epsilon,blocks=None):
@@ -83,6 +83,7 @@
Pj = Pj[blocks] #expand block aggregates into constituent dofs
Pp = B.indptr
Px = B.data
+
else:
S = sa_strong_connections(A,epsilon)
@@ -140,7 +141,18 @@
return Q,coarse_candidates
-
+def sa_smoothed_prolongator(A,T,epsilon,omega,blocks=None):
+ A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
+
+ D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
+ D_inv_A = D_inv * A_filtered
+ D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
+
+ # smooth tentative prolongator T
+ P = T - (D_inv_A*T)
+
+ return P
+
def sa_interpolation(A,candidates,epsilon=0.0,omega=4.0/3.0,blocks=None,AggOp=None):
if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
@@ -148,7 +160,7 @@
AggOp = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
else:
if not isspmatrix_csr(AggOp):
- raise TypeError,'aggregation is specified by a list of csr_matrix objects'
+ raise TypeError,'expected csr_matrix for argument AggOp'
if A.shape[1] != AggOp.shape[0]:
raise ValueError,'incompatible aggregation operator'
@@ -156,14 +168,19 @@
A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
- D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
- D_inv_A = D_inv * A_filtered
- D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
+ P = sa_smoothed_prolongator(A,T,epsilon,omega,blocks)
- # smooth tentative prolongator T
- P = T - (D_inv_A*T)
-
- return P,coarse_candidates
+## D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
+## D_inv_A = D_inv * A_filtered
+## D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
+##
+## # smooth tentative prolongator T
+## P = T - (D_inv_A*T)
+ if blocks is not None:
+ blocks = arange(AggOp.shape[1]).repeat(len(candidates))
+
+ return P,coarse_candidates,blocks
+
Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-17 18:28:32 UTC (rev 3443)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-17 20:29:14 UTC (rev 3444)
@@ -12,7 +12,8 @@
set_package_path()
import scipy.sandbox.multigrid
from scipy.sandbox.multigrid.sa import sa_strong_connections, sa_constant_interpolation, \
- sa_interpolation, sa_fit_candidates
+ sa_interpolation, sa_fit_candidates, \
+ sa_smoothed_prolongator
from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D, \
smoothed_aggregation_solver
from scipy.sandbox.multigrid.utils import diag_sparse
@@ -68,6 +69,15 @@
S_result = sa_constant_interpolation(A,epsilon,blocks=arange(A.shape[0]))
assert_array_equal(S_result.todense(),S_expected.todense())
+ # two aggregates in 1D
+ A = poisson_problem1D(6)
+ AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),dims=(6,2))
+ candidates = [ones(6)]
+
+ T_result,coarse_candidates_result = sa_fit_candidates(AggOp,candidates)
+ T_expected = csr_matrix((sqrt(1.0/3.0)*ones(6),array([0,0,0,1,1,1]),arange(7)),dims=(6,2))
+ assert_almost_equal(T_result.todense(),T_expected.todense())
+
#check simple block examples
A = csr_matrix(arange(16).reshape(4,4))
A = A + A.T
@@ -85,6 +95,32 @@
S_expected = matrix([[1,0],[1,0],[0,1],[0,1]])
assert_array_equal(S_result.todense(),S_expected)
+
+ def check_user_aggregation(self):
+ """check that the sa_interpolation accepts user-defined aggregates"""
+
+ user_cases = []
+
+ #simple 1d example w/ two aggregates
+ A = poisson_problem1D(6)
+ AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),dims=(6,2))
+ candidates = [ones(6)]
+ user_cases.append((A,AggOp,candidates))
+
+ #simple 1d example w/ two aggregates (not all nodes are aggregated)
+ A = poisson_problem1D(6)
+ AggOp = csr_matrix((ones(4),array([0,0,1,1]),array([0,1,1,2,3,3,4])),dims=(6,2))
+ candidates = [ones(6)]
+ user_cases.append((A,AggOp,candidates))
+
+ for A,AggOp,candidates in user_cases:
+ T,coarse_candidates_result = sa_fit_candidates(AggOp,candidates)
+
+ P_result = sa_interpolation(A,candidates,omega=4.0/3.0,AggOp=AggOp)[0]
+ P_expected = sa_smoothed_prolongator(A,T,epsilon=0.0,omega=4.0/3.0)
+
+ assert_almost_equal(P_result.todense(),P_expected.todense())
+
class TestFitCandidates(NumpyTestCase):
@@ -149,9 +185,7 @@
-
-
-class TestSASolver(NumpyTestCase):
+class TestSASolverPerformance(NumpyTestCase):
def setUp(self):
self.cases = []
More information about the Scipy-svn
mailing list