[Scipy-svn] r3445 - trunk/scipy/sandbox/multigrid
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Oct 18 04:06:14 EDT 2007
Author: wnbell
Date: 2007-10-18 03:06:05 -0500 (Thu, 18 Oct 2007)
New Revision: 3445
Added:
trunk/scipy/sandbox/multigrid/dec_test.py
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
fixed bug in expansion of AggOp
added dec example, incomplete
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-18 08:06:05 UTC (rev 3445)
@@ -294,7 +294,7 @@
print "solving"
if True:
- x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
+ x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=25,tol=1e-8,return_residuals=True)
else:
residuals = []
def add_resid(x):
Added: trunk/scipy/sandbox/multigrid/dec_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/dec_test.py 2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/dec_test.py 2007-10-18 08:06:05 UTC (rev 3445)
@@ -0,0 +1,174 @@
+
+from scipy 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.utils import expand_into_blocks
+
+
+## Load mesh from file
+mesh_path = '../../../../../hirani_group/wnbell/meshes/'
+#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['vertices'],mesh['elements'] = loop_subdivision(mesh['vertices'],mesh['elements'])
+cmplx = simplicial_complex(mesh['vertices'],mesh['elements'])
+
+## Construct mesh manually
+#bitmap = ones((60,60),dtype='bool')
+#bitmap[1::10,1::10] = False
+#bitmap[100:150,100:400] = False
+#cmplx = regular_cube_complex(regular_cube_mesh(bitmap))
+
+
+
+def whitney_innerproduct_cache(cmplx,k):
+ h = hash(cmplx.vertices.tostring()) ^ hash(cmplx.simplices.tostring()) ^ hash(k)
+
+ filename = "/home/nathan/.pydec/cache/whitney_" + str(h) + ".mtx"
+
+ try:
+ import pydec
+ M = pydec.io.read_array(filename)
+ except:
+ import pydec
+ M = whitney_innerproduct(cmplx,k)
+ pydec.io.write_array(filename,M)
+
+ return M
+
+
+
+def cube_innerproduct_cache(cmplx,k):
+ h = hash(cmplx.mesh.bitmap.tostring()) ^ hash(cmplx.mesh.bitmap.shape) ^ hash(k)
+
+ filename = "/home/nathan/.pydec/cache/cube_" + str(h) + ".mtx"
+
+ try:
+ import pydec
+ M = pydec.io.read_array(filename)
+ except:
+ import pydec
+ M = regular_cube_innerproduct(cmplx,k)
+ pydec.io.write_array(filename,M)
+
+ return M
+
+
+
+#solve d_k d_k problem for all reasonable k
+#from pylab import semilogy,show,xlabel,ylabel,legend,ylim,xlim
+#from matplotlib.font_manager import fontManager, FontProperties
+
+cochain_complex = cmplx.cochain_complex()
+
+for i in [1]: #range(len(cochain_complex)-1):
+ print "computing mass matrix"
+
+ if isinstance(cmplx,simplicial_complex):
+ Mi = whitney_innerproduct_cache(cmplx,i+1)
+ else:
+ Mi = regular_cube_innerproduct(cmplx,i+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:
+ candidates = None
+ 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] )
+
+ K = len(candidates)
+
+ 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)
+ #ml = smoothed_aggregation_solver(A,candidates)
+
+ x = rand(A.shape[0])
+ b = zeros_like(x)
+ #b = A*rand(A.shape[0])
+
+ if True:
+ x_sol,residuals = ml.solve(b,x0=x,maxiter=50,tol=1e-12,return_residuals=True)
+ else:
+ residuals = []
+ 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]
+
+
+ residuals = array(residuals)/residuals[0]
+ avg_convergence_ratio = residuals[-1]**(1.0/len(residuals))
+ print "average convergence ratio",avg_convergence_ratio
+ print "last convergence ratio",residuals[-1]/residuals[-2]
+
+ print residuals
+
+
+
+
+##candidates = None
+##blocks = None
+##
+##
+##
+##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)
+##
+##ml = smoothed_aggregation_solver(A,candidates,blocks=blocks,epsilon=0,max_coarse=10,max_levels=10)
+###ml = ruge_stuben_solver(A)
+##
+##x = rand(A.shape[0])
+###b = zeros_like(x)
+##b = A*rand(A.shape[0])
+##
+##if True:
+## x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-12,return_residuals=True)
+##else:
+## residuals = []
+## def add_resid(x):
+## residuals.append(linalg.norm(b - A*x))
+## A.psolve = ml.psolve
+## x_sol = linalg.cg(A,b,x0=x,maxiter=25,tol=1e-12,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
+##print "last convergence ratio",residuals[-1]/residuals[-2]
+##
+##print residuals
+##
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-18 08:06:05 UTC (rev 3445)
@@ -87,7 +87,7 @@
else:
#use user-defined aggregation
for AggOp in aggregation:
- P,candidates = sa_interpolation(A,candidates,omega=omega,AggOp=AggOp)
+ P,candidates,blocks = sa_interpolation(A,candidates,omega=omega,AggOp=AggOp)
A = (P.T.tocsr() * A) * P #galerkin operator
@@ -210,11 +210,11 @@
#ml = ruge_stuben_solver(A)
x = rand(A.shape[0])
- #b = zeros_like(x)
- b = A*rand(A.shape[0])
+ b = zeros_like(x)
+ #b = A*rand(A.shape[0])
if True:
- x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-12,return_residuals=True)
+ x_sol,residuals = ml.solve(b,x0=x,maxiter=30,tol=1e-8,return_residuals=True)
else:
residuals = []
def add_resid(x):
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/sa.py 2007-10-18 08:06:05 UTC (rev 3445)
@@ -2,28 +2,53 @@
import numpy
from numpy import array,arange,ones,zeros,sqrt,isinf,asarray,empty,diff,\
ascontiguousarray
-from scipy.sparse import csr_matrix,isspmatrix_csr
+from scipy.sparse import csr_matrix,isspmatrix_csr,spidentity
-from utils import diag_sparse,approximate_spectral_radius
+from utils import diag_sparse,approximate_spectral_radius,expand_into_blocks
import multigridtools
__all__ = ['sa_filtered_matrix','sa_strong_connections','sa_constant_interpolation',
'sa_interpolation','sa_smoothed_prolongator','sa_fit_candidates']
+## nnz = A.nnz
+##
+## indptr = A.indptr
+## indices = A.indices
+## data = A.data
+##
+## if n != 1:
+## # expand horizontally
+## indptr = n*A.indptr
+## indices = (n*indices).repeat(n) + tile(arange(n),nnz)
+##
+## if m != 1:
+## #expand vertically
+## indptr = concatenate( (array([0]), cumsum(diff(A).repeat(m))) )
+## indices = indices.repeat(m)
+##
+## if m != 1 or n != 1:
+## data = A.data.repeat(m*n)
+
+
def sa_filtered_matrix(A,epsilon,blocks=None):
+ """The filtered matrix is obtained from A by lumping all weak off-diagonal
+ entries onto the diagonal. Weak off-diagonals are determined by
+ the standard strength of connection measure using the parameter epsilon.
+
+ In the case epsilon = 0.0, (i.e. no weak connections) A is returned.
+ """
+
if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
if epsilon == 0:
A_filtered = A
-
else:
if blocks is None:
Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
A_filtered = csr_matrix((Sx,Sj,Sp),A.shape)
else:
A_filtered = A #TODO subtract weak blocks from diagonal blocks?
-
## num_dofs = A.shape[0]
## num_blocks = blocks.max() + 1
##
@@ -50,7 +75,6 @@
return A_filtered
-
def sa_strong_connections(A,epsilon):
if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
@@ -101,7 +125,9 @@
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))
+ #TODO fix this and add unittests
+ #AggOp = csr_matrix((AggOp.data.repeat(K),AggOp.indices.repeat(K),arange(K*N_fine + 1)),dims=(K*N_fine,N_coarse))
+ AggOp = expand_into_blocks(AggOp,K,1).tocsr()
N_fine = K*N_fine
#TODO convert this to list of coarse candidates
@@ -142,6 +168,20 @@
return Q,coarse_candidates
def sa_smoothed_prolongator(A,T,epsilon,omega,blocks=None):
+ """For a given matrix A and tentative prolongator T return the
+ smoothed prolongator P
+
+ P = (I - omega/rho(S) S) * T
+
+ where S is a Jacobi smoothing operator defined as follows:
+
+ omega - damping parameter
+ rho(S) - spectral radius of S (estimated)
+ S - inv(diag(A_filtered)) * A_filtered (Jacobi smoother)
+ A_filtered - sa_filtered_matrix(A,epsilon)
+ """
+
+
A_filtered = sa_filtered_matrix(A,epsilon,blocks) #use filtered matrix for anisotropic problems
D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
@@ -150,6 +190,9 @@
# smooth tentative prolongator T
P = T - (D_inv_A*T)
+
+ #S = (spidentity(A.shape[0]).tocsr() - D_inv_A) #TODO drop this?
+ #P = S * ( S * T)
return P
@@ -163,20 +206,15 @@
raise TypeError,'expected csr_matrix for argument AggOp'
if A.shape[1] != AggOp.shape[0]:
raise ValueError,'incompatible aggregation operator'
+
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)
-## D_inv = diag_sparse(1.0/diag_sparse(A_filtered))
-## D_inv_A = D_inv * A_filtered
-## D_inv_A *= omega/approximate_spectral_radius(D_inv_A)
-##
-## # smooth tentative prolongator T
-## P = T - (D_inv_A*T)
-
if blocks is not None:
blocks = arange(AggOp.shape[1]).repeat(len(candidates))
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-10-17 20:29:14 UTC (rev 3444)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-10-18 08:06:05 UTC (rev 3445)
@@ -1,9 +1,9 @@
__all__ =['approximate_spectral_radius','infinity_norm','diag_sparse',
- 'hstack_csr','vstack_csr']
+ 'hstack_csr','vstack_csr','expand_into_blocks']
import numpy
import scipy
-from numpy import ravel,arange,concatenate
+from numpy import ravel,arange,concatenate,tile
from scipy.linalg import norm
from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
csr_matrix,csc_matrix,extract_diagonal, \
@@ -67,3 +67,47 @@
V = concatenate((A.data,B.data))
return coo_matrix((V,(I,J)),dims=(A.shape[0]+B.shape[0],A.shape[1])).tocsr()
+
+def expand_into_blocks(A,m,n):
+ """Expand each element in a sparse matrix A into an m-by-n block.
+
+ Example:
+ >>> A.todense()
+ matrix([[ 1., 2.],
+ [ 4., 5.]])
+
+ >>> expand_into_blocks(A,2,2).todense()
+ matrix([[ 1., 1., 2., 2.],
+ [ 1., 1., 2., 2.],
+ [ 4., 4., 5., 5.],
+ [ 4., 4., 5., 5.]])
+
+ """
+ #TODO EXPLAIN MORE
+
+ if n is None:
+ n = m
+
+ if m == 1 and n == 1:
+ return A #nothing to do
+
+ A = A.tocoo()
+
+ # expand 1x1 -> mxn
+ row = ( m*A.row ).repeat(m*n).reshape(-1,m,n)
+ col = ( n*A.col ).repeat(m*n).reshape(-1,m,n)
+
+ # increment indices
+ row += tile(arange(m).reshape(-1,1),(1,n))
+ col += tile(arange(n).reshape(1,-1),(m,1))
+
+ # flatten
+ row = row.reshape(-1)
+ col = col.reshape(-1)
+
+ data = A.data.repeat(m*n)
+
+ return coo_matrix((data,(row,col)),dims=(m*A.shape[0],n*A.shape[1]))
+
+
+
More information about the Scipy-svn
mailing list