[Scipy-svn] r3427 - in trunk/scipy/sandbox/multigrid: . multigridtools tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Oct 9 17:30:34 EDT 2007
Author: wnbell
Date: 2007-10-09 16:30:30 -0500 (Tue, 09 Oct 2007)
New Revision: 3427
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/tests/test_sa.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
fixed bug in approximate_spectral_radius where A was assumed to be symmetric
adaptive SA works for symmetric diagonal scaled case now
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-09 21:30:30 UTC (rev 3427)
@@ -4,7 +4,7 @@
from relaxation import gauss_seidel
from multilevel import multilevel_solver
-from sa import sa_constant_interpolation
+from sa import sa_constant_interpolation,sa_fit_candidates
#from utils import infinity_norm
from utils import approximate_spectral_radius
@@ -144,7 +144,7 @@
Ps = []
for W in Ws:
- P,x = fit_candidates(W,x)
+ P,x = sa_fit_candidates(W,x)
I = smoothed_prolongator(P,A)
A = I.T.tocsr() * A * I
As.append(A)
@@ -301,7 +301,8 @@
while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
#W_l = sa_constant_interpolation(A_l,epsilon=0.08*0.5**(len(AggOps)-1)) #step 4b #TEST
W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b
- P_l,x = fit_candidate(W_l,x) #step 4c
+ P_l,x = sa_fit_candidates(W_l,[x]) #step 4c
+ x = x[0] #TODO make sa_fit_candidates accept a single x
I_l = smoothed_prolongator(P_l,A_l) #step 4d
A_l = I_l.T.tocsr() * A_l * I_l #step 4e
@@ -337,15 +338,15 @@
from scipy import *
from utils import diag_sparse
from multilevel import poisson_problem1D,poisson_problem2D
-A = poisson_problem2D(50)
+A = poisson_problem2D(200)
#A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
#A = io.mmread("laplacian_40_3dcube.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 = A*A
-#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-#A = D * A * D
+D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+A = D * A * D
#A = io.mmread("nos2.mtx").tocsr()
asa = adaptive_sa_solver(A,max_candidates=1)
#x = arange(A.shape[0]).astype('d') + 1
Modified: trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h 2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h 2007-10-09 21:30:30 UTC (rev 3427)
@@ -42,18 +42,26 @@
T eps_Aii = epsilon*epsilon*diags[i];
+ T weak_sum = 0.0;
+
for(int jj = row_start; jj < row_end; jj++){
const int j = Aj[jj];
const T Aij = Ax[jj];
- if(i == j){continue;}
+ if(i == j){continue;} //skip diagonal until end of row
// |A(i,j)| < epsilon * sqrt(|A(i,i)|*|A(j,j)|)
if(Aij*Aij >= std::abs(eps_Aii * diags[j])){
Sj->push_back(j);
Sx->push_back(Aij);
+ } else {
+ weak_sum += Aij;
}
}
+ //Add modified diagonal entry
+ Sj->push_back(i);
+ Sx->push_back(diags[i] + weak_sum); //filtered matrix
+
Sp->push_back(Sj->size());
}
}
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-09 21:30:30 UTC (rev 3427)
@@ -4,8 +4,9 @@
import scipy
import numpy
-from numpy import zeros,zeros_like,array
+from numpy import ones,zeros,zeros_like,array
from numpy.linalg import norm
+from scipy.linsolve import spsolve
from sa import sa_interpolation
from rs import rs_interpolation
@@ -19,8 +20,8 @@
with standard 3-point finite difference stencil on a
grid with N points.
"""
- D = 2*numpy.ones(N)
- O = -numpy.ones(N)
+ D = 2*ones(N)
+ O = -ones(N)
return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate zeros
def poisson_problem2D(N):
@@ -29,9 +30,9 @@
with standard 5-point finite difference stencil on a
square N-by-N grid.
"""
- D = 4*numpy.ones(N*N)
- T = -numpy.ones(N*N)
- O = -numpy.ones(N*N)
+ D = 4*ones(N*N)
+ T = -ones(N*N)
+ O = -ones(N*N)
T[N-1::N] = 0
return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate zeros
@@ -130,12 +131,12 @@
#TODO change use of tol (relative tolerance) to agree with other iterative solvers
A = self.As[0]
- residuals = [scipy.linalg.norm(b-A*x)]
+ residuals = [ norm(b-A*x) ]
while len(residuals) <= maxiter and residuals[-1]/residuals[0] > tol:
self.__solve(0,x,b)
- residuals.append(scipy.linalg.norm(b-A*x))
+ residuals.append( norm(b-A*x) )
if callback is not None:
callback(x)
@@ -150,7 +151,7 @@
A = self.As[lvl]
if len(self.As) == 1:
- x[:] = scipy.linsolve.spsolve(A,b)
+ x[:] = spsolve(A,b)
return
self.presmoother(A,x,b)
@@ -162,7 +163,7 @@
if lvl == len(self.As) - 2:
#use direct solver on coarsest level
- coarse_x[:] = scipy.linsolve.spsolve(self.As[-1],coarse_b)
+ coarse_x[:] = spsolve(self.As[-1],coarse_b)
#coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0]
#print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
else:
Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-09 21:30:30 UTC (rev 3427)
@@ -1,43 +1,60 @@
from numpy.testing import *
-from numpy import sqrt,empty,ones,arange,array_split,eye,array,zeros,diag
+from numpy import sqrt,empty,ones,arange,array_split,eye,array, \
+ zeros,diag,zeros_like
+from numpy.linalg import norm
from scipy import rand
-from scipy.sparse import spdiags,csr_matrix,lil_matrix
+from scipy.sparse import spdiags,csr_matrix,lil_matrix, \
+ isspmatrix_csr,isspmatrix_csc,isspmatrix_coo, \
+ isspmatrix_lil
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
+from scipy.sandbox.multigrid.multilevel import poisson_problem1D,poisson_problem2D, \
+ smoothed_aggregation_solver
+from scipy.sandbox.multigrid.utils import diag_sparse
restore_path()
+#def sparsity(A):
+# A = A.copy()
+#
+# if isspmatrix_csr(A) or isspmatrix_csc(A) or isspmatrix_coo(A):
+# A.data[:] = 1
+# elif isspmatrix_lil:
+# for row in A.data:
+# row[:] = [1]*len(row)
+# else:
+# raise ValueError,'expected csr,csc,coo, or lil'
+#
+# return A
+
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 i == j or abs(v) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
+ S[i,j] += v
+ else:
+ S[i,i] += v
- if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
- S[i,j] = v
+ return S
- 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_simple(self):
+# N = 4
+# A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+# assert_array_equal(sa_strong_connections(A,0.50).todense(),A.todense()) #all connections are strong
+# assert_array_equal(sa_strong_connections(A,0.51).todense(),0*A.todense()) #no connections are strong
+#
+# N = 100
+# A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+# assert_array_equal(sa_strong_connections(A,0.50).todense(),A.todense()) #all connections are strong
+# assert_array_equal(sa_strong_connections(A,0.51).todense(),0*A.todense()) #no connections are strong
def check_random(self):
numpy.random.seed(0)
@@ -47,7 +64,8 @@
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())
+ assert_almost_equal(S_result.todense(),S_expected.todense())
+ #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
def check_poisson1D(self):
for N in [2,3,5,7,10,11,19]:
@@ -56,6 +74,7 @@
S_result = sa_strong_connections(A,epsilon)
S_expected = reference_sa_strong_connections(A,epsilon)
assert_array_equal(S_result.todense(),S_expected.todense())
+ #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
def check_poisson2D(self):
for N in [2,3,5,7,10,11]:
@@ -64,6 +83,7 @@
S_result = sa_strong_connections(A,epsilon)
S_expected = reference_sa_strong_connections(A,epsilon)
assert_array_equal(S_result.todense(),S_expected.todense())
+ #assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
@@ -209,7 +229,89 @@
+class TestSASolver(NumpyTestCase):
+ def setUp(self):
+ self.cases = []
+ #self.cases.append((poisson_problem1D(10),None))
+
+ self.cases.append((poisson_problem1D(500),None))
+ self.cases.append((poisson_problem2D(50),None))
+
+ def check_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)
+
+ numpy.random.seed(0) #make tests repeatable
+
+ x = rand(A.shape[0])
+ b = A*rand(A.shape[0]) #zeros_like(x)
+
+ x_sol,residuals = ml.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
+
+ avg_convergence_ratio = (residuals[-1]/residuals[0])**(1.0/len(residuals))
+
+ assert(avg_convergence_ratio < 0.5)
+
+ def check_DAD(self):
+
+ for A,candidates in self.cases:
+
+ x = rand(A.shape[0])
+ b = A*rand(A.shape[0]) #zeros_like(x)
+
+ D = diag_sparse(rand(A.shape[0]))
+ D_inv = diag_sparse(1.0/D.data)
+ DAD = D*A*D
+
+ if candidates is None:
+ candidates = [ ones(A.shape[0]) ]
+
+ DAD_candidates = [ (D_inv * c) for c in candidates ]
+
+ #ml = smoothed_aggregation_solver(A,candidates,max_coarse=1,max_levels=2)
+
+ ml = smoothed_aggregation_solver(DAD,DAD_candidates,max_coarse=100,max_levels=2)
+
+ #print (D_inv*ml.Ps[0]).todense()
+
+ x_sol,residuals = ml.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
+
+ avg_convergence_ratio = (residuals[-1]/residuals[0])**(1.0/len(residuals))
+ print avg_convergence_ratio
+
+ assert(avg_convergence_ratio < 0.5)
+
+## def check_DAD(self):
+## """check that method is invariant to symmetric diagonal scaling (A -> DAD)"""
+##
+## for A,A_candidates in self.cases:
+## numpy.random.seed(0) #make tests repeatable
+##
+## x = rand(A.shape[0])
+## b = zeros_like(x)
+##
+## D = diag_sparse(rand(A.shape[0]))
+## D_inv = diag_sparse(1.0/D.data)
+## DAD = D*A*D
+##
+##
+## if A_candidates is None:
+## A_candidates = [ ones(A.shape[0]) ]
+##
+## DAD_candidates = [ (D_inv * c) for c in A_candidates ]
+##
+## ml_A = smoothed_aggregation_solver(A, A_candidates, max_coarse=10, max_levels=10, epsilon=0.0)
+## x_sol_A = ml_A.solve(b, x0=x, maxiter=1, tol=1e-12)
+##
+## ml_DAD = smoothed_aggregation_solver(DAD, DAD_candidates, max_coarse=10, max_levels=10, epsilon=0.0)
+## x_sol_DAD = ml_DAD.solve(b, x0=D*x, maxiter=1, tol=1e-12)
+##
+## assert_almost_equal(x_sol_A, D_inv * x_sol_DAD)
+
+
if __name__ == '__main__':
NumpyTest().run()
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-10-08 22:54:53 UTC (rev 3426)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-10-09 21:30:30 UTC (rev 3427)
@@ -10,8 +10,8 @@
"""
Approximate the spectral radius of a symmetric matrix using ARPACK
"""
- from scipy.sandbox.arpack import eigen_symmetric
- return eigen_symmetric(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0]
+ from scipy.sandbox.arpack import eigen
+ return eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0]
More information about the Scipy-svn
mailing list