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

scipy-svn@scip... scipy-svn@scip...
Sun Oct 21 20:18:52 CDT 2007


Author: wnbell
Date: 2007-10-21 20:18:47 -0500 (Sun, 21 Oct 2007)
New Revision: 3452

Modified:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h
   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:
made relaxation method ignore rows with zero diagonals
changed SA candidates to use dense 2D array


Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-22 01:18:47 UTC (rev 3452)
@@ -66,13 +66,13 @@
  
 
 
-def sa_hierarchy(A,Ws,x):
+def sa_hierarchy(A,Ws,B):
     """
     Construct multilevel hierarchy using Smoothed Aggregation
         Inputs:
           A  - matrix
           Is - list of constant prolongators
-          x  - "candidate" basis function to be approximated
+          B  - "candidate" basis function to be approximated
         Ouputs:
           (As,Is,Ps) - tuple of lists
                   - As - [A, Ps[0].T*A*Ps[0], Ps[1].T*A*Ps[1], ... ]
@@ -84,7 +84,7 @@
     Ps = []
 
     for W in Ws:
-        P,x = sa_fit_candidates(W,x)
+        P,B = sa_fit_candidates(W,B)
         I   = smoothed_prolongator(P,A)  
         A   = I.T.tocsr() * A * I
         As.append(A)
@@ -109,19 +109,21 @@
             raise ValueError,'small matrices not handled yet'
        
         #first candidate 
-        x,AggOps = self.__initialization_stage(A, blocks=blocks,\
-                                               max_levels=max_levels, max_coarse=max_coarse,\
-                                               mu=mu, epsilon=epsilon) 
+        x,AggOps = self.__initialization_stage(A, blocks = blocks, \
+                                               max_levels = max_levels, \
+                                               max_coarse = max_coarse, \
+                                               mu = mu, epsilon = epsilon) 
         
         Ws = AggOps
 
         self.candidates = [x]
 
         #create SA using x here
-        As,Is,Ps = sa_hierarchy(A,Ws,self.candidates)
+        As,Is,Ps, = sa_hierarchy(A,AggOps,self.candidates)
 
         for i in range(max_candidates - 1):
-            x = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps,mu=mu)    
+            #x = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps,mu=mu)    
+            x = self.__develop_candidate(As,Is,Ps,AggOps,self.candidates,mu=mu)    
             
             self.candidates.append(x)
             
@@ -189,9 +191,7 @@
         return x,AggOps  #first candidate,aggregation
 
 
-
-
-    def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps,mu):
+    def __develop_candidate(self,As,Is,Ps,AggOps,candidates,mu):
         #scipy.random.seed(0)  #TEST
         x = scipy.rand(A.shape[0])
         b = zeros_like(x)
@@ -201,7 +201,9 @@
         x = solver.solve(b, x0=x, tol=1e-10, maxiter=mu)
     
         #TEST FOR CONVERGENCE HERE
- 
+
+        #TODO augment candiates each time, then call fit_candidates?
+
         A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
         
         temp_Is = []
@@ -229,6 +231,46 @@
             x = I*x
         
         return x
+
+
+##    def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps,mu):
+##        #scipy.random.seed(0)  #TEST
+##        x = scipy.rand(A.shape[0])
+##        b = zeros_like(x)
+##
+##        solver = multilevel_solver(As,Is)
+##
+##        x = solver.solve(b, x0=x, tol=1e-10, maxiter=mu)
+##    
+##        #TEST FOR CONVERGENCE HERE
+## 
+##        A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
+##        
+##        temp_Is = []
+##        for i in range(len(As) - 2):
+##            P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])    
+## 
+##            I_l_new = smoothed_prolongator(P_l_new,A_l)
+##            A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
+##            bridge = make_bridge(Is[i+1],A_m_new.shape[0]) 
+## 
+##            temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Is[i+2:] )
+## 
+##            for n in range(mu):
+##                x_m = temp_solver.solve(zeros_like(x_m), x0=x_m, tol=1e-8, maxiter=1)
+## 
+##            temp_Is.append(I_l_new)
+## 
+##            W_l = vstack_csr(Ws[i+1],W_m_new)  #prepare for next iteration
+##            A_l = A_m_new
+##            x_l = x_m
+##            P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
+## 
+##        x = x_l
+##        for I in reversed(temp_Is):
+##            x = I*x
+##        
+##        return x
            
 
     def __augment_cycle(self,A,As,Ps,Ws,AggOps,x):
@@ -286,15 +328,16 @@
 A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
 blocks = arange(A.shape[0]/2).repeat(2)
 
-asa = adaptive_sa_solver(A,max_candidates=4,mu=12)
+asa = adaptive_sa_solver(A,max_candidates=4,mu=10)
 scipy.random.seed(0)  #make tests repeatable
 x = rand(A.shape[0])
-b = A*rand(A.shape[0])
+#b = A*rand(A.shape[0])
+b = zeros(A.shape[0])
 
 
 print "solving"
 if True:
-    x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=25,tol=1e-8,return_residuals=True)
+    x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=25,tol=1e-7,return_residuals=True)
 else:
     residuals = []
     def add_resid(x):

Modified: trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/multigridtools/relaxation.h	2007-10-22 01:18:47 UTC (rev 3452)
@@ -29,9 +29,10 @@
                 rsum += Ax[jj]*x[j];
         }
         
-        assert(diag != 0);
-        
-        x[i] = (b[i] - rsum)/diag;
+        //TODO raise error? inform user?
+        if (diag != 0){
+            x[i] = (b[i] - rsum)/diag;
+        }
     }
 }
 
@@ -64,9 +65,10 @@
                 rsum += Ax[jj]*temp[j];
         }
         
-        assert(diag != 0);
-        
-        x[i] = (1 - omega) * temp[i] + omega * ((b[i] - rsum)/diag);
+        //TODO raise error? inform user?
+        if (diag != 0){ 
+            x[i] = (1 - omega) * temp[i] + omega * ((b[i] - rsum)/diag);
+        }
     }
 }
 

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-10-22 01:18:47 UTC (rev 3452)
@@ -61,20 +61,53 @@
     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):
-    """
-    Create a multilevel solver using Smoothed Aggregation (SA)
+    """Create a multilevel solver using Smoothed Aggregation (SA)
 
-        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
+    *Parameters*:
+        
+        A : {csr_matrix}
+            NxN matrix in CSR 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))
+        blocks : {None, array_like} : optional
+            Array of length N that groups the variables into 'superblocks'.
+            For example, in a 2d vector-valued problem where the even 
+            variables [0,2,4,...N-2] correspond to the x-components and the 
+            odd variables [1,3,5,...,N-1] correspond to the y-components then
+            blocks=[0,0,1,1,2,2,...,N/2,N/2] is expected.  The default 
+            value blocks=None is equivalent to blocks=[0,1,2,..,N] which 
+            implies that each variable should be aggregated seperately.
+            The default is appropriate for scalar valued problems.
+        aggregation: {None, list of csr_matrix} : optional
+            List of csr_matrix objects that describe a user-defined 
+            multilevel aggregation of the variables.
+            TODO ELABORATE
+        max_levels: {integer} : optional
+            Maximum number of levels to be used in the multilevel solver.
+        max_coarse: {integer} : optional
+            Maximum number of variables permitted on the coarse grid.
+        epsilon: {float} : optional
+            Strength of connection parameter used in aggregation.
+        omega: {float} : optional
+            Damping parameter used in prolongator smoothing (0 < omega < 2)
+
+    *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
     
     """
     As = [A]
     Ps = []
     
     if candidates is None:
-        candidates = [ ones(A.shape[0]) ] # use constant vector
+        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:
@@ -203,10 +236,10 @@
     
     A = io.mmread('tests/sample_data/elas30_A.mtx').tocsr()
     candidates = io.mmread('tests/sample_data/elas30_nullspace.mtx')
-    candidates = [ array(candidates[:,x]) for x in range(candidates.shape[1]) ]
+    #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,epsilon=0,max_coarse=10,max_levels=10)
+    ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=100,max_levels=2)
     #ml = ruge_stuben_solver(A)
 
     x = rand(A.shape[0])

Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/sa.py	2007-10-22 01:18:47 UTC (rev 3452)
@@ -11,26 +11,7 @@
            'sa_interpolation','sa_smoothed_prolongator','sa_fit_candidates']
 
 
-##    nnz = A.nnz
-##
-##    indptr  = A.indptr
-##    indices = A.indices
-##    data    = A.data
-##
-##    if n != 1:
-##        # expand horizontally
-##        indptr = n*A.indptr
-##        indices = (n*indices).repeat(n) + tile(arange(n),nnz) 
-##
-##    if m != 1:
-##        #expand vertically
-##        indptr  = concatenate( (array([0]), cumsum(diff(A).repeat(m))) )
-##        indices = indices.repeat(m)
-##
-##    if m != 1 or n != 1:
-##        data = A.data.repeat(m*n)
 
-
 def sa_filtered_matrix(A,epsilon,blocks=None):
     """The filtered matrix is obtained from A by lumping all weak off-diagonal 
     entries onto the diagonal.  Weak off-diagonals are determined by
@@ -48,7 +29,8 @@
             Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
             A_filtered = csr_matrix((Sx,Sj,Sp),A.shape)
         else:
-            A_filtered = A  #TODO subtract weak blocks from diagonal blocks?
+            raise NotImplementedError,'blocks not handled yet'
+##            #TODO subtract weak blocks from diagonal blocks?
 ##            num_dofs   = A.shape[0]
 ##            num_blocks = blocks.max() + 1
 ##            
@@ -99,6 +81,7 @@
        
         B  = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
         #1-norms of blocks entries of A
+        #TODO figure out what to do for blocks here
         Block_A = B.T.tocsr() * csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape) * B 
 
         S = sa_strong_connections(Block_A,epsilon)
@@ -119,22 +102,22 @@
 
 
 def sa_fit_candidates(AggOp,candidates):
-    K = len(candidates)
+    
+    K = candidates.shape[1] # num candidates
 
     N_fine,N_coarse = AggOp.shape
 
-    if K > 1 and len(candidates[0]) == K*N_fine:
+    if K > 1 and candidates.shape[0] == K*N_fine:
         #see if fine space has been expanded (all levels except for first)
-        #TODO fix this and add unittests
-        #AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
         AggOp = expand_into_blocks(AggOp,K,1).tocsr()
         N_fine = K*N_fine
 
-    #TODO convert this to list of coarse candidates
     R = zeros((K*N_coarse,K)) #storage for coarse candidates
 
     candidate_matrices = []
-    for i,c in enumerate(candidates):
+
+    for i in range(K):
+        c = candidates[:,i]
         c = c[diff(AggOp.indptr) == 1]  #eliminate DOFs that aggregation misses
         X = csr_matrix((c,AggOp.indices,AggOp.indptr),dims=AggOp.shape)
        
@@ -154,6 +137,7 @@
 
         candidate_matrices.append(X)
 
+    # expand AggOp blocks horizontally 
     Q_indptr  = K*AggOp.indptr
     Q_indices = (K*AggOp.indices).repeat(K)
     for i in range(K):
@@ -163,9 +147,7 @@
         Q_data[i::K] = X.data
     Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
 
-    coarse_candidates = [ascontiguousarray(R[:,i]) for i in range(K)]
-
-    return Q,coarse_candidates
+    return Q,R
     
 def sa_smoothed_prolongator(A,T,epsilon,omega,blocks=None):
     """For a given matrix A and tentative prolongator T return the 
@@ -216,7 +198,7 @@
     P = sa_smoothed_prolongator(A,T,epsilon,omega,blocks)
 
     if blocks is not None:
-        blocks = arange(AggOp.shape[1]).repeat(len(candidates))
+        blocks = arange(AggOp.shape[1]).repeat(candidates.shape[1])
           
     return P,coarse_candidates,blocks
 

Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py	2007-10-22 01:18:47 UTC (rev 3452)
@@ -1,7 +1,7 @@
 from numpy.testing import *
 
 from numpy import sqrt,empty,ones,arange,array_split,eye,array, \
-                  zeros,diag,zeros_like,diff,matrix
+                  zeros,diag,zeros_like,diff,matrix,hstack,vstack
 from numpy.linalg import norm                  
 from scipy import rand
 from scipy.sparse import spdiags,csr_matrix,lil_matrix, \
@@ -72,7 +72,7 @@
         # 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)]
+        candidates = ones((6,1))
 
         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))
@@ -104,13 +104,13 @@
         #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)]
+        candidates = ones((6,1))
         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)]
+        candidates = ones((6,1))
         user_cases.append((A,AggOp,candidates))
 
         for A,AggOp,candidates in user_cases:
@@ -129,29 +129,28 @@
 
         ## tests where AggOp includes all DOFs
         #one candidate
-        self.cases.append((csr_matrix((ones(5),array([0,0,0,1,1]),arange(6)),dims=(5,2)),[ones(5)]))
-        self.cases.append((csr_matrix((ones(5),array([1,1,0,0,0]),arange(6)),dims=(5,2)),[ones(5)]))
-        self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)),[ones(9)]))
-        self.cases.append((csr_matrix((ones(9),array([2,1,0,0,1,2,1,0,2]),arange(10)),dims=(9,3)),[arange(9)]))
+        self.cases.append((csr_matrix((ones(5),array([0,0,0,1,1]),arange(6)),dims=(5,2)), ones((5,1)) ))
+        self.cases.append((csr_matrix((ones(5),array([1,1,0,0,0]),arange(6)),dims=(5,2)), ones((5,1)) ))
+        self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)), ones((9,1)) ))
+        self.cases.append((csr_matrix((ones(9),array([2,1,0,0,1,2,1,0,2]),arange(10)),dims=(9,3)), arange(9).reshape(9,1) ))
         #two candidates
-        self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),arange(5)),dims=(4,2)),[ones(4),arange(4)]))
-        self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)),[ones(9),arange(9)]))
-        self.cases.append((csr_matrix((ones(9),array([0,0,1,1,2,2,3,3,3]),arange(10)),dims=(9,4)),[ones(9),arange(9)]))
+        self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),arange(5)),dims=(4,2)), vstack((ones(4),arange(4))).T ))
+        self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)), vstack((ones(9),arange(9))).T ))
+        self.cases.append((csr_matrix((ones(9),array([0,0,1,1,2,2,3,3,3]),arange(10)),dims=(9,4)), vstack((ones(9),arange(9))).T ))
         #block candidates
-        self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)),[array([1]*9 + [0]*9),arange(2*9)]))
+        self.cases.append((csr_matrix((ones(9),array([0,0,0,1,1,1,2,2,2]),arange(10)),dims=(9,3)), vstack((array([1]*9 + [0]*9),arange(2*9))).T ))
         
         ## tests where AggOp excludes some DOFs
-        self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)),[ones(5)]))
-        self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)),[ones(5),arange(5)]))
-        self.cases.append((csr_matrix((ones(6),array([1,3,0,2,1,0]),array([0,0,1,2,2,3,4,5,5,6])),dims=(9,4)),[ones(9),arange(9)]))
+        self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), ones((5,1)) ))
+        self.cases.append((csr_matrix((ones(4),array([0,0,1,1]),array([0,1,2,2,3,4])),dims=(5,2)), vstack((ones(5),arange(5))).T ))
+        self.cases.append((csr_matrix((ones(6),array([1,3,0,2,1,0]),array([0,0,1,2,2,3,4,5,5,6])),dims=(9,4)), vstack((ones(9),arange(9))).T ))
 
     def check_all_cases(self):
         """Test case where aggregation includes all fine nodes"""
         
         def mask_candidate(AggOp,candidates):
             #mask out all DOFs that are not included in the aggregation
-            for c in candidates:
-                c[diff(AggOp.indptr) == 0] = 0
+            candidates[diff(AggOp.indptr) == 0,:] = 0
 
         for AggOp,fine_candidates in self.cases:
 
@@ -159,32 +158,25 @@
 
             Q,coarse_candidates = sa_fit_candidates(AggOp,fine_candidates)
 
-            assert_equal(len(coarse_candidates),len(fine_candidates))
-
             #each fine level candidate should be fit exactly
-            for fine,coarse in zip(fine_candidates,coarse_candidates):
-                assert_almost_equal(fine,Q*coarse)
-                assert_almost_equal(Q*(Q.T*fine),fine)
+            assert_almost_equal(fine_candidates,Q*coarse_candidates)
+            assert_almost_equal(Q*(Q.T*fine_candidates),fine_candidates)
 
-            #aggregate one more level (to a single aggregate)
-            K = len(coarse_candidates)
-            N = K*AggOp.shape[1]
-            AggOp = csr_matrix((ones(N),zeros(N),arange(N + 1)),dims=(N,1)) #aggregate to a single point
-            fine_candidates = coarse_candidates
-            
-            #mask_candidate(AggOp,fine_candidates)  #not needed
-            
-            #now check the coarser problem
-            Q,coarse_candidates = sa_fit_candidates(AggOp,fine_candidates)
+##            #aggregate one more level (to a single aggregate)
+##            K = coarse_candidates.shape[1]
+##            N = K*AggOp.shape[1]
+##            AggOp = csr_matrix((ones(N),zeros(N),arange(N + 1)),dims=(N,1)) #aggregate to a single point
+##            fine_candidates = coarse_candidates
+##            
+##            #mask_candidate(AggOp,fine_candidates)  #not needed
+##            
+##            #now check the coarser problem
+##            Q,coarse_candidates = sa_fit_candidates(AggOp,fine_candidates)
+##
+##            assert_almost_equal(fine_candidates,Q*coarse_candidates)
+##            assert_almost_equal(Q*(Q.T*fine_candidates),fine_candidates)
 
-            assert_equal(len(coarse_candidates),len(fine_candidates))
 
-            for fine,coarse in zip(fine_candidates,coarse_candidates):
-                assert_almost_equal(fine,Q*coarse)
-                assert_almost_equal(Q*(Q.T*fine),fine)
-
-
-
 class TestSASolverPerformance(NumpyTestCase):
     def setUp(self):
         self.cases = []
@@ -223,9 +215,9 @@
             DAD   = D*A*D
            
             if candidates is None:
-                candidates = [ ones(A.shape[0]) ]
+                candidates = ones((A.shape[0],1))
            
-            DAD_candidates = [ (D_inv * c) for c in candidates ]
+            DAD_candidates = D_inv * candidates
            
             #TODO force 2 level method and check that result is the same
 

Modified: trunk/scipy/sandbox/multigrid/tests/test_utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_utils.py	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/tests/test_utils.py	2007-10-22 01:18:47 UTC (rev 3452)
@@ -2,12 +2,13 @@
 
 import numpy
 import scipy
-from scipy import matrix,array,diag
+from scipy import matrix,array,diag,zeros
 from scipy.sparse import csr_matrix
 
 
 set_package_path()
-from scipy.sandbox.multigrid.utils import infinity_norm,diag_sparse
+from scipy.sandbox.multigrid.utils import infinity_norm, diag_sparse, \
+                                          expand_into_blocks
 restore_path()
 
 
@@ -52,9 +53,29 @@
         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_expand_into_blocks(self):
+        cases = []
+        cases.append( ( matrix([[1]]), (1,2) ) )
+        cases.append( ( matrix([[1]]), (2,1) ) )
+        cases.append( ( matrix([[1]]), (2,2) ) )
+        cases.append( ( matrix([[1,2]]), (1,2) ) )
+        cases.append( ( matrix([[1,2],[3,4]]), (2,2) ) )
+        cases.append( ( matrix([[0,0],[0,0]]), (3,1) ) )
+        cases.append( ( matrix([[0,1,0],[0,2,3]]), (3,2) ) )
+        cases.append( ( matrix([[1,0,0],[2,0,3]]), (2,5) ) )
 
+        for A,dims in cases:
+            m,n = dims
+            result = expand_into_blocks(csr_matrix(A),m,n).todense()
 
+            expected = zeros((m*A.shape[0],n*A.shape[1]))
+            for i in range(m):
+                for j in range(n):
+                    expected[i::m,j::n] = A
 
+            assert_equal(expected,result)
+
+
 if __name__ == '__main__':
     NumpyTest().run()
 

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-10-22 00:21:14 UTC (rev 3451)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-10-22 01:18:47 UTC (rev 3452)
@@ -3,22 +3,21 @@
 
 import numpy
 import scipy
-from numpy import ravel,arange,concatenate,tile
+from numpy import ravel,arange,concatenate,tile,asarray
 from scipy.linalg import norm
 from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
                         csr_matrix,csc_matrix,extract_diagonal, \
                         coo_matrix
 
 
-def approximate_spectral_radius(A,tol=0.1,maxiter=20):
+def approximate_spectral_radius(A,tol=0.1,maxiter=None):
     """
     Approximate the spectral radius of a symmetric matrix using ARPACK
     """
     from scipy.sandbox.arpack import eigen
-    return norm(eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0])
+    return norm(eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False))
 
 
-
 def infinity_norm(A):
     """
     Infinity norm of a sparse matrix (maximum absolute row sum).  This serves 
@@ -44,12 +43,16 @@
     if isspmatrix(A):
         return extract_diagonal(A)
     else:
-        return csr_matrix((A,arange(len(A)),arange(len(A)+1)),(len(A),len(A)))
+        return csr_matrix((asarray(A),arange(len(A)),arange(len(A)+1)),(len(A),len(A)))
 
 
 def hstack_csr(A,B):
-    #TODO OPTIMIZE THIS
-    assert(A.shape[0] == B.shape[0])
+    if not isspmatrix(A) or not isspmatrix(B):
+        raise TypeError,'expected sparse matrix'
+
+    if A.shape[0] != B.shape[0]:
+        raise ValueError,'row dimensions must agree'
+
     A = A.tocoo()
     B = B.tocoo()
     I = concatenate((A.row,B.row))
@@ -59,7 +62,12 @@
 
 def vstack_csr(A,B):
     #TODO OPTIMIZE THIS
-    assert(A.shape[1] == B.shape[1])
+    if not isspmatrix(A) or not isspmatrix(B):
+        raise TypeError,'expected sparse matrix'
+    
+    if A.shape[0] != B.shape[0]:
+        raise ValueError,'row dimensions must agree'
+
     A = A.tocoo()
     B = B.tocoo()
     I = concatenate((A.row,B.row+A.shape[0]))
@@ -84,6 +92,7 @@
               
     """
     #TODO EXPLAIN MORE
+    #TODO use spkron instead, time for compairson
 
     if n is None:
         n = m



More information about the Scipy-svn mailing list