[Scipy-svn] r3852 - in trunk/scipy/sandbox/multigrid: . examples tests
scipy-svn@scip...
scipy-svn@scip...
Sat Jan 19 13:43:07 CST 2008
Author: wnbell
Date: 2008-01-19 13:43:04 -0600 (Sat, 19 Jan 2008)
New Revision: 3852
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/examples/adaptive.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
updated aSA
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2008-01-19 19:43:04 UTC (rev 3852)
@@ -6,7 +6,8 @@
from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
asarray, hstack, ascontiguousarray, isinf, dot
from numpy.random import randn, rand
-from scipy.sparse import csr_matrix,coo_matrix
+from scipy.linalg import norm
+from scipy.sparse import csr_matrix, coo_matrix, bsr_matrix
from relaxation import gauss_seidel
from multilevel import multilevel_solver
@@ -40,7 +41,7 @@
for AggOp in AggOps:
P,B = sa_fit_candidates(AggOp,B)
I = sa_smoothed_prolongator(A,P)
- A = I.T.tocsr() * A * I
+ A = I.T.asformat(I.format) * A * I
As.append(A)
Ts.append(P)
Ps.append(I)
@@ -99,18 +100,22 @@
"""
if A.shape[0] <= max_coarse:
- raise ValueError,'small matrices not handled yet'
+ return multilevel_solver( [A], [] )
#first candidate
x,AggOps = asa_initial_setup_stage(A, max_levels = max_levels, \
max_coarse = max_coarse, mu = mu, epsilon = epsilon, \
aggregation = aggregation )
+ #TODO make sa_fit_candidates work for small Bs
+ x /= norm(x)
+
#create SA using x here
As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)
for i in range(max_candidates - 1):
x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu=mu)
+ x /= norm(x)
B = hstack((Bs[0],x))
As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
@@ -122,7 +127,7 @@
for i in range(max_candidates):
B = B[:,1:]
As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
- x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu=mu)
+ x = asa_general_setup_stage(As,Ps,Ts,Bs,AggOps,mu)
B = hstack((B,x))
As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
@@ -159,7 +164,7 @@
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.tocsr() * A_l * I_l #step 4e
+ A_l = I_l.T.asformat(I_l.format) * A_l * I_l #step 4e
AggOps.append(W_l)
Ps.append(I_l)
@@ -189,7 +194,7 @@
def asa_general_setup_stage(As, Ps, Ts, Bs, AggOps, mu):
A = As[0]
- x = randn(A.shape[0],1)
+ x = rand(A.shape[0],1)
b = zeros_like(x)
x = multilevel_solver(As,Ps).solve(b, x0=x, tol=1e-10, maxiter=mu)
@@ -204,15 +209,20 @@
K = P.blocksize[0]
bnnz = P.indptr[-1]
data = zeros( (bnnz, K+1, K), dtype=P.dtype )
- data[:,:-1,:-1] = P.data
+ data[:,:-1,:] = P.data
return bsr_matrix( (data, P.indices, P.indptr), shape=( (K+1)*(M/K), N) )
for i in range(len(As) - 2):
- B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
- T,R = sa_fit_candidates(AggOp,B)
+ B_old = Bs[i]
+ B = zeros( (x.shape[0], B_old.shape[1] + 1), dtype=x.dtype)
+ B[:B_old.shape[0],:B_old.shape[1]] = B_old
+ B[:,-1] = x.reshape(-1)
+
+ T,R = sa_fit_candidates(AggOps[i],B)
+
P = sa_smoothed_prolongator(A,T)
- A = P.T.tocsr() * A * P
+ A = P.T.asformat(P.format) * A * P
temp_Ps.append(P)
temp_As.append(A)
@@ -242,7 +252,7 @@
# T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
#
# P = sa_smoothed_prolongator(A,T)
-# A = P.T.tocsr() * A * P
+# A = P.T.asformat(P.format) * A * P
#
# new_As.append(A)
# new_Ps.append(P)
Modified: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py 2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py 2008-01-19 19:43:04 UTC (rev 3852)
@@ -3,21 +3,25 @@
from scipy.sandbox.multigrid.utils import diag_sparse
from scipy.sandbox.multigrid.gallery import poisson, linear_elasticity
from scipy.sandbox.multigrid.adaptive import adaptive_sa_solver
+from scipy.sandbox.multigrid import smoothed_aggregation_solver
-#A,B = linear_elasticity( (100,100) )
+A,B = linear_elasticity( (20,20) )
-A = poisson( (200,200), format='csr' )
+#A = poisson( (200,200), format='csr' )
#A = poisson( (200,200), spacing=(1,1e-2) ) #anisotropic
-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
from time import clock; start = clock()
-asa = adaptive_sa_solver(A, max_levels=2, max_candidates=1, mu=10)
+#asa = smoothed_aggregation_solver(A,B,max_levels=2)
+#Bs = [B]
+asa,Bs = adaptive_sa_solver(A, max_levels=2, max_coarse=10, max_candidates=3, mu=15)
+
print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
random.seed(0) #make tests repeatable
@@ -27,14 +31,14 @@
print "solving"
-if True:
- x_sol,residuals = asa.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
+if False:
+ x_sol,residuals = asa.solve(b,x0=x,maxiter=20,tol=1e-10,return_residuals=True)
else:
residuals = []
def add_resid(x):
residuals.append(linalg.norm(b - A*x))
A.psolve = asa.psolve
- x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
+ x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-10,callback=add_resid)[0]
residuals = array(residuals)/residuals[0]
@@ -69,10 +73,10 @@
show()
-#for c in asa.Bs[0].T:
-# plot2d(c)
-# #plot2d_arrows(c)
-# print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+for c in Bs[0].T:
+ #plot2d(c)
+ plot2d_arrows(c)
+ print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
#plot2d(x_sol)
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/sa.py 2008-01-19 19:43:04 UTC (rev 3852)
@@ -103,6 +103,13 @@
if candidates.dtype != 'float32':
candidates = asarray(candidates,dtype='float64')
+ if len(candidates.shape) != 2:
+ raise ValueError,'expected rank 2 array for argument B'
+
+ if candidates.shape[0] % AggOp.shape[0] != 0:
+ raise ValueError,'dimensions of AggOp %s and B %s are incompatible' % (AggOp.shape, B.shape)
+
+
K = candidates.shape[1] # number of near-nullspace candidates
blocksize = candidates.shape[0] / AggOp.shape[0]
@@ -167,8 +174,11 @@
A_filtered = sa_filtered_matrix(A,epsilon) #use filtered matrix for anisotropic problems
# TODO use scale_rows()
- D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
- D_inv_A = D_inv * A_filtered
+ D = diag_sparse(A_filtered)
+ D_inv = 1.0 / D
+ D_inv[D == 0] = 0
+
+ D_inv_A = diag_sparse(D_inv) * A_filtered
D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
# smooth tentative prolongator T
Modified: trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_adaptive.py 2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/tests/test_adaptive.py 2008-01-19 19:43:04 UTC (rev 3852)
@@ -1,18 +1,56 @@
from scipy.testing import *
+from numpy import arange, ones, zeros, array, eye, vstack, diff
+from scipy import rand
from scipy.sparse import csr_matrix
-from numpy import arange,ones,zeros,array,eye,vstack,diff
from scipy.sandbox.multigrid.sa import sa_fit_candidates
+from scipy.sandbox.multigrid import smoothed_aggregation_solver
#from scipy.sandbox.multigrid.adaptive import augment_candidates
+from scipy.sandbox.multigrid.gallery import *
+from scipy.sandbox.multigrid.adaptive import *
class TestAdaptiveSA(TestCase):
def setUp(self):
- pass
+ from numpy.random import seed
+ seed(0)
+ def test_poisson(self):
+ A = poisson( (100,100), format='csr' )
+
+ asa = adaptive_sa_solver(A, max_candidates = 1)
+ sa = smoothed_aggregation_solver(A, B = ones((A.shape[0],1)) )
+
+ b = rand(A.shape[0])
+
+ sol0,residuals0 = asa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+ sol1,residuals1 = sa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+
+ conv_asa = (residuals0[-1]/residuals0[0])**(1.0/len(residuals0))
+ conv_sa = (residuals1[-1]/residuals1[0])**(1.0/len(residuals1))
+
+ assert( conv_asa < 1.1 * conv_sa ) #aSA shouldn't be any worse than SA
+
+# def test_elasticity(self):
+# A,B = linear_elasticity( (100,100), format='bsr' )
+#
+# asa = adaptive_sa_solver(A, max_candidates = 3)
+# sa = smoothed_aggregation_solver(A, B=B )
+#
+# b = rand(A.shape[0])
+#
+# sol0,residuals0 = asa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+# sol1,residuals1 = sa.solve(b, maxiter=20, tol=1e-10, return_residuals=True)
+#
+# conv_asa = (residuals0[-1]/residuals0[0])**(1.0/len(residuals0))
+# conv_sa = (residuals1[-1]/residuals1[0])**(1.0/len(residuals1))
+#
+# print "ASA convergence",conv_asa
+# assert( conv_asa < 1.1 * conv_sa ) #aSA shouldn't be any worse than SA
+
#class TestAugmentCandidates(TestCase):
# def setUp(self):
# self.cases = []
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2008-01-19 15:39:51 UTC (rev 3851)
+++ trunk/scipy/sandbox/multigrid/utils.py 2008-01-19 19:43:04 UTC (rev 3852)
@@ -46,6 +46,8 @@
if A.shape[0] != A.shape[1]:
raise ValueError,'expected square matrix'
+ maxiter = min(A.shape[0],maxiter)
+
#TODO make method adaptive
numpy.random.seed(0) #make results deterministic
More information about the Scipy-svn
mailing list