[Scipy-svn] r3458 - trunk/scipy/sandbox/multigrid
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Oct 23 00:08:26 EDT 2007
Author: wnbell
Date: 2007-10-22 23:08:10 -0500 (Mon, 22 Oct 2007)
New Revision: 3458
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/relaxation.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
continued work on adaptiveSA
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-23 04:08:10 UTC (rev 3458)
@@ -6,12 +6,66 @@
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
+from utils import approximate_spectral_radius,hstack_csr,vstack_csr,expand_into_blocks
+
+def augment_candidates(AggOp, old_Q, old_R, new_candidate, level):
+ K = old_R.shape[1]
+ #determine blocksizes
+ if level == 0:
+ old_bs = (1,K)
+ new_bs = (1,K+1)
+ else:
+ old_bs = (K,K)
+ new_bs = (K+1,K+1)
+
+ AggOp = expand_into_blocks(AggOp,new_bs[0],1).tocsr() #TODO switch to block matrix
+ # tentative prolongator
+ #TODO USE BSR
+ Q_indptr = (K+1)*AggOp.indptr
+ Q_indices = ((K+1)*AggOp.indices).repeat(K+1)
+ for i in range(K+1):
+ Q_indices[i::K+1] += i
+ Q_data = zeros((AggOp.shape[0]/new_bs[0],) + new_bs)
+ Q_data[:,:new_bs[0],:new_bs[1]] = old_Q.data.reshape((-1,) + old_bs) #TODO BSR change
+ # coarse candidates
+ R = zeros(((K+1)*AggOp.shape[1],K+1))
+ for i in range(K):
+ R[i::K+1,:K] = old_R[i::K,:]
+
+ c = new_candidate.reshape(-1)[diff(AggOp.indptr) == 1] #eliminate DOFs that aggregation misses
+ X = csr_matrix((c,AggOp.indices,AggOp.indptr),dims=AggOp.shape)
+
+ #orthogonalize X against previous
+ for i in range(K):
+ old_c = ascontiguousarray(Q_data[:,:,i].reshape(-1))
+ D_AtX = csr_matrix((old_c*X.data,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of A.T * X
+ R[i::K+1,-1] = D_AtX
+ X.data -= D_AtX[X.indices] * old_c
+
+ #normalize X
+ D_XtX = csr_matrix((X.data**2,X.indices,X.indptr),dims=X.shape).sum(axis=0).A.flatten() #same as diagonal of X.T * X
+ col_norms = sqrt(D_XtX)
+ R[K::K+1,-1] = col_norms
+
+ col_norms = 1.0/col_norms
+ col_norms[isinf(col_norms)] = 0
+ X.data *= col_norms[X.indices]
+ last_column = Q_data[:,:,-1].reshape(-1)
+ last_column = X.data.reshape(-1)
+ Q_data = Q_data.reshape(-1) #TODO BSR change
+
+ Q = csr_matrix((Q_data,Q_indices,Q_indptr),dims=(AggOp.shape[0],(K+1)*AggOp.shape[1]))
+
+ return Q,R
+
+
+
+
def orthonormalize_prolongator(P_l,x_l,W_l,W_m):
"""
@@ -72,26 +126,28 @@
Construct multilevel hierarchy using Smoothed Aggregation
Inputs:
A - matrix
- Is - list of constant prolongators
+ Ps - list of constant prolongators
B - "candidate" basis function to be approximated
Ouputs:
- (As,Is,Ps) - tuple of lists
- - As - [A, Ps[0].T*A*Ps[0], Ps[1].T*A*Ps[1], ... ]
- - Is - smoothed prolongators
- - Ps - tentative prolongators
+ (As,Ps,Ts) - tuple of lists
+ - As - [A, Ts[0].T*A*Ts[0], Ts[1].T*A*Ts[1], ... ]
+ - Ps - smoothed prolongators
+ - Ts - tentative prolongators
"""
As = [A]
- Is = []
Ps = []
+ Ts = []
+ Bs = [B]
for W in Ws:
P,B = sa_fit_candidates(W,B)
I = smoothed_prolongator(P,A)
A = I.T.tocsr() * A * I
As.append(A)
- Ps.append(P)
- Is.append(I)
- return As,Is,Ps
+ Ts.append(P)
+ Ps.append(I)
+ Bs.append(B)
+ return As,Ps,Ts,Bs
def make_bridge(I,N):
tail = I.indptr[-1].repeat(N - I.shape[0])
@@ -117,29 +173,30 @@
Ws = AggOps
- self.candidates = x
+ fine_candidates = x
#create SA using x here
- As,Is,Ps = sa_hierarchy(A,AggOps,self.candidates)
+ As,Ps,Ts,self.candidates = sa_hierarchy(A,AggOps,fine_candidates)
+ #self.candidates = [x]
for i in range(max_candidates - 1):
- #x = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps,mu=mu)
- x = self.__develop_candidate(As,Is,Ps,AggOps,self.candidates,mu=mu)
+ #x = self.__develop_candidate_OLD(A,As,Ps,Ts,Ws,AggOps,mu=mu)
+ #As,Ps,Ts,Ws = self.__augment_cycle(A,As,Ts,Ws,AggOps,x)
+ #self.candidates.append(x)
- self.candidates = hstack((self.candidates,x))
-
#TODO which is faster?
- #As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)
- As,Is,Ps = sa_hierarchy(A,AggOps,self.candidates)
+ x = self.__develop_candidate(As,Ps,Ts,AggOps,self.candidates,mu=mu)
+ fine_candidates = hstack((fine_candidates,x))
+ As,Ps,Ts,self.candidates = sa_hierarchy(A,AggOps,fine_candidates)
- self.Ps = Ps
- self.solver = multilevel_solver(As,Is)
+ self.Ts = Ts
+ self.solver = multilevel_solver(As,Ps)
self.AggOps = AggOps
def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
AggOps = []
- Is = []
+ Ps = []
# aSA parameters
# mu - number of test relaxation iterations
@@ -166,7 +223,7 @@
print 'A_l.shape',A_l.shape
AggOps.append(W_l)
- Is.append(I_l)
+ Ps.append(I_l)
As.append(A_l)
if A_l.shape <= max_coarse: break
@@ -183,7 +240,7 @@
x = x_hat #need to restore x
#update fine-level candidate
- for A_l,I in reversed(zip(As[1:],Is)):
+ for A_l,I in reversed(zip(As[1:],Ps)):
gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric') #TEST
x = I * x
gauss_seidel(A,x,zeros_like(x),iterations=mu,sweep='symmetric') #TEST
@@ -191,93 +248,113 @@
return x,AggOps #first candidate,aggregation
- def __develop_candidate(self,As,Is,Ps,AggOps,candidates,mu):
- #scipy.random.seed(0) #TEST
- x = scipy.rand(A.shape[0])
+ def __develop_candidate(self,As,Ps,Ts,AggOps,candidates,mu):
+ A = As[0]
+
+ x = A*scipy.rand(A.shape[0],1)
b = zeros_like(x)
- solver = multilevel_solver(As,Is)
-
- x = solver.solve(b, x0=x, tol=1e-10, maxiter=mu)
-
+ x = multilevel_solver(As,Ps).solve(b, x0=x, tol=1e-10, maxiter=mu)
+
#TEST FOR CONVERGENCE HERE
- #TODO augment candiates each time, then call fit_candidates?
+ temp_Ps = []
+ temp_As = [A]
+
+ def make_bridge(P,K):
+ indptr = P.indptr[:-1].reshape(-1,K-1)
+ indptr = hstack((indptr,indptr[:,-1].reshape(-1,1)))
+ indptr = indptr.reshape(-1)
+ indptr = hstack((indptr,indptr[-1:])) #duplicate last element
+ return csr_matrix((P.data,P.indices,indptr),dims=(K*P.shape[0]/(K-1),P.shape[1]))
- A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
-
- temp_Is = []
- for i in range(len(As) - 2):
- P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])
+ for i in range(len(As) - 2):
+ #TODO test augment_candidates against fit candidates
+ if i == 0:
+ temp = hstack((candidates[0],x))
+ else:
+ K = candidates[i].shape[1]
+ temp = zeros((x.shape[0]/(K+1),K+1,K + 1))
+ temp[:,:-1,:-1] = candidates[i].reshape(-1,K,K)
+ temp[:,:,-1] = x.reshape(-1,K+1,1)
+ temp = temp.reshape(-1,K+1)
+ T_,R_ = sa_fit_candidates(AggOps[i],temp)
+ #print "T - T_",(T - T_).data.max()
+ #assert((T - T_).data.max() < 1e-10)
+ #assert((R - R_).data.max() < 1e-10)
+ T,R = T_,R_
+ #TODO end test
+
+ #T,R = augment_candidates(AggOps[i], Ts[i], candidates[i+1], x, i)
+ P = smoothed_prolongator(T,A)
+ A = P.T.tocsr() * A * P
+
+ temp_Ps.append(P)
+ temp_As.append(A)
+
+ #TODO USE BSR (K,K) -> (K,K-1)
+ bridge = make_bridge(Ps[i+1],R.shape[1])
- I_l_new = smoothed_prolongator(P_l_new,A_l)
- A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
- bridge = make_bridge(Is[i+1],A_m_new.shape[0])
+ solver = multilevel_solver( [A] + As[i+2:], [bridge] + Ps[i+2:] )
+
+ x = R[:,-1].reshape(-1,1)
+ x = solver.solve(zeros_like(x), x0=x, tol=1e-8, maxiter=mu)
- temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Is[i+2:] )
-
- for n in range(mu):
- x_m = temp_solver.solve(zeros_like(x_m), x0=x_m, tol=1e-8, maxiter=1)
-
- temp_Is.append(I_l_new)
-
- W_l = vstack_csr(Ws[i+1],W_m_new) #prepare for next iteration
- A_l = A_m_new
- x_l = x_m
- P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
-
- x = x_l
- for I in reversed(temp_Is):
- x = I*x
+ for A,P in reversed(zip(temp_As,temp_Ps)):
+ x = P * x
+ gauss_seidel(A,x,zeros_like(x),iterations=mu,sweep='symmetric')
+ #for P in reversed(temp_Ps):
+ # x = P*x
+
return x
-## def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps,mu):
+## def __develop_candidate_OLD(self,A,As,Ps,Ts,Ws,AggOps,mu):
## #scipy.random.seed(0) #TEST
## x = scipy.rand(A.shape[0])
## b = zeros_like(x)
##
-## solver = multilevel_solver(As,Is)
+## solver = multilevel_solver(As,Ps)
##
## x = solver.solve(b, x0=x, tol=1e-10, maxiter=mu)
##
## #TEST FOR CONVERGENCE HERE
##
-## A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
+## A_l,P_l,W_l,x_l = As[0],Ts[0],Ws[0],x
##
-## temp_Is = []
+## temp_Ps = []
## for i in range(len(As) - 2):
## P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])
##
## I_l_new = smoothed_prolongator(P_l_new,A_l)
## A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
-## bridge = make_bridge(Is[i+1],A_m_new.shape[0])
+## bridge = make_bridge(Ps[i+1],A_m_new.shape[0])
##
-## temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Is[i+2:] )
+## temp_solver = multilevel_solver( [A_m_new] + As[i+2:], [bridge] + Ps[i+2:] )
##
## for n in range(mu):
## x_m = temp_solver.solve(zeros_like(x_m), x0=x_m, tol=1e-8, maxiter=1)
##
-## temp_Is.append(I_l_new)
+## temp_Ps.append(I_l_new)
##
## W_l = vstack_csr(Ws[i+1],W_m_new) #prepare for next iteration
## A_l = A_m_new
## x_l = x_m
-## P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
+## P_l = make_bridge(Ts[i+1],A_m_new.shape[0])
##
## x = x_l
-## for I in reversed(temp_Is):
+## for I in reversed(temp_Ps):
## x = I*x
##
## return x
- def __augment_cycle(self,A,As,Ps,Ws,AggOps,x):
+ def __augment_cycle(self,A,As,Ts,Ws,AggOps,x):
#make a new cycle using the new candidate
- A_l,P_l,W_l,x_l = As[0],Ps[0],AggOps[0],x
+ A_l,P_l,W_l,x_l = As[0],Ts[0],AggOps[0],x
- new_As,new_Is,new_Ps,new_Ws = [A],[],[],[AggOps[0]]
+ new_As,new_Ps,new_Ts,new_Ws = [A],[],[],[AggOps[0]]
for i in range(len(As) - 2):
P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])
@@ -288,24 +365,24 @@
new_As.append(A_m_new)
new_Ws.append(W_m_new)
- new_Is.append(I_l_new)
- new_Ps.append(P_l_new)
+ new_Ps.append(I_l_new)
+ new_Ts.append(P_l_new)
#prepare for next iteration
W_l = W_m_new
A_l = A_m_new
x_l = x_m
- P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
+ P_l = make_bridge(Ts[i+1],A_m_new.shape[0])
P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, csr_matrix((P_l.shape[1],1)))
I_l_new = smoothed_prolongator(P_l_new,A_l)
A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
new_As.append(A_m_new)
- new_Is.append(I_l_new)
- new_Ps.append(P_l_new)
+ new_Ps.append(I_l_new)
+ new_Ts.append(P_l_new)
- return new_As,new_Is,new_Ps,new_Ws
+ return new_As,new_Ps,new_Ts,new_Ws
@@ -325,10 +402,10 @@
#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()
-#blocks = arange(A.shape[0]/2).repeat(2)
+A = io.mmread("tests/sample_data/elas30_A.mtx").tocsr()
+blocks = arange(A.shape[0]/2).repeat(2)
-asa = adaptive_sa_solver(A,max_candidates=1,mu=10)
+asa = adaptive_sa_solver(A,max_candidates=4,mu=5)
scipy.random.seed(0) #make tests repeatable
x = rand(A.shape[0])
#b = A*rand(A.shape[0])
@@ -356,7 +433,7 @@
print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
-for c in asa.candidates.T:
+for c in asa.candidates[0].T:
print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
@@ -364,6 +441,19 @@
##W = asa.AggOps[0]*asa.AggOps[1]
##pcolor((W * rand(W.shape[1])).reshape((200,200)))
+def plot2d_arrows(x):
+ from pylab import quiver
+ x = x.reshape(-1)
+ N = (len(x)/2)**0.5
+ assert(2 * N * N == len(x))
+ X = linspace(-1,1,N).reshape(1,N).repeat(N,0).reshape(-1)
+ Y = linspace(-1,1,N).reshape(1,N).repeat(N,0).T.reshape(-1)
+
+ dX = x[0::2]
+ dY = x[1::2]
+
+ quiver(X,Y,dX,dY)
+
def plot2d(x):
from pylab import pcolor
pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-10-23 04:08:10 UTC (rev 3458)
@@ -168,7 +168,7 @@
if x0 is None:
x = zeros_like(b)
else:
- x = array(x0)
+ x = array(x0) #copy
#TODO change use of tol (relative tolerance) to agree with other iterative solvers
@@ -193,19 +193,21 @@
A = self.As[lvl]
if len(self.As) == 1:
- x[:] = spsolve(A,b)
+ #TODO make spsolve preserve dimensions
+ x[:] = spsolve(A,b).reshape(x.shape)
return
self.presmoother(A,x,b)
residual = b - A*x
- coarse_x = zeros((self.As[lvl+1].shape[0]))
coarse_b = self.Ps[lvl].T * residual
+ coarse_x = zeros_like(coarse_b)
if lvl == len(self.As) - 2:
#use direct solver on coarsest level
- coarse_x[:] = spsolve(self.As[-1],coarse_b) #TODO reuse factors for efficiency?
+ #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]
#print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
else:
Modified: trunk/scipy/sandbox/multigrid/relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/relaxation.py 2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/relaxation.py 2007-10-23 04:08:10 UTC (rev 3458)
@@ -77,6 +77,9 @@
iterations - number of iterations to perform (default: 1)
omega - damping parameter (default: 1.0)
"""
+ x = asarray(x).reshape(-1)
+ b = asarray(b).reshape(-1)
+
sweep = slice(None)
(row_start,row_stop,row_step) = sweep.indices(A.shape[0])
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-10-22 20:13:01 UTC (rev 3457)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-10-23 04:08:10 UTC (rev 3458)
@@ -65,8 +65,8 @@
if not isspmatrix(A) or not isspmatrix(B):
raise TypeError,'expected sparse matrix'
- if A.shape[0] != B.shape[0]:
- raise ValueError,'row dimensions must agree'
+ if A.shape[1] != B.shape[1]:
+ raise ValueError,'column dimensions must agree'
A = A.tocoo()
B = B.tocoo()
@@ -94,9 +94,6 @@
#TODO EXPLAIN MORE
#TODO use spkron instead, time for compairson
- if n is None:
- n = m
-
if m == 1 and n == 1:
return A #nothing to do
More information about the Scipy-svn
mailing list