[Scipy-svn] r3395 - in trunk/scipy/sandbox/multigrid: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Oct 3 17:38:00 EDT 2007
Author: wnbell
Date: 2007-10-03 16:37:54 -0500 (Wed, 03 Oct 2007)
New Revision: 3395
Added:
trunk/scipy/sandbox/multigrid/tests/test_sa.py
Removed:
trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/coarsen.py
trunk/scipy/sandbox/multigrid/multilevel.py
Log:
separated SA and RS codes
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-03 21:37:54 UTC (rev 3395)
@@ -4,36 +4,36 @@
from relaxation import gauss_seidel
from multilevel import multilevel_solver
-from coarsen import sa_constant_interpolation
+from sa import sa_constant_interpolation
#from utils import infinity_norm
from utils import approximate_spectral_radius
-def fit_candidate(I,x):
- """
- For each aggregate in I (i.e. each column of I) compute vector R and
- sparse matrix Q (having the sparsity of I) such that the following holds:
+##def fit_candidate(I,x):
+## """
+## For each aggregate in I (i.e. each column of I) compute vector R and
+## sparse matrix Q (having the sparsity of I) such that the following holds:
+##
+## Q*R = x and Q^T*Q = I
+##
+## In otherwords, find a prolongator Q with orthonormal columns so that
+## x is represented exactly on the coarser level by R.
+## """
+## x = asarray(x)
+## Q = csr_matrix((x.copy(),I.indices,I.indptr),dims=I.shape,check=False)
+## R = sqrt(ravel(csr_matrix((x*x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0))) #column 2-norms
+##
+## Q.data *= (1.0/R)[Q.indices] #normalize columns of Q
+##
+## #print "norm(R)",scipy.linalg.norm(R)
+## #print "min(R),max(R)",min(R),max(R)
+## #print "infinity_norm(Q.T*Q - I) ",infinity_norm((Q.T.tocsr() * Q - scipy.sparse.spidentity(Q.shape[1])))
+## #print "norm(Q*R - x)",scipy.linalg.norm(Q*R - x)
+## #print "norm(x - Q*Q.Tx)",scipy.linalg.norm(x - Q*(Q.T*x))
+## return Q,R
- Q*R = x and Q^T*Q = I
- In otherwords, find a prolongator Q with orthonormal columns so that
- x is represented exactly on the coarser level by R.
- """
- x = asarray(x)
- Q = csr_matrix((x.copy(),I.indices,I.indptr),dims=I.shape,check=False)
- R = sqrt(ravel(csr_matrix((x*x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0))) #column 2-norms
- Q.data *= (1.0/R)[Q.indices] #normalize columns of Q
-
- #print "norm(R)",scipy.linalg.norm(R)
- #print "min(R),max(R)",min(R),max(R)
- #print "infinity_norm(Q.T*Q - I) ",infinity_norm((Q.T.tocsr() * Q - scipy.sparse.spidentity(Q.shape[1])))
- #print "norm(Q*R - x)",scipy.linalg.norm(Q*R - x)
- #print "norm(x - Q*Q.Tx)",scipy.linalg.norm(x - Q*(Q.T*x))
- return Q,R
-
-
-
##def orthonormalize_candidate(I,x,basis):
## Px = csr_matrix((x,I.indices,I.indptr),dims=I.shape,check=False)
## Rs = []
@@ -62,7 +62,6 @@
V = concatenate((A.data,B.data))
return coo_matrix((V,(I,J)),dims=(A.shape[0],A.shape[1]+B.shape[1])).tocsr()
-
def vstack_csr(A,B):
#OPTIMIZE THIS
assert(A.shape[1] == B.shape[1])
@@ -145,7 +144,6 @@
Ps = []
for W in Ws:
- #P,x = fit_candidate(W,x)
P,x = fit_candidates(W,x)
I = smoothed_prolongator(P,A)
A = I.T.tocsr() * A * I
Modified: trunk/scipy/sandbox/multigrid/coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/coarsen.py 2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/coarsen.py 2007-10-03 21:37:54 UTC (rev 3395)
@@ -1,163 +1,163 @@
-import scipy
-import numpy
-
-from numpy import arange,ones,zeros,sqrt,isinf,asarray,empty
-from scipy.sparse import csr_matrix,isspmatrix_csr
-
-from utils import diag_sparse,approximate_spectral_radius
-import multigridtools
-
-
-def rs_strong_connections(A,theta):
- """
- Return a strength of connection matrix using the method of Ruge and Stuben
-
- An off-diagonal entry A[i.j] is a strong connection iff
-
- -A[i,j] >= theta * max( -A[i,k] ) where k != i
- """
- if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
- Sp,Sj,Sx = multigridtools.rs_strong_connections(A.shape[0],theta,A.indptr,A.indices,A.data)
- return csr_matrix((Sx,Sj,Sp),A.shape)
-
-
-def rs_interpolation(A,theta=0.25):
- if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
- S = rs_strong_connections(A,theta)
-
- T = S.T.tocsr() #transpose S for efficient column access
-
- Ip,Ij,Ix = multigridtools.rs_interpolation(A.shape[0],\
- A.indptr,A.indices,A.data,\
- S.indptr,S.indices,S.data,\
- T.indptr,T.indices,T.data)
-
- return csr_matrix((Ix,Ij,Ip))
-
-
-def sa_strong_connections(A,epsilon):
- if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
- Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
-
- return csr_matrix((Sx,Sj,Sp),A.shape)
-
-def sa_constant_interpolation(A,epsilon,blocks=None):
- if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
- #handle epsilon = 0 case without creating strength of connection matrix?
-
- if blocks is not None:
- num_dofs = A.shape[0]
- num_blocks = blocks.max()
-
- if num_dofs != len(blocks):
- raise ValueError,'improper block specification'
-
- # for non-scalar problems, use pre-defined blocks in aggregation
- # the strength of connection matrix is based on the Frobenius norms of the blocks
-
- B = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
- Block_Frob = B.T.tocsr() * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B #Frobenius norms of blocks entries of A
-
- S = sa_strong_connections(Block_Frob,epsilon)
-
- Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
- Pj = Pj[blocks] #expand block aggregates into constituent dofs
- Pp = B.indptr
- Px = B.data
- else:
- 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))
-
- return csr_matrix((Px,Pj,Pp))
-
-
-def fit_candidates(AggOp,candidates):
- K = len(candidates)
-
- N_fine,N_coarse = AggOp.shape
-
- if K > 1 and len(candidates[0]) == K*N_fine:
- #see if fine space has been expanded (all levels except for first)
- AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
- N_fine = K*N_fine
-
- R = zeros((K*N_coarse,K))
-
- candidate_matrices = []
- for i,c in enumerate(candidates):
- X = csr_matrix((c.copy(),AggOp.indices,AggOp.indptr),dims=AggOp.shape)
-
- #TODO optimize this
-
- #orthogonalize X against previous
- for j,A in enumerate(candidate_matrices):
- D_AtX = csr_matrix((A.data*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
- R[j::K,i] = D_AtX
- X.data -= D_AtX[X.indices] * A.data
-
- #AtX = csr_matrix(A.T.tocsr() * X
- #R[j::K,i] = AtX.data
- #X = X - A * AtX
-
- #normalize X
- XtX = X.T.tocsr() * X
- col_norms = sqrt(asarray(XtX.sum(axis=0)).flatten())
- R[i::K,i] = col_norms
- col_norms = 1.0/col_norms
- col_norms[isinf(col_norms)] = 0
- X.data *= col_norms[X.indices]
-
- candidate_matrices.append(X)
-
-
- Q_indptr = K*AggOp.indptr
- Q_indices = (K*AggOp.indices).repeat(K)
- for i in range(K):
- Q_indices[i::K] += i
- Q_data = empty(N_fine * K)
- for i,X in enumerate(candidate_matrices):
- Q_data[i::K] = X.data
- Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
-
- coarse_candidates = [R[:,i] for i in range(K)]
-
- return Q,coarse_candidates
-
-## S = sa_strong_connections(A,epsilon)
+##import scipy
+##import numpy
##
-## #tentative (non-smooth) interpolation operator I
-## Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
-## Pp = numpy.arange(len(Pj)+1)
-## Px = numpy.ones(len(Pj))
+##from numpy import arange,ones,zeros,sqrt,isinf,asarray,empty
+##from scipy.sparse import csr_matrix,isspmatrix_csr
+##
+##from utils import diag_sparse,approximate_spectral_radius
+##import multigridtools
+##
+##
+##def rs_strong_connections(A,theta):
+## """
+## Return a strength of connection matrix using the method of Ruge and Stuben
+##
+## An off-diagonal entry A[i.j] is a strong connection iff
+##
+## -A[i,j] >= theta * max( -A[i,k] ) where k != i
+## """
+## if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##
+## Sp,Sj,Sx = multigridtools.rs_strong_connections(A.shape[0],theta,A.indptr,A.indices,A.data)
+## return csr_matrix((Sx,Sj,Sp),A.shape)
+##
+##
+##def rs_interpolation(A,theta=0.25):
+## if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
##
-## return scipy.sparse.csr_matrix((Px,Pj,Pp))
-
-##def sa_smoother(A,S,omega):
-## Bp,Bj,Bx = multigridtools.sa_smoother(A.shape[0],omega,A.indptr,A.indices,A.data,S.indptr,S.indices,S.data)
+## S = rs_strong_connections(A,theta)
##
-## return csr_matrix((Bx,Bj,Bp),dims=A.shape)
-
-def sa_interpolation(A,candidates,epsilon,omega=4.0/3.0,blocks=None):
- if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
-
- AggOp = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
- T,coarse_candidates = fit_candidates(AggOp,candidates)
-
- D_inv = diag_sparse(1.0/diag_sparse(A))
-
- D_inv_A = D_inv * A
- D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
-
- P = T - (D_inv_A*T) #same as I=S*P, (faster?)
-
- return P,coarse_candidates
-
-
-
+## T = S.T.tocsr() #transpose S for efficient column access
+##
+## Ip,Ij,Ix = multigridtools.rs_interpolation(A.shape[0],\
+## A.indptr,A.indices,A.data,\
+## S.indptr,S.indices,S.data,\
+## T.indptr,T.indices,T.data)
+##
+## return csr_matrix((Ix,Ij,Ip))
+##
+##
+##def sa_strong_connections(A,epsilon):
+## if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##
+## Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+##
+## return csr_matrix((Sx,Sj,Sp),A.shape)
+##
+##def sa_constant_interpolation(A,epsilon,blocks=None):
+## if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##
+## #handle epsilon = 0 case without creating strength of connection matrix?
+##
+## if blocks is not None:
+## num_dofs = A.shape[0]
+## num_blocks = blocks.max()
+##
+## if num_dofs != len(blocks):
+## raise ValueError,'improper block specification'
+##
+## # for non-scalar problems, use pre-defined blocks in aggregation
+## # the strength of connection matrix is based on the Frobenius norms of the blocks
+##
+## B = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
+## Block_Frob = B.T.tocsr() * csr_matrix((A.data**2,A.indices,A.indptr),dims=A.shape) * B #Frobenius norms of blocks entries of A
+##
+## S = sa_strong_connections(Block_Frob,epsilon)
+##
+## Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+## Pj = Pj[blocks] #expand block aggregates into constituent dofs
+## Pp = B.indptr
+## Px = B.data
+## else:
+## 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))
+##
+## return csr_matrix((Px,Pj,Pp))
+##
+##
+##def fit_candidates(AggOp,candidates):
+## K = len(candidates)
+##
+## N_fine,N_coarse = AggOp.shape
+##
+## if K > 1 and len(candidates[0]) == K*N_fine:
+## #see if fine space has been expanded (all levels except for first)
+## AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
+## N_fine = K*N_fine
+##
+## R = zeros((K*N_coarse,K))
+##
+## candidate_matrices = []
+## for i,c in enumerate(candidates):
+## X = csr_matrix((c.copy(),AggOp.indices,AggOp.indptr),dims=AggOp.shape)
+##
+## #TODO optimize this
+##
+## #orthogonalize X against previous
+## for j,A in enumerate(candidate_matrices):
+## D_AtX = csr_matrix((A.data*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
+## R[j::K,i] = D_AtX
+## X.data -= D_AtX[X.indices] * A.data
+##
+## #AtX = csr_matrix(A.T.tocsr() * X
+## #R[j::K,i] = AtX.data
+## #X = X - A * AtX
+##
+## #normalize X
+## XtX = X.T.tocsr() * X
+## col_norms = sqrt(asarray(XtX.sum(axis=0)).flatten())
+## R[i::K,i] = col_norms
+## col_norms = 1.0/col_norms
+## col_norms[isinf(col_norms)] = 0
+## X.data *= col_norms[X.indices]
+##
+## candidate_matrices.append(X)
+##
+##
+## Q_indptr = K*AggOp.indptr
+## Q_indices = (K*AggOp.indices).repeat(K)
+## for i in range(K):
+## Q_indices[i::K] += i
+## Q_data = empty(N_fine * K)
+## for i,X in enumerate(candidate_matrices):
+## Q_data[i::K] = X.data
+## Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(N_fine,K*N_coarse))
+##
+## coarse_candidates = [R[:,i] for i in range(K)]
+##
+## return Q,coarse_candidates
+##
+#### S = sa_strong_connections(A,epsilon)
+####
+#### #tentative (non-smooth) interpolation operator I
+#### Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+#### Pp = numpy.arange(len(Pj)+1)
+#### Px = numpy.ones(len(Pj))
+####
+#### return scipy.sparse.csr_matrix((Px,Pj,Pp))
+##
+####def sa_smoother(A,S,omega):
+#### Bp,Bj,Bx = multigridtools.sa_smoother(A.shape[0],omega,A.indptr,A.indices,A.data,S.indptr,S.indices,S.data)
+####
+#### return csr_matrix((Bx,Bj,Bp),dims=A.shape)
+##
+##def sa_interpolation(A,candidates,epsilon,omega=4.0/3.0,blocks=None):
+## if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
+##
+## AggOp = sa_constant_interpolation(A,epsilon=epsilon,blocks=blocks)
+## T,coarse_candidates = fit_candidates(AggOp,candidates)
+##
+## D_inv = diag_sparse(1.0/diag_sparse(A))
+##
+## D_inv_A = D_inv * A
+## D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
+##
+## P = T - (D_inv_A*T) #same as I=S*P, (faster?)
+##
+## return P,coarse_candidates
+##
+##
+##
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-03 21:37:54 UTC (rev 3395)
@@ -2,12 +2,13 @@
'ruge_stuben_solver','smoothed_aggregation_solver',
'multilevel_solver']
-from numpy.linalg import norm
-from numpy import zeros,zeros_like,array
import scipy
import numpy
+from numpy import zeros,zeros_like,array
+from numpy.linalg import norm
-from coarsen import sa_interpolation,rs_interpolation
+from sa import sa_interpolation
+from rs import rs_interpolation
from relaxation import gauss_seidel,jacobi,sor
from utils import infinity_norm
@@ -188,11 +189,12 @@
#A = io.mmread("9pt-100x100.mtx").tocsr()
#A = io.mmread("/home/nathan/Desktop/9pt/9pt-100x100.mtx").tocsr()
#A = io.mmread("/home/nathan/Desktop/BasisShift_W_EnergyMin_Luke/9pt-5x5.mtx").tocsr()
+
#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]) ]
- ml = smoothed_aggregation_solver(A,candidates,max_coarse=100,max_levels=3)
+ ml = smoothed_aggregation_solver(A,candidates,max_coarse=10,max_levels=10)
#ml = ruge_stuben_solver(A)
x = rand(A.shape[0])
Deleted: trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_coarsen.py 2007-10-03 18:42:23 UTC (rev 3394)
+++ trunk/scipy/sandbox/multigrid/tests/test_coarsen.py 2007-10-03 21:37:54 UTC (rev 3395)
@@ -1,161 +0,0 @@
-from numpy.testing import *
-
-from numpy import sqrt,empty,ones,arange,array_split
-from scipy import rand
-from scipy.sparse import spdiags,csr_matrix,lil_matrix
-import numpy
-
-set_package_path()
-import scipy.sandbox.multigrid
-from scipy.sandbox.multigrid.coarsen import sa_strong_connections,sa_constant_interpolation
-from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D
-restore_path()
-
-
-def reference_sa_strong_connections(A,epsilon):
- A_coo = A.tocoo()
- S = lil_matrix(A.shape)
-
- for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
- if i == j: continue #skip diagonal
-
- if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
- S[i,j] = v
-
- return S.tocsr()
-
-
-# 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])
-
- n = A.shape[0]
-
- R = set(range(n))
- j = 0
-
- aggregates = empty(n,dtype=A.indices.dtype)
- aggregates[:] = -1
-
- # Pass #1
- for i,row in enumerate(S):
- Ni = set(row) | set([i])
-
- if Ni.issubset(R):
- R -= Ni
- for x in Ni:
- aggregates[x] = j
- j += 1
-
- # Pass #2
- Old_R = R.copy()
- for i,row in enumerate(S):
- if i not in R: continue
-
- for x in row:
- if x not in Old_R:
- aggregates[i] = aggregates[x]
- R.remove(i)
- break
-
-
- # Pass #3
- for i,row in enumerate(S):
- if i not in R: continue
- Ni = set(row) | set([i])
-
- for x in Ni:
- if x in R:
- aggregates[x] = j
- j += 1
-
- assert(len(R) == 0)
-
- Pj = aggregates
- Pp = arange(n+1)
- Px = ones(n)
-
- return csr_matrix((Px,Pj,Pp))
-
-class TestSaStrongConnections(NumpyTestCase):
- def check_simple(self):
- N = 4
- A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
- S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
- assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense()) #all connections are strong
- assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
-
- N = 100
- A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
- S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
- assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense()) #all connections are strong
- assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
-
- def check_random(self):
- numpy.random.seed(0)
-
- for N in [2,3,5,10]:
- A = csr_matrix(rand(N,N))
- for epsilon in [0.0,0.1,0.5,0.8,1.0,10.0]:
- S_result = sa_strong_connections(A,epsilon)
- S_expected = reference_sa_strong_connections(A,epsilon)
- assert_array_equal(S_result.todense(),S_expected.todense())
-
- def check_poisson1D(self):
- for N in [2,3,5,7,10,11,19]:
- A = poisson_problem1D(N)
- for epsilon in [0.0,0.1,0.5,0.8,1.0]:
- S_result = sa_strong_connections(A,epsilon)
- S_expected = reference_sa_strong_connections(A,epsilon)
- assert_array_equal(S_result.todense(),S_expected.todense())
-
- def check_poisson2D(self):
- for N in [2,3,5,7,10,11,19]:
- A = poisson_problem2D(N)
- for epsilon in [0.0,0.1,0.5,0.8,1.0]:
- S_result = sa_strong_connections(A,epsilon)
- S_expected = reference_sa_strong_connections(A,epsilon)
- assert_array_equal(S_result.todense(),S_expected.todense())
-
-
-class TestSaConstantInterpolation(NumpyTestCase):
- def check_random(self):
- numpy.random.seed(0)
- for N in [2,3,5,10]:
- A = csr_matrix(rand(N,N))
- for epsilon in [0.0,0.1,0.5,0.8,1.0]:
- S_result = sa_constant_interpolation(A,epsilon)
- S_expected = reference_sa_constant_interpolation(A,epsilon)
- assert_array_equal(S_result.todense(),S_expected.todense())
-
- def check_poisson1D(self):
- for N in [2,3,5,7,10,11,20,21,29,30]:
- A = poisson_problem1D(N)
- for epsilon in [0.0,0.1,0.5,0.8,1.0]:
- S_result = sa_constant_interpolation(A,epsilon)
- S_expected = reference_sa_constant_interpolation(A,epsilon)
- assert_array_equal(S_result.todense(),S_expected.todense())
-
- def check_poisson2D(self):
- for N in [2,3,5,7,10,11,20,21,29,30]:
- A = poisson_problem2D(N)
- for epsilon in [0.0,0.1,0.5,0.8,1.0]:
- S_result = sa_constant_interpolation(A,epsilon)
- S_expected = reference_sa_constant_interpolation(A,epsilon)
- assert_array_equal(S_result.todense(),S_expected.todense())
-
-## def check_sample_data(self):
-## from examples import all_examples,read_matrix
-##
-## for filename in all_examples:
-## A = read_matrix(filename)
-## for epsilon in [0.0,0.08,0.51,1.0]:
-## S_result = sa_constant_interpolation(A,epsilon)
-## S_expected = reference_sa_constant_interpolation(A,epsilon)
-## assert_array_equal((S_result - S_expected).nnz,0)
-
-if __name__ == '__main__':
- NumpyTest().run()
-
Copied: trunk/scipy/sandbox/multigrid/tests/test_sa.py (from rev 3393, trunk/scipy/sandbox/multigrid/tests/test_coarsen.py)
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_coarsen.py 2007-10-03 06:19:08 UTC (rev 3393)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-03 21:37:54 UTC (rev 3395)
@@ -0,0 +1,160 @@
+from numpy.testing import *
+
+from numpy import sqrt,empty,ones,arange,array_split
+from scipy import rand
+from scipy.sparse import spdiags,csr_matrix,lil_matrix
+import numpy
+
+set_package_path()
+import scipy.sandbox.multigrid
+from scipy.sandbox.multigrid.sa import sa_strong_connections, sa_constant_interpolation, \
+ sa_interpolation, sa_fit_candidates
+from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D
+restore_path()
+
+
+def reference_sa_strong_connections(A,epsilon):
+ A_coo = A.tocoo()
+ S = lil_matrix(A.shape)
+
+ for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
+ if i == j: continue #skip diagonal
+
+ if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
+ S[i,j] = v
+
+ return S.tocsr()
+
+class TestSAStrongConnections(NumpyTestCase):
+ def check_simple(self):
+ N = 4
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+ S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+ assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense()) #all connections are strong
+ assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+
+ N = 100
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+ S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+ assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense()) #all connections are strong
+ assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+
+ def check_random(self):
+ numpy.random.seed(0)
+
+ for N in [2,3,5]:
+ A = csr_matrix(rand(N,N))
+ for epsilon in [0.0,0.1,0.5,1.0,10.0]:
+ S_result = sa_strong_connections(A,epsilon)
+ S_expected = reference_sa_strong_connections(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson1D(self):
+ for N in [2,3,5,7,10,11,19]:
+ A = poisson_problem1D(N)
+ for epsilon in [0.0,0.1,0.5,1.0]:
+ S_result = sa_strong_connections(A,epsilon)
+ S_expected = reference_sa_strong_connections(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson2D(self):
+ for N in [2,3,5,7,10,11]:
+ A = poisson_problem2D(N)
+ for epsilon in [0.0,0.1,0.5,1.0]:
+ S_result = sa_strong_connections(A,epsilon)
+ S_expected = reference_sa_strong_connections(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+
+# 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])
+
+ n = A.shape[0]
+
+ R = set(range(n))
+ j = 0
+
+ aggregates = empty(n,dtype=A.indices.dtype)
+ aggregates[:] = -1
+
+ # Pass #1
+ for i,row in enumerate(S):
+ Ni = set(row) | set([i])
+
+ if Ni.issubset(R):
+ R -= Ni
+ for x in Ni:
+ aggregates[x] = j
+ j += 1
+
+ # Pass #2
+ Old_R = R.copy()
+ for i,row in enumerate(S):
+ if i not in R: continue
+
+ for x in row:
+ if x not in Old_R:
+ aggregates[i] = aggregates[x]
+ R.remove(i)
+ break
+
+ # Pass #3
+ for i,row in enumerate(S):
+ if i not in R: continue
+ Ni = set(row) | set([i])
+
+ for x in Ni:
+ if x in R:
+ aggregates[x] = j
+ j += 1
+
+ assert(len(R) == 0)
+
+ Pj = aggregates
+ Pp = arange(n+1)
+ Px = ones(n)
+
+ return csr_matrix((Px,Pj,Pp))
+
+class TestSAConstantInterpolation(NumpyTestCase):
+ def check_random(self):
+ numpy.random.seed(0)
+ for N in [2,3,5,10]:
+ A = csr_matrix(rand(N,N))
+ for epsilon in [0.0,0.1,0.5,1.0]:
+ S_result = sa_constant_interpolation(A,epsilon)
+ S_expected = reference_sa_constant_interpolation(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson1D(self):
+ for N in [2,3,5,7,10,11,20,21,29,30]:
+ A = poisson_problem1D(N)
+ for epsilon in [0.0,0.1,0.5,1.0]:
+ S_result = sa_constant_interpolation(A,epsilon)
+ S_expected = reference_sa_constant_interpolation(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson2D(self):
+ for N in [2,3,5,7,10,11]:
+ A = poisson_problem2D(N)
+ for epsilon in [0.0,0.1,0.5,1.0]:
+ S_result = sa_constant_interpolation(A,epsilon)
+ S_expected = reference_sa_constant_interpolation(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+## def check_sample_data(self):
+## from examples import all_examples,read_matrix
+##
+## for filename in all_examples:
+## A = read_matrix(filename)
+## for epsilon in [0.0,0.08,0.51,1.0]:
+## S_result = sa_constant_interpolation(A,epsilon)
+## S_expected = reference_sa_constant_interpolation(A,epsilon)
+## assert_array_equal((S_result - S_expected).nnz,0)
+
+if __name__ == '__main__':
+ NumpyTest().run()
+
More information about the Scipy-svn
mailing list