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

scipy-svn@scip... scipy-svn@scip...
Wed Oct 17 15:29:17 CDT 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 @@
residuals.append(linalg.norm(b - A*x))
A.psolve = ml.psolve

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
+
+    # 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
+    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
+##
+##    # 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 = []

```