[Scipy-svn] r3452 - in trunk/scipy/sandbox/multigrid: . multigridtools tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Oct 21 21:18:52 EDT 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