[Scipy-svn] r3838 - in trunk/scipy/sandbox/multigrid: . examples tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Jan 15 11:15:17 EST 2008
Author: wnbell
Date: 2008-01-15 10:15:12 -0600 (Tue, 15 Jan 2008)
New Revision: 3838
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
Log:
updated adaptive SA to use BSR
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2008-01-15 16:15:12 UTC (rev 3838)
@@ -1,3 +1,6 @@
+__all__ = ['adaptive_sa_solver']
+
+
import numpy,scipy,scipy.sparse
from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
@@ -2,9 +5,10 @@
asarray, hstack, ascontiguousarray, isinf, dot
-from numpy.random import randn
+from numpy.random import randn, rand
from scipy.sparse import csr_matrix,coo_matrix
from relaxation import gauss_seidel
from multilevel import multilevel_solver
-from sa import sa_constant_interpolation,sa_fit_candidates
from utils import approximate_spectral_radius,hstack_csr,vstack_csr,diag_sparse
+from sa import sa_constant_interpolation, sa_smoothed_prolongator, \
+ sa_fit_candidates
@@ -19,13 +23,14 @@
Construct multilevel hierarchy using Smoothed Aggregation
Inputs:
A - matrix
- Ps - list of constant prolongators
- B - "candidate" basis function to be approximated
+ B - fine-level near nullspace candidates be approximated
+
Ouputs:
- (As,Ps,Ts) - tuple of lists
- - As - [A, Ts[0].T*A*Ts[0], Ts[1].T*A*Ts[1], ... ]
+ (As,Ps,Ts,Bs) - tuple of lists
+ - As -
- Ps - smoothed prolongators
- Ts - tentative prolongators
+ - Bs - near nullspace candidates
"""
As = [A]
Ps = []
@@ -34,7 +39,7 @@
for AggOp in AggOps:
P,B = sa_fit_candidates(AggOp,B)
- I = sa_smoothed_prolongator(P,A)
+ I = sa_smoothed_prolongator(A,P)
A = I.T.tocsr() * A * I
As.append(A)
Ts.append(P)
@@ -44,19 +49,17 @@
class adaptive_sa_solver:
- def __init__(self, A, blocks=None, aggregation=None, max_levels=10, max_coarse=100,\
+ def __init__(self, A, aggregation=None, max_levels=10, max_coarse=100,\
max_candidates=1, mu=5, epsilon=0.1):
-
- self.A = A
-
+
+ self.A = A
self.Rs = []
if self.A.shape[0] <= max_coarse:
raise ValueError,'small matrices not handled yet'
#first candidate
- x,AggOps = self.__initialization_stage(A, blocks = blocks, \
- max_levels = max_levels, \
+ x,AggOps = self.__initialization_stage(A, max_levels = max_levels, \
max_coarse = max_coarse, \
mu = mu, epsilon = epsilon, \
aggregation = aggregation )
@@ -67,9 +70,6 @@
for i in range(max_candidates - 1):
x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
- #TODO which is faster?
- #As,Ps,Ts,Bs = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
-
B = hstack((Bs[0],x))
As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
@@ -89,7 +89,7 @@
self.AggOps = AggOps
self.Bs = Bs
- def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon,aggregation):
+ def __initialization_stage(self,A,max_levels,max_coarse,mu,epsilon,aggregation):
if aggregation is not None:
max_coarse = 0
max_levels = len(aggregation) + 1
@@ -100,31 +100,28 @@
#step 1
A_l = A
- x = randn(A_l.shape[0],1)
+ x = rand(A_l.shape[0],1) # TODO see why randn() fails here
skip_f_to_i = False
#step 2
- gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')
+ gauss_seidel(A_l, x, zeros_like(x), iterations=mu, sweep='symmetric')
#step 3
#TODO test convergence rate here
- As = [A]
+ As = [A]
+ Ps = []
AggOps = []
- Ps = []
while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
if aggregation is None:
- W_l = sa_constant_interpolation(A_l,epsilon=0,blocks=blocks) #step 4b
+ W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b
else:
- W_l = aggregation[len(AggOps)]
- P_l,x = sa_fit_candidates(W_l,x) #step 4c
- I_l = sa_smoothed_prolongator(P_l,A_l) #step 4d
- A_l = I_l.T.tocsr() * A_l * I_l #step 4e
+ 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
- blocks = None #not needed on subsequent levels
-
- print 'A_l.shape',A_l.shape
AggOps.append(W_l)
Ps.append(I_l)
As.append(A_l)
@@ -132,7 +129,6 @@
if A_l.shape <= max_coarse: break
if not skip_f_to_i:
- print "."
x_hat = x.copy() #step 4g
gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric') #step 4h
x_A_x = dot(x.T,A_l*x)
@@ -176,7 +172,7 @@
B = zeros( (x.shape[0], Bs[i+1].shape[1] + 1), dtype=x.dtype)
T,R = sa_fit_candidates(AggOp,B)
- P = sa_smoothed_prolongator(T,A)
+ P = sa_smoothed_prolongator(A,T)
A = P.T.tocsr() * A * P
temp_Ps.append(P)
@@ -206,7 +202,7 @@
# for i in range(len(As) - 1):
# T,R = augment_candidates(AggOps[i], Ts[i], Bs[i+1], x)
#
-# P = sa_smoothed_prolongator(T,A)
+# P = sa_smoothed_prolongator(A,T)
# A = P.T.tocsr() * A * P
#
# new_As.append(A)
@@ -219,9 +215,6 @@
# return new_As,new_Ps,new_Ts,new_Bs
-
-
-
#def augment_candidates(AggOp, old_Q, old_R, new_candidate):
# #TODO update P and A also
#
Modified: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py 2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py 2008-01-15 16:15:12 UTC (rev 3838)
@@ -2,36 +2,32 @@
from scipy import *
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
+#A,B = linear_elasticity( (100,100) )
-#A = poisson( (200,200), spacing=(1,1e-2)
-#aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ]
+A = poisson( (200,200), format='csr' )
-#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 = 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
-#A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
-A,B = linear_elasticity( (100,100) )
from time import clock; start = clock()
-asa = adaptive_sa_solver(A,max_candidates=3,mu=5,blocks=blocks,aggregation=aggregation)
+
+asa = adaptive_sa_solver(A, max_levels=2, max_candidates=1, mu=10)
+
print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
-scipy.random.seed(0) #make tests repeatable
+random.seed(0) #make tests repeatable
x = randn(A.shape[0])
-b = A*randn(A.shape[0])
-#b = zeros(A.shape[0])
+#b = A*randn(A.shape[0])
+b = zeros(A.shape[0])
print "solving"
-if False:
+if True:
x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
else:
residuals = []
@@ -74,10 +70,11 @@
for c in asa.Bs[0].T:
- #plot2d(c)
- plot2d_arrows(c)
+ plot2d(c)
+ #plot2d_arrows(c)
print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+#plot2d(x_sol)
##W = asa.AggOps[0]*asa.AggOps[1]
##pcolor((W * rand(W.shape[1])).reshape((200,200)))
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/sa.py 2008-01-15 16:15:12 UTC (rev 3838)
@@ -97,6 +97,9 @@
sa_constant_interpolation(csr_matrix(A),epsilon)
def sa_fit_candidates(AggOp,candidates,tol=1e-10):
+ if not isspmatrix_csr(AggOp):
+ raise TypeError,'expected csr_matrix for argument AggOp'
+
if candidates.dtype != 'float32':
candidates = asarray(candidates,dtype='float64')
@@ -147,7 +150,7 @@
return Q,R
-def sa_smoothed_prolongator(A,T,epsilon,omega):
+def sa_smoothed_prolongator(A,T,epsilon=0.0,omega=4.0/3.0):
"""For a given matrix A and tentative prolongator T return the
smoothed prolongator P
@@ -161,7 +164,6 @@
A_filtered - sa_filtered_matrix(A,epsilon)
"""
-
A_filtered = sa_filtered_matrix(A,epsilon) #use filtered matrix for anisotropic problems
# TODO use scale_rows()
Modified: trunk/scipy/sandbox/multigrid/tests/test_adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_adaptive.py 2008-01-15 13:54:53 UTC (rev 3837)
+++ trunk/scipy/sandbox/multigrid/tests/test_adaptive.py 2008-01-15 16:15:12 UTC (rev 3838)
@@ -5,7 +5,7 @@
from scipy.sandbox.multigrid.sa import sa_fit_candidates
-from scipy.sandbox.multigrid.adaptive import augment_candidates
+#from scipy.sandbox.multigrid.adaptive import augment_candidates
More information about the Scipy-svn
mailing list