[Scipy-svn] r3478 - trunk/scipy/sandbox/multigrid
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Oct 31 21:31:40 EDT 2007
Author: wnbell
Date: 2007-10-31 20:31:37 -0500 (Wed, 31 Oct 2007)
New Revision: 3478
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/dec_test.py
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
Log:
added new method for curl-curl problem
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-11-01 01:31:37 UTC (rev 3478)
@@ -1,18 +1,18 @@
import numpy,scipy,scipy.sparse
from numpy import sqrt, ravel, diff, zeros, zeros_like, inner, concatenate, \
- asarray, hstack, ascontiguousarray, isinf
+ asarray, hstack, ascontiguousarray, isinf, dot
from numpy.random import randn
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,expand_into_blocks
+from utils import approximate_spectral_radius,hstack_csr,vstack_csr,expand_into_blocks,diag_sparse
def augment_candidates(AggOp, old_Q, old_R, new_candidate):
- #TODO update P also
+ #TODO update P and A also
K = old_R.shape[1]
@@ -124,7 +124,7 @@
class adaptive_sa_solver:
- def __init__(self, A, blocks=None, max_levels=10, max_coarse=100,\
+ def __init__(self, A, blocks=None, aggregation=None, max_levels=10, max_coarse=100,\
max_candidates=1, mu=5, epsilon=0.1):
self.A = A
@@ -138,8 +138,9 @@
x,AggOps = self.__initialization_stage(A, blocks = blocks, \
max_levels = max_levels, \
max_coarse = max_coarse, \
- mu = mu, epsilon = epsilon)
-
+ mu = mu, epsilon = epsilon, \
+ aggregation = aggregation )
+
#create SA using x here
As,Ps,Ts,Bs = sa_hierarchy(A,x,AggOps)
@@ -147,8 +148,20 @@
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 = self.__augment_cycle(As,Ps,Ts,Bs,AggOps,x)
+
+ #B = hstack((Bs[0],x))
+ #As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+
+ #improve candidates?
+ if True:
+ print "improving candidates"
+ B = Bs[0]
+ for i in range(max_candidates):
+ B = B[:,1:]
+ As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
+ x = self.__develop_new_candidate(As,Ps,Ts,Bs,AggOps,mu=mu)
+ B = hstack((B,x))
As,Ps,Ts,Bs = sa_hierarchy(A,B,AggOps)
self.Ts = Ts
@@ -156,12 +169,12 @@
self.AggOps = AggOps
self.Bs = Bs
+ def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon,aggregation):
+ if aggregation is not None:
+ max_coarse = 0
+ max_levels = len(aggregation) + 1
- def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
- AggOps = []
- Ps = []
-
- # aSA parameters
+ # aSA parameters
# mu - number of test relaxation iterations
# epsilon - minimum acceptable relaxation convergence factor
@@ -177,9 +190,14 @@
#TODO test convergence rate here
As = [A]
+ AggOps = []
+ Ps = []
- while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
- W_l = sa_constant_interpolation(A_l,epsilon=0,blocks=blocks) #step 4b
+ 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
+ else:
+ W_l = aggregation[len(AggOps)]
P_l,x = sa_fit_candidates(W_l,x) #step 4c
I_l = smoothed_prolongator(P_l,A_l) #step 4d
A_l = I_l.T.tocsr() * A_l * I_l #step 4e
@@ -286,8 +304,11 @@
from multilevel import poisson_problem1D,poisson_problem2D
blocks = None
+ aggregation = None
- A = poisson_problem2D(100)
+ #A = poisson_problem2D(200,1e-2)
+ #aggregation = [ sa_constant_interpolation(A*A*A,epsilon=0.0) ]
+
#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()
@@ -301,17 +322,17 @@
blocks = arange(A.shape[0]/2).repeat(2)
from time import clock; start = clock()
- asa = adaptive_sa_solver(A,max_candidates=5,mu=6,blocks=blocks)
+ asa = adaptive_sa_solver(A,max_candidates=3,mu=5,blocks=blocks,aggregation=aggregation)
print "Adaptive Solver Construction: %s seconds" % (clock() - start); del start
- #scipy.random.seed(0) #make tests repeatable
+ scipy.random.seed(0) #make tests repeatable
x = randn(A.shape[0])
b = A*randn(A.shape[0])
#b = zeros(A.shape[0])
print "solving"
- if True:
+ if False:
x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
else:
residuals = []
@@ -352,11 +373,12 @@
pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
show()
+
for c in asa.Bs[0].T:
+ #plot2d(c)
plot2d_arrows(c)
print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
-
##W = asa.AggOps[0]*asa.AggOps[1]
##pcolor((W * rand(W.shape[1])).reshape((200,200)))
Modified: trunk/scipy/sandbox/multigrid/dec_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/dec_test.py 2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/dec_test.py 2007-11-01 01:31:37 UTC (rev 3478)
@@ -1,10 +1,12 @@
from scipy import *
+from scipy.sparse import *
from pydec import *
-from pydec.multigrid import *
from pydec.multigrid.discrete_laplacian import boundary_hierarchy, discrete_laplacian_solver, hodge_solver
-from scipy.sandbox.multigrid import smoothed_aggregation_solver
+from scipy.sandbox.multigrid import smoothed_aggregation_solver,multigridtools,multilevel_solver
+from scipy.sandbox.multigrid.adaptive import adaptive_sa_solver
+from scipy.sandbox.multigrid.sa import sa_smoothed_prolongator
from scipy.sandbox.multigrid.utils import expand_into_blocks
@@ -13,8 +15,9 @@
#mesh = read_mesh(mesh_path + 'rocket/rocket.xml')
#mesh = read_mesh(mesh_path + 'genus3/genus3_168k.xml')
#mesh = read_mesh(mesh_path + 'genus3/genus3_455k.xml')
-mesh = read_mesh(mesh_path + '/torus/torus.xml')
-for i in range(3):
+#mesh = read_mesh(mesh_path + '/torus/torus.xml')
+mesh = read_mesh(mesh_path + '/sq14tri/sq14tri.xml')
+for i in range(5):
mesh['vertices'],mesh['elements'] = loop_subdivision(mesh['vertices'],mesh['elements'])
cmplx = simplicial_complex(mesh['vertices'],mesh['elements'])
@@ -24,8 +27,34 @@
#bitmap[100:150,100:400] = False
#cmplx = regular_cube_complex(regular_cube_mesh(bitmap))
+def curl_curl_prolongator(D_nodal,vertices):
+ if not isspmatrix_csr(D_nodal):
+ raise TypeError('expected csr_matrix')
+
+ A = D_nodal.T.tocsr() * D_nodal
+ aggs = multigridtools.sa_get_aggregates(A.shape[0],A.indptr,A.indices)
+
+ num_edges = D_nodal.shape[0]
+ num_basis = vertices.shape[1]
+ num_aggs = aggs.max() + 1
+ # replace with CSR + eliminate duplicates
+ #indptr = (2*num_basis) * arange(num_edges+1)
+ ## same same
+ #csr_matrix((data,indices,indptr),dims=(num_edges,num_aggs))
+ row = arange(num_edges).repeat(2*num_basis)
+ col = (num_basis*aggs[D_nodal.indices]).repeat(num_basis)
+ col = col.reshape(-1,num_basis) + arange(num_basis)
+ col = col.reshape(-1)
+ data = tile(0.5 * (D_nodal*vertices),(1,2)).reshape(-1)
+
+ return coo_matrix((data,(row,col)),dims=(num_edges,num_basis*num_aggs)).tocsr()
+
+
+
+
+
def whitney_innerproduct_cache(cmplx,k):
h = hash(cmplx.vertices.tostring()) ^ hash(cmplx.simplices.tostring()) ^ hash(k)
@@ -73,47 +102,83 @@
Mi = whitney_innerproduct_cache(cmplx,i+1)
else:
Mi = regular_cube_innerproduct(cmplx,i+1)
+
+
+ dimension = mesh['vertices'].shape[1]
- ##print "constructing solver"
- ##ss = discrete_laplacian_solver(cochain_complex,len(cochain_complex)-i-1,innerproduct=Mi)
- ##print ss
- ##
- ##print "solving"
- ##x,res = ss.solve(b=zeros(ss.A.shape[0]),x0=rand(ss.A.shape[0]),return_residuals=True)
-
- bh = boundary_hierarchy(cochain_complex)
- while len(bh) < 3:
- bh.coarsen()
- print repr(bh)
-
- N = len(cochain_complex) - 1
-
- B = bh[0][N - i].B
-
- A = B.T.tocsr() * B
- #A = B.T.tocsr() * Mi * B
-
- constant_prolongators = [lvl[N - i].I for lvl in bh[:-1]]
-
- if i == 0:
+ if True:
+
+ d0 = cmplx[0].d
+ d1 = cmplx[1].d
+
+ #A = (d1.T.tocsr() * d1 + d0 * d0.T.tocsr()).astype('d')
+ A = (d1.T.tocsr() * d1).astype('d')
+
+ P = curl_curl_prolongator(d0,mesh['vertices'])
+
+ num_blocks = P.shape[1]/dimension
+ blocks = arange(num_blocks).repeat(dimension)
+
+ P = sa_smoothed_prolongator(A,P,epsilon=0,omega=4.0/3.0)
+
+ PAP = P.T.tocsr() * A * P
+
candidates = None
+ candidates = zeros((num_blocks,dimension,dimension))
+ for n in range(dimension):
+ candidates[:,n,n] = 1.0
+ candidates = candidates.reshape(-1,dimension)
+
+ ml = smoothed_aggregation_solver(PAP,epsilon=0.0,candidates=candidates,blocks=blocks)
+ #A = PAP
+ ml = multilevel_solver([A] + ml.As, [P] + ml.Ps)
else:
- #candidates = [ones(A.shape[0])]
- #TODO test
- candidates = []
- for coord in range(mesh['vertices'].shape[1]):
- candidates.append( bh[0][N-i+1].B * mesh['vertices'][:,coord] )
+ bh = boundary_hierarchy(cochain_complex)
+ while len(bh) < 3:
+ bh.coarsen()
+ print repr(bh)
+
+ N = len(cochain_complex) - 1
+
+ B = bh[0][N - i].B
+
+ A = (B.T.tocsr() * B).astype('d')
+ #A = B.T.tocsr() * Mi * B
+
+ constant_prolongators = [lvl[N - i].I for lvl in bh[:-1]]
+
+ method = 'aSA'
- K = len(candidates)
+ if method == 'RS':
+ As = [A]
+ Ps = []
+ for T in constant_prolongators:
+ Ps.append( sa_smoothed_prolongator(As[-1],T,epsilon=0.0,omega=4.0/3.0) )
+ As.append(Ps[-1].T.tocsr() * As[-1] * Ps[-1])
+ ml = multilevel_solver(As,Ps)
+
+ else:
+ if method == 'BSA':
+ if i == 0:
+ candidates = None
+ else:
+ candidates = cmplx[0].d * mesh['vertices']
+ K = candidates.shape[1]
+
+ constant_prolongators = [constant_prolongators[0]] + \
+ [expand_into_blocks(T,K,1).tocsr() for T in constant_prolongators[1:] ]
- constant_prolongators = [constant_prolongators[0]] + \
- [expand_into_blocks(T,K,1).tocsr() for T in constant_prolongators[1:] ]
+ ml = smoothed_aggregation_solver(A,candidates,aggregation=constant_prolongators)
+ elif method == 'aSA':
+ asa = adaptive_sa_solver(A,aggregation=constant_prolongators,max_candidates=dimension,epsilon=0.0)
+ ml = asa.solver
+ else:
+ raise ValuerError,'unknown method'
+
+ #ml = smoothed_aggregation_solver(A,candidates)
-
- ml = smoothed_aggregation_solver(A,candidates,aggregation=constant_prolongators)
- #ml = smoothed_aggregation_solver(A,candidates)
-
+ #x = d0 * mesh['vertices'][:,0]
x = rand(A.shape[0])
b = zeros_like(x)
#b = A*rand(A.shape[0])
@@ -125,9 +190,11 @@
def add_resid(x):
residuals.append(linalg.norm(b - A*x))
A.psolve = ml.psolve
- x_sol = linalg.cg(A,b,x0=x,maxiter=30,tol=1e-12,callback=add_resid)[0]
+ from pydec import cg
+ x_sol = cg(A,b,x0=x,maxiter=40,tol=1e-8,callback=add_resid)[0]
+
residuals = array(residuals)/residuals[0]
avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
print "average convergence ratio",avg_convergence_ratio
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-11-01 01:31:37 UTC (rev 3478)
@@ -24,14 +24,14 @@
O = -ones(N)
return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate explicit zeros
-def poisson_problem2D(N):
+def poisson_problem2D(N,epsilon=1.0):
"""
Return a sparse CSR matrix for the 2d poisson problem
with standard 5-point finite difference stencil on a
square N-by-N grid.
"""
- D = 4*ones(N*N)
- T = -ones(N*N)
+ D = (2 + 2*epsilon)*ones(N*N)
+ T = -epsilon * 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 explicit zeros
@@ -208,7 +208,9 @@
#use direct solver on coarsest level
#TODO reuse factors for efficiency?
coarse_x[:] = spsolve(self.As[-1],coarse_b).reshape(coarse_x.shape)
- #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0]
+ #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0].reshape(coarse_x.shape)
+ #A_inv = asarray(scipy.linalg.pinv2(self.As[-1].todense()))
+ #coarse_x[:] = scipy.dot(A_inv,coarse_b)
#print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
else:
self.__solve(lvl+1,coarse_x,coarse_b)
@@ -230,18 +232,17 @@
from scipy import *
candidates = None
blocks = None
- #A = poisson_problem2D(100)
+ A = poisson_problem2D(40,1e-2)
#A = io.mmread("rocker_arm_surface.mtx").tocsr()
#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]) ]
- blocks = arange(A.shape[0]/2).repeat(2)
+ #A = io.mmread('tests/sample_data/elas30_A.mtx').tocsr()
+ #candidates = io.mmread('tests/sample_data/elas30_nullspace.mtx')
+ #blocks = arange(A.shape[0]/2).repeat(2)
- ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=100,max_levels=2)
+ ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0.08,max_coarse=100,max_levels=10)
#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-31 21:46:29 UTC (rev 3477)
+++ trunk/scipy/sandbox/multigrid/sa.py 2007-11-01 01:31:37 UTC (rev 3478)
@@ -182,7 +182,7 @@
P = T - (D_inv_A*T)
#S = (spidentity(A.shape[0]).tocsr() - D_inv_A) #TODO drop this?
- #P = S * ( S * T)
+ #P = S *(S * ( S * T))
return P
@@ -201,8 +201,6 @@
T,coarse_candidates = sa_fit_candidates(AggOp,candidates)
#T = AggOp #TODO test
- A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
-
P = sa_smoothed_prolongator(A,T,epsilon,omega,blocks)
if blocks is not None:
More information about the Scipy-svn
mailing list