[Scipy-svn] r3854 - in trunk/scipy: sandbox/multigrid sparse
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Jan 22 10:06:30 EST 2008
Author: wnbell
Date: 2008-01-22 09:06:22 -0600 (Tue, 22 Jan 2008)
New Revision: 3854
Modified:
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sparse/bsr.py
Log:
minor change to SA
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2008-01-21 20:40:25 UTC (rev 3853)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2008-01-22 15:06:22 UTC (rev 3854)
@@ -102,35 +102,45 @@
As = [A]
Ps = []
+ Rs = []
if aggregation is None:
while len(As) < max_levels and A.shape[0] > max_coarse:
P,B = sa_interpolation(A,B,epsilon*0.5**(len(As)-1),omega=omega)
+ R = P.T.asformat(P.format)
- A = (P.T.asformat(P.format) * A) * P #galerkin operator
+ A = R * A * P #galerkin operator
As.append(A)
+ Rs.append(R)
Ps.append(P)
else:
#use user-defined aggregation
for AggOp in aggregation:
P,B = sa_interpolation(A,B,omega=omega,AggOp=AggOp)
+ R = P.T.asformat(P.format)
- A = (P.T.tocsr() * A) * P #galerkin operator
+ A = R * A * P #galerkin operator
As.append(A)
+ Rs.append(R)
Ps.append(P)
- return multilevel_solver(As,Ps,preprocess=pre,postprocess=post)
+ return multilevel_solver(As,Ps,Rs=Rs,preprocess=pre,postprocess=post)
class multilevel_solver:
- def __init__(self,As,Ps,preprocess=None,postprocess=None):
+ def __init__(self,As,Ps,Rs=None,preprocess=None,postprocess=None):
self.As = As
self.Ps = Ps
self.preprocess = preprocess
self.postprocess = postprocess
+ if Rs is None:
+ self.Rs = [P.T for P in self.Ps]
+ else:
+ self.Rs = Rs
+
def __repr__(self):
output = 'multilevel_solver\n'
output += 'Number of Levels: %d\n' % len(self.As)
@@ -202,7 +212,7 @@
residual = b - A*x
- coarse_b = self.Ps[lvl].T * residual
+ coarse_b = self.Rs[lvl] * residual
coarse_x = zeros_like(coarse_b)
if lvl == len(self.As) - 2:
@@ -212,7 +222,7 @@
#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)
+ print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
else:
self.__solve(lvl+1,coarse_x,coarse_b)
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2008-01-21 20:40:25 UTC (rev 3853)
+++ trunk/scipy/sandbox/multigrid/sa.py 2008-01-22 15:06:22 UTC (rev 3854)
@@ -96,38 +96,38 @@
else:
sa_constant_interpolation(csr_matrix(A),epsilon)
-def sa_fit_candidates(AggOp,candidates,tol=1e-10):
+def sa_fit_candidates(AggOp,B,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')
+ if B.dtype != 'float32':
+ B = asarray(B,dtype='float64')
- if len(candidates.shape) != 2:
+ if len(B.shape) != 2:
raise ValueError,'expected rank 2 array for argument B'
- if candidates.shape[0] % AggOp.shape[0] != 0:
+ if B.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]
+ K = B.shape[1] # number of near-nullspace candidates
+ blocksize = B.shape[0] / AggOp.shape[0]
N_fine,N_coarse = AggOp.shape
- R = zeros((N_coarse,K,K),dtype=candidates.dtype) #storage for coarse candidates
+ R = zeros((N_coarse,K,K), dtype=B.dtype) #storage for coarse candidates
candidate_matrices = []
- threshold = tol * abs(candidates).max() # cutoff for small basis functions
-
for i in range(K):
- c = candidates[:,i]
+ c = B[:,i]
c = c.reshape(-1,blocksize,1)[diff(AggOp.indptr) == 1] # eliminate DOFs that aggregation misses
X = bsr_matrix( (c, AggOp.indices, AggOp.indptr), \
shape=(blocksize*N_fine, N_coarse) )
+ col_thresholds = tol * bsr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten()
+
#orthogonalize X against previous
for j,A in enumerate(candidate_matrices):
D_AtX = bsr_matrix((A.data*X.data,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
@@ -135,9 +135,10 @@
X.data -= scale_columns(A,D_AtX).data
#normalize X
- D_XtX = bsr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
- col_norms = sqrt(D_XtX)
- mask = col_norms < threshold # set small basis functions to 0
+ col_norms = bsr_matrix((X.data**2,X.indices,X.indptr),shape=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
+ mask = col_norms < col_thresholds # set small basis functions to 0
+
+ col_norms = sqrt(col_norms)
col_norms[mask] = 0
R[:,i,i] = col_norms
col_norms = 1.0/col_norms
@@ -186,7 +187,7 @@
return P
-def sa_interpolation(A,candidates,epsilon=0.0,omega=4.0/3.0,AggOp=None):
+def sa_interpolation(A,B,epsilon=0.0,omega=4.0/3.0,AggOp=None):
if not (isspmatrix_csr(A) or isspmatrix_bsr(A)):
A = csr_matrix(A)
@@ -200,7 +201,7 @@
raise ValueError,'incompatible aggregation operator'
- T,coarse_candidates = sa_fit_candidates(AggOp,candidates)
+ T,coarse_candidates = sa_fit_candidates(AggOp,B)
P = sa_smoothed_prolongator(A,T,epsilon,omega)
return P,coarse_candidates
Modified: trunk/scipy/sparse/bsr.py
===================================================================
--- trunk/scipy/sparse/bsr.py 2008-01-21 20:40:25 UTC (rev 3853)
+++ trunk/scipy/sparse/bsr.py 2008-01-22 15:06:22 UTC (rev 3854)
@@ -166,7 +166,7 @@
else:
self.shape = shape
- self.check_format()
+ self.check_format(full_check=False)
def check_format(self, full_check=True):
"""check whether the matrix format is valid
More information about the Scipy-svn
mailing list