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

scipy-svn@scip... scipy-svn@scip...
Thu Jan 24 15:11:20 CST 2008


Author: wnbell
Date: 2008-01-24 15:11:16 -0600 (Thu, 24 Jan 2008)
New Revision: 3859

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/sa.py
   trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
   trunk/scipy/sandbox/multigrid/tests/test_sa.py
Log:
reworked SA code
made substeps accept user-defined functions


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-23 22:37:53 UTC (rev 3858)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2008-01-24 21:11:16 UTC (rev 3859)
@@ -12,8 +12,8 @@
 from relaxation import gauss_seidel
 from multilevel import multilevel_solver
 from utils import approximate_spectral_radius,hstack_csr,vstack_csr,diag_sparse
-from sa import sa_constant_interpolation, sa_smoothed_prolongator, \
-        sa_fit_candidates
+from sa import sa_standard_aggregation, sa_smoothed_prolongator, \
+        sa_fit_candidates, sa_strong_connections
 
 
 
@@ -40,7 +40,7 @@
 
     for AggOp in AggOps:
         P,B = sa_fit_candidates(AggOp,B)
-        I   = sa_smoothed_prolongator(A,P)
+        I   = sa_smoothed_prolongator(A,A,P)
         A   = I.T.asformat(I.format) * A * I
         As.append(A)
         Ts.append(P)
@@ -159,12 +159,14 @@
 
     while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
         if aggregation is None:
-            W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b
+            C_l = sa_strong_connections(A_l,epsilon)
+            W_l = sa_standard_aggregation(C_l) #step 4b
         else:
             W_l = aggregation[len(AggOps)]
         P_l,x = sa_fit_candidates(W_l,x)                   #step 4c
-        I_l   = sa_smoothed_prolongator(A_l,P_l)           #step 4d
-        A_l   = I_l.T.asformat(I_l.format) * A_l * I_l                  #step 4e
+        I_l   = sa_smoothed_prolongator(A_l,A_l,P_l)       #step 4d
+        A_l   = I_l.T.asformat(I_l.format) * A_l * I_l     #step 4e
+        #TODO change variable names I_l -> P, P_l -> T
 
         AggOps.append(W_l)
         Ps.append(I_l)

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2008-01-23 22:37:53 UTC (rev 3858)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2008-01-24 21:11:16 UTC (rev 3859)
@@ -1,5 +1,4 @@
-__all__ = ['ruge_stuben_solver','smoothed_aggregation_solver',
-           'multilevel_solver']
+__all__ = ['ruge_stuben_solver','multilevel_solver']
 
 import scipy
 import numpy
@@ -8,7 +7,6 @@
 from scipy.linsolve import spsolve
 from scipy.sparse import dia_matrix
 
-from sa import sa_interpolation
 from rs import rs_interpolation
 from relaxation import gauss_seidel,jacobi,sor
 from utils import symmetric_rescaling, diag_sparse
@@ -41,94 +39,8 @@
 
 
 
-def smoothed_aggregation_solver(A, B=None, max_levels=10, max_coarse=500, \
-        epsilon=0.0, omega=4.0/3.0, symmetric=True, rescale=True, \
-        aggregation=None):
-    """Create a multilevel solver using Smoothed Aggregation (SA)
 
-    *Parameters*:
 
-        A : {csr_matrix}
-            NxN matrix in CSR or BSR format
-        B : {None, array_like} : optional
-            Near-nullspace candidates stored in the columns of an NxK array.
-            The default value B=None is equivalent to B=ones((N,1))
-        max_levels: {integer} : default 10
-            Maximum number of levels to be used in the multilevel solver.
-        max_coarse: {integer} : default 500
-            Maximum number of variables permitted on the coarse grid.
-        epsilon: {float} : default 0.0
-            Strength of connection parameter used in aggregation.
-        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
-        aggregation: {None, list of csr_matrix} : optional
-            List of csr_matrix objects that describe a user-defined
-            multilevel aggregation of the variables.
-            TODO ELABORATE
-
-    *Example*:
-        TODO
-
-    *References*:
-        "Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems",
-            Petr Vanek and Jan Mandel and Marian Brezina
-            http://citeseer.ist.psu.edu/vanek96algebraic.html
-
-    """
-
-    A = A.asfptype()
-
-    if B is None:
-        B = ones((A.shape[0],1),dtype=A.dtype) # use constant vector
-    else:
-        B = asarray(B,dtype=A.dtype)
-
-    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 = []
-    Rs = []
-
-    if aggregation is None:
-        while len(As) < max_levels and A.shape[0] > max_coarse:
-            P,B = sa_interpolation(A,B,epsilon*0.5**(len(As)-1),omega=omega)
-            R = P.T.asformat(P.format)
-
-            A = R * A * P     #galerkin operator
-
-            As.append(A)
-            Rs.append(R)
-            Ps.append(P)
-    else:
-        #use user-defined aggregation
-        for AggOp in aggregation:
-            P,B = sa_interpolation(A,B,omega=omega,AggOp=AggOp)
-            R = P.T.asformat(P.format)
-
-            A = R * A * P     #galerkin operator
-
-            As.append(A)
-            Rs.append(R)
-            Ps.append(P)
-
-    return multilevel_solver(As,Ps,Rs=Rs,preprocess=pre,postprocess=post)
-
-
 class multilevel_solver:
     def __init__(self,As,Ps,Rs=None,preprocess=None,postprocess=None):
         self.As = As
@@ -222,7 +134,7 @@
             #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0].reshape(coarse_x.shape)
             #A_inv = asarray(scipy.linalg.pinv2(self.As[-1].todense()))
             #coarse_x[:] = scipy.dot(A_inv,coarse_b)
-            print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
+            #print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
         else:
             self.__solve(lvl+1,coarse_x,coarse_b)
 

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2008-01-23 22:37:53 UTC (rev 3858)
+++ trunk/scipy/sandbox/multigrid/sa.py	2008-01-24 21:11:16 UTC (rev 3859)
@@ -1,15 +1,14 @@
-import scipy
-import numpy
-from numpy import array,arange,ones,zeros,sqrt,isinf,asarray,empty,diff,\
-                  ascontiguousarray
+from numpy import array, arange, ones, zeros, sqrt, asarray, empty, diff
 from scipy.sparse import csr_matrix, isspmatrix_csr, bsr_matrix, isspmatrix_bsr
 
+import multigridtools
+from multilevel import multilevel_solver
 from utils import diag_sparse, approximate_spectral_radius, \
                   symmetric_rescaling, scale_columns
-import multigridtools
 
-__all__ = ['sa_filtered_matrix','sa_strong_connections','sa_constant_interpolation',
-           'sa_interpolation','sa_smoothed_prolongator','sa_fit_candidates']
+__all__ = ['smoothed_aggregation_solver',
+        'sa_filtered_matrix','sa_strong_connections','sa_standard_aggregation',
+        'sa_smoothed_prolongator','sa_fit_candidates']
 
 
 
@@ -58,43 +57,54 @@
 
     return A_filtered
 
-def sa_strong_connections(A,epsilon):
-    if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+def sa_strong_connections(A,epsilon=0):
+    """Compute a strength of connection matrix C
 
-    Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+        C[i,j] = 1 if abs(A[i,j]) >= epsilon * abs(A[i,i] * A[j,j])
+        C[i,j] = 0 otherwise
 
-    return csr_matrix((Sx,Sj,Sp),A.shape)
-
-
-def sa_constant_interpolation(A,epsilon):
-    """Compute the sparsity pattern of the tentative prolongator
     """
-    if isspmatrix_csr(A): 
-        S = sa_strong_connections(A,epsilon)
 
-        Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
-        Pp = numpy.arange(len(Pj)+1)
-        Px = numpy.ones(len(Pj)) #TODO replace this with something else?
-        
-        return csr_matrix((Px,Pj,Pp))
+    if isspmatrix_csr(A):
+        if epsilon == 0:
+            return A
+        else:
+            fn = multigridtools.sa_strong_connections
+            Sp,Sj,Sx = fn(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+            return csr_matrix((Sx,Sj,Sp),A.shape)
 
     elif isspmatrix_bsr(A):
-
-        # the strength of connection matrix is based on the Frobenius norms of the blocks
         M,N = A.shape
         R,C = A.blocksize
 
         if R != C:
             raise ValueError,'matrix must have square blocks'
 
-        f_norms = (A.data*A.data).reshape(-1,R*C).sum(axis=1) #Frobenius norm of each block
+        if epsilon == 0:
+            data = ones( len(A.indices), dtype=A.dtype )
+            return csr_matrix((data,A.indices,A.indptr),shape=(M/R,N/C))
+        else:
+            # the strength of connection matrix is based on the 
+            # Frobenius norms of the blocks
+            data = (A.data*A.data).reshape(-1,R*C).sum(axis=1) 
+            A = csr_matrix((data,A.indices,A.indptr),shape=(M/R,N/C))
+            return sa_strong_connections(A,epsilon)
+    else:
+        raise TypeError('expected csr_matrix or bsr_matrix') 
+
+def sa_standard_aggregation(C):
+    """Compute the sparsity pattern of the tentative prolongator 
+    from a strength of connection matrix C
+    """
+    if isspmatrix_csr(C): 
         
-        A = csr_matrix((f_norms,A.indices,A.indptr),shape=(M/R,N/C))
-
-        return sa_constant_interpolation(A,epsilon)
-
+        Pj = multigridtools.sa_get_aggregates(C.shape[0],C.indptr,C.indices)
+        Pp = arange(len(Pj)+1)
+        Px = ones(len(Pj)) #TODO replace this with something else?
+        
+        return csr_matrix((Px,Pj,Pp))
     else:
-        sa_constant_interpolation(csr_matrix(A),epsilon)
+        raise TypeError('expected csr_matrix') 
 
 def sa_fit_candidates(AggOp,B,tol=1e-10):
     if not isspmatrix_csr(AggOp):
@@ -158,7 +168,7 @@
 
     return Q,R
 
-def sa_smoothed_prolongator(A,T,epsilon=0.0,omega=4.0/3.0):
+def sa_smoothed_prolongator(A,C,T,epsilon=0.0,omega=4.0/3.0):
     """For a given matrix A and tentative prolongator T return the
     smoothed prolongator P
 
@@ -187,21 +197,110 @@
 
     return P
 
-def sa_interpolation(A,B,epsilon=0.0,omega=4.0/3.0,AggOp=None):
 
+
+def smoothed_aggregation_solver(A, B=None, 
+        max_levels = 10, 
+        max_coarse = 500,
+        strength   = sa_strong_connections, 
+        aggregate  = sa_standard_aggregation,
+        tentative  = sa_fit_candidates,
+        smooth     = sa_smoothed_prolongator):
+    """Create a multilevel solver using Smoothed Aggregation (SA)
+
+    *Parameters*:
+
+        A : {csr_matrix, bsr_matrix}
+            Square matrix in CSR or BSR format
+        B : {None, array_like}
+            Near-nullspace candidates stored in the columns of an NxK array.
+            The default value B=None is equivalent to B=ones((N,1))
+        max_levels: {integer} : default 10
+            Maximum number of levels to be used in the multilevel solver.
+        max_coarse: {integer} : default 500
+            Maximum number of variables permitted on the coarse grid.
+        strength :
+            Function that computes the strength of connection matrix C
+                strength(A) -> C
+        aggregate : 
+            Function that computes an aggregation operator
+                aggregate(C) -> AggOp
+        tentative:
+            Function that computes a tentative prolongator
+                tentative(AggOp,B) -> T,B_coarse
+        smooth :
+            Function that smooths the tentative prolongator
+                smooth(A,C,T) -> P
+
+    Unused Parameters
+        epsilon: {float} : default 0.0
+            Strength of connection parameter used in aggregation.
+        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
+        aggregation: {None, list of csr_matrix} : optional
+            List of csr_matrix objects that describe a user-defined
+            multilevel aggregation of the variables.
+            TODO ELABORATE
+
+    *Example*:
+        TODO
+
+    *References*:
+        "Algebraic Multigrid by Smoothed Aggregation for Second and Fourth Order Elliptic Problems",
+            Petr Vanek and Jan Mandel and Marian Brezina
+            http://citeseer.ist.psu.edu/vanek96algebraic.html
+
+    """
+
+    A = A.asfptype()
+
     if not (isspmatrix_csr(A) or isspmatrix_bsr(A)):
-        A = csr_matrix(A)
+        raise TypeError('argument A must have type csr_matrix or bsr_matrix')
 
-    if AggOp is None:
-        AggOp = sa_constant_interpolation(A,epsilon=epsilon)
+    if A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    if B is None:
+        B = ones((A.shape[0],1),dtype=A.dtype) # use constant vector
     else:
-        if not isspmatrix_csr(AggOp):
-            AggOp = csr_matrix(AggOp)
-        if A.shape[1] != AggOp.shape[0]:
-            raise ValueError,'incompatible aggregation operator'
+        B = asarray(B,dtype=A.dtype)
 
+    pre,post = None,None   #preprocess/postprocess
 
-    T,coarse_candidates = sa_fit_candidates(AggOp,B)
-    P = sa_smoothed_prolongator(A,T,epsilon,omega)
-    return P,coarse_candidates
+    #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 = []
+    Rs = []
+
+    while len(As) < max_levels and A.shape[0] > max_coarse:
+        C     = strength(A)
+        AggOp = aggregate(C)
+        T,B   = tentative(AggOp,B)
+        P     = smooth(A,C,T)
+
+        R = P.T.asformat(P.format)
+
+        A = R * A * P     #galerkin operator
+
+        As.append(A)
+        Rs.append(R)
+        Ps.append(P)
+
+
+    return multilevel_solver(As,Ps,Rs=Rs,preprocess=pre,postprocess=post)
+
+

Modified: trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_adaptive.py	2008-01-23 22:37:53 UTC (rev 3858)
+++ trunk/scipy/sandbox/multigrid/tests/test_adaptive.py	2008-01-24 21:11:16 UTC (rev 3859)
@@ -5,8 +5,7 @@
 from scipy.sparse import csr_matrix
 
 
-from scipy.sandbox.multigrid.sa import sa_fit_candidates
-from scipy.sandbox.multigrid import smoothed_aggregation_solver
+from scipy.sandbox.multigrid.sa import sa_fit_candidates, smoothed_aggregation_solver
 #from scipy.sandbox.multigrid.adaptive import augment_candidates
 
 from scipy.sandbox.multigrid.gallery import *

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2008-01-23 22:37:53 UTC (rev 3858)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2008-01-24 21:11:16 UTC (rev 3859)
@@ -16,12 +16,8 @@
 
 
 import scipy.sandbox.multigrid
-from scipy.sandbox.multigrid.sa import sa_strong_connections, sa_constant_interpolation, \
-                                        sa_interpolation, sa_fit_candidates, \
-                                        sa_smoothed_prolongator
-from scipy.sandbox.multigrid.multilevel import smoothed_aggregation_solver
+from scipy.sandbox.multigrid.sa import *
 from scipy.sandbox.multigrid.utils import diag_sparse
-
 from scipy.sandbox.multigrid.gallery import poisson
 
 #def sparsity(A):
@@ -61,72 +57,71 @@
                 assert_almost_equal(S_result.todense(),S_expected.todense())
                 #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
 
-    def test_sa_constant_interpolation(self):
-        for A in self.cases:
-            for epsilon in [0.0,0.1,0.5,1.0]:
-                S_expected = reference_sa_constant_interpolation(A,epsilon)
+        ## two aggregates in 1D
+        #A = poisson( (6,), format='csr')
+        #AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),shape=(6,2))
+        #candidates = ones((6,1))
 
-                S_result   = sa_constant_interpolation(A,epsilon)
-                assert_array_equal(S_result.todense(),S_expected.todense())
+        #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)),shape=(6,2))
+        #assert_almost_equal(T_result.todense(),T_expected.todense())
 
-                #A = A.tobsr( blocksize=(1,1) )
-                #S_result   = sa_constant_interpolation(A,epsilon)
-                #assert_array_equal(S_result.todense(),S_expected.todense())
+        ##check simple block examples
+        #A = csr_matrix(arange(16).reshape(4,4))
+        #A = A + A.T
+        #A = A.tobsr(blocksize=(2,2))
 
-        # two aggregates in 1D
-        A = poisson( (6,), format='csr')
-        AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),shape=(6,2))
-        candidates = ones((6,1))
+        #S_result   = sa_standard_aggregation(A,epsilon=0.0)
+        #S_expected = matrix([1,1]).T
+        #assert_array_equal(S_result.todense(),S_expected)
 
-        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)),shape=(6,2))
-        assert_almost_equal(T_result.todense(),T_expected.todense())
+        #S_result   = sa_standard_aggregation(A,epsilon=0.5)
+        #S_expected = matrix([1,1]).T
+        #assert_array_equal(S_result.todense(),S_expected)
 
-        #check simple block examples
-        A = csr_matrix(arange(16).reshape(4,4))
-        A = A + A.T
-        A = A.tobsr(blocksize=(2,2))
+        #S_result   = sa_standard_aggregation(A,epsilon=2.0)
+        #S_expected = matrix([[1,0],[0,1]])
+        #assert_array_equal(S_result.todense(),S_expected)
 
-        S_result   = sa_constant_interpolation(A,epsilon=0.0)
-        S_expected = matrix([1,1]).T
-        assert_array_equal(S_result.todense(),S_expected)
+    def test_sa_standard_aggregation(self):
+        for C in self.cases:
+            S_expected = reference_sa_standard_aggregation(C)
 
-        S_result   = sa_constant_interpolation(A,epsilon=0.5)
-        S_expected = matrix([1,1]).T
-        assert_array_equal(S_result.todense(),S_expected)
+            S_result   = sa_standard_aggregation(C)
+            assert_array_equal(S_result.todense(),S_expected.todense())
 
-        S_result   = sa_constant_interpolation(A,epsilon=2.0)
-        S_expected = matrix([[1,0],[0,1]])
-        assert_array_equal(S_result.todense(),S_expected)
+            #A = A.tobsr( blocksize=(1,1) )
+            #S_result   = sa_constant_interpolation(A,epsilon)
+            #assert_array_equal(S_result.todense(),S_expected.todense())
 
 
-    def test_user_aggregation(self):
-        """check that the sa_interpolation accepts user-defined aggregates"""
+#    def test_user_aggregation(self):
+#        """check that the sa_interpolation accepts user-defined aggregates"""
+#
+#        user_cases = []
+#
+#        #simple 1d example w/ two aggregates
+#        A = poisson( (6,), format='csr')
+#        AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),shape=(6,2))
+#        candidates = ones((6,1))
+#        user_cases.append((A,AggOp,candidates))
+#
+#        #simple 1d example w/ two aggregates (not all nodes are aggregated)
+#        A = poisson( (6,), format='csr')
+#        AggOp = csr_matrix((ones(4),array([0,0,1,1]),array([0,1,1,2,3,3,4])),shape=(6,2))
+#        candidates = ones((6,1))
+#        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())
 
-        user_cases = []
 
-        #simple 1d example w/ two aggregates
-        A = poisson( (6,), format='csr')
-        AggOp = csr_matrix((ones(6),array([0,0,0,1,1,1]),arange(7)),shape=(6,2))
-        candidates = ones((6,1))
-        user_cases.append((A,AggOp,candidates))
 
-        #simple 1d example w/ two aggregates (not all nodes are aggregated)
-        A = poisson( (6,), format='csr')
-        AggOp = csr_matrix((ones(4),array([0,0,1,1]),array([0,1,1,2,3,3,4])),shape=(6,2))
-        candidates = ones((6,1))
-        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(TestCase):
     def setUp(self):
         self.cases = []
@@ -191,8 +186,8 @@
     def test_basic(self):
         """check that method converges at a reasonable rate"""
 
-        for A,candidates in self.cases:
-            ml = smoothed_aggregation_solver(A,candidates,max_coarse=10,max_levels=10)
+        for A,B in self.cases:
+            ml = smoothed_aggregation_solver(A,B,max_coarse=10,max_levels=10)
 
             numpy.random.seed(0) #make tests repeatable
 
@@ -221,7 +216,7 @@
         #TODO force 2 level method and check that result is the same
  
         #sa1 = smoothed_aggregation_solver(A, B, max_levels=2, rescale=False)
-        sa2 = smoothed_aggregation_solver(D*A*D, D_inv * B, max_levels=2, rescale=True)
+        sa2 = smoothed_aggregation_solver(D*A*D, D_inv * B, max_levels=2)
  
         #assert_almost_equal( sa2.Ps[0], sa1.Ps[0] 
         x_sol,residuals = sa2.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
@@ -251,16 +246,15 @@
 
 # note that this method only tests the current implementation, not
 # all possible implementations
-def reference_sa_constant_interpolation(A,epsilon):
-    S = sa_strong_connections(A,epsilon)
-    S = array_split(S.indices,S.indptr[1:-1])
+def reference_sa_standard_aggregation(C):
+    S = array_split(C.indices,C.indptr[1:-1])
 
-    n = A.shape[0]
+    n = C.shape[0]
 
     R = set(range(n))
     j = 0
 
-    aggregates = empty(n,dtype=A.indices.dtype)
+    aggregates = empty(n,dtype=C.indices.dtype)
     aggregates[:] = -1
 
     # Pass #1



More information about the Scipy-svn mailing list