[Scipy-svn] r3307 - in trunk/scipy/sandbox/multigrid: . multigridtools tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Sep 7 15:56:01 EDT 2007
Author: wnbell
Date: 2007-09-07 14:55:57 -0500 (Fri, 07 Sep 2007)
New Revision: 3307
Added:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
trunk/scipy/sandbox/multigrid/tests/test_utils.py
Modified:
trunk/scipy/sandbox/multigrid/coarsen.py
trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
trunk/scipy/sandbox/multigrid/multilevel.py
trunk/scipy/sandbox/multigrid/relaxation.py
trunk/scipy/sandbox/multigrid/tests/test_relaxation.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
added test cases for utils and coarsening
added preliminary adaptive SA code
Added: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -0,0 +1,371 @@
+import numpy,scipy,scipy.sparse
+from numpy import sqrt,ravel,diff,zeros,zeros_like,inner,concatenate
+from scipy.sparse import csr_matrix,coo_matrix
+
+from relaxation import gauss_seidel
+from multilevel import multilevel_solver
+from coarsen import sa_constant_interpolation
+from utils import infinity_norm
+
+
+def fit_candidate(I,x):
+ """
+ For each aggregate in I (i.e. each column of I) compute vector R and
+ sparse matrix Q (having the sparsity of I) such that the following holds:
+
+ Q*R = x and Q^T*Q = I
+
+ In otherwords, find a prolongator Q with orthonormal columns so that
+ x is represented exactly on the coarser level by R.
+ """
+ Q = csr_matrix((x.copy(),I.indices,I.indptr),dims=I.shape,check=False)
+ R = sqrt(ravel(csr_matrix((x*x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0))) #column 2-norms
+ Q.data *= (1.0/R)[Q.indices]
+
+ #print "norm(R)",scipy.linalg.norm(R)
+ #print "min(R),max(R)",min(R),max(R)
+ #print "infinity_norm(Q.T*Q - I) ",infinity_norm((Q.T.tocsr() * Q - scipy.sparse.spidentity(Q.shape[1])))
+ #print "norm(Q*R - x)",scipy.linalg.norm(Q*R - x)
+ #print "norm(x - Q*Q.Tx)",scipy.linalg.norm(x - Q*(Q.T*x))
+ return Q,R
+
+
+##def orthonormalize_candidate(I,x,basis):
+## Px = csr_matrix((x,I.indices,I.indptr),dims=I.shape,check=False)
+## Rs = []
+## #othogonalize columns of Px against other candidates
+## for b in basis:
+## Pb = csr_matrix((b,I.indices,I.indptr),dims=I.shape,check=False)
+## R = ravel(csr_matrix((Pb.data*Px.data,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0)) # columnwise projection of Px on Pb
+## Px.data -= R[I.indices] * Pb.data #subtract component in b direction
+## Rs.append(R)
+##
+## #filter columns here, set unused cols to 0, add to mask
+##
+## #normalize columns of Px
+## R = ravel(csr_matrix((x**x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0))
+## Px.data *= (1.0/R)[I.indices]
+## Rs.append(R.reshape(-1,1))
+## return Rs
+
+def hstack_csr(A,B):
+ #OPTIMIZE THIS
+ assert(A.shape[0] == B.shape[0])
+ A = A.tocoo()
+ B = B.tocoo()
+ I = concatenate((A.row,B.row))
+ J = concatenate((A.col,B.col+A.shape[1]))
+ V = concatenate((A.data,B.data))
+ return coo_matrix((V,(I,J)),dims=(A.shape[0],A.shape[1]+B.shape[1])).tocsr()
+
+
+def vstack_csr(A,B):
+ #OPTIMIZE THIS
+ assert(A.shape[1] == B.shape[1])
+ A = A.tocoo()
+ B = B.tocoo()
+ I = concatenate((A.row,B.row+A.shape[0]))
+ J = concatenate((A.col,B.col))
+ V = concatenate((A.data,B.data))
+ return coo_matrix((V,(I,J)),dims=(A.shape[0]+B.shape[0],A.shape[1])).tocsr()
+
+
+
+def orthonormalize_prolongator(P_l,x_l,W_l,W_m):
+ """
+
+ """
+ X = csr_matrix((x_l,W_l.indices,W_l.indptr),dims=W_l.shape,check=False) #candidate prolongator (assumes every value from x is used)
+
+ R = (P_l.T.tocsr() * X) # R has at most 1 nz per row
+ X = X - P_l*R # othogonalize X against P_l
+
+ #DROP REDUNDANT COLUMNS FROM P (AND R?) HERE (NULL OUT R ACCORDINGLY?)
+ #REMOVE CORRESPONDING COLUMNS FROM W_l AND ROWS FROM A_m ALSO
+ W_l_new = W_l
+ W_m_new = W_m
+
+ #normalize surviving columns of X
+ col_norms = ravel(sqrt(csr_matrix((X.data*X.data,X.indices,X.indptr),dims=X.shape,check=False).sum(axis=0)))
+ print "zero cols",sum(col_norms == 0)
+ print "small cols",sum(col_norms < 1e-8)
+ Xcopy = X.copy()
+ X.data *= (1.0/col_norms)[X.indices]
+
+ P_l_new = hstack_csr(P_l,X)
+
+
+ #check orthonormality
+ print "norm(P.T*P - I) ",scipy.linalg.norm((P_l_new.T * P_l_new - scipy.sparse.spidentity(P_l_new.shape[1])).data)
+ #assert(scipy.linalg.norm((P_l_new.T * P_l_new - scipy.sparse.spidentity(P_l_new.shape[1])).data)<1e-8)
+
+ x_m = zeros(P_l_new.shape[1],dtype=x_l.dtype)
+ x_m[:P_l.shape[1]][diff(R.indptr).astype('bool')] = R.data
+ x_m[P_l.shape[1]:] = col_norms
+
+ print "||x_l - P_l*x_m||",scipy.linalg.norm(P_l_new* x_m - x_l) #see if x_l is represented exactly
+
+ return P_l_new,x_m,W_l,W_m
+
+
+def smoothed_prolongator(P,A):
+ #just use Richardson for now
+ #omega = 4.0/(3.0*infinity_norm(A))
+ #return P - omega*(A*P)
+ #return P
+ D = diag_sparse(A)
+ D_inv_A = diag_sparse(1.0/D)*A
+ omega = 4.0/(3.0*infinity_norm(D_inv_A))
+ D_inv_A *= omega
+ return P - D_inv_A*P
+
+
+def sa_hierarchy(A,Ws,x):
+ """
+ Construct multilevel hierarchy using Smoothed Aggregation
+ Inputs:
+ A - matrix
+ Is - list of constant prolongators
+ x - "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 = [A]
+ Is = []
+ Ps = []
+
+ for W in Ws:
+ P,x = fit_candidate(W,x)
+ I = smoothed_prolongator(P,A)
+ A = I.T.tocsr() * A * I
+ As.append(A)
+ Ps.append(P)
+ Is.append(I)
+ return As,Is,Ps
+
+def make_bridge(I,N):
+ tail = I.indptr[-1].repeat(N - I.shape[0])
+ ptr = concatenate((I.indptr,tail))
+ return csr_matrix((I.data,I.indices,ptr),dims=(N,I.shape[1]),check=False)
+
+class adaptive_sa_solver:
+ def __init__(self,A,options=None):
+ self.A = A
+
+ self.Rs = []
+ self.__construct_hierarchy(A)
+
+ def __construct_hierarchy(self,A):
+ #if self.A.shape[0] <= self.opts['coarse: max size']:
+ # raise ValueError,'small matrices not handled yet'
+
+ x,AggOps = self.__initialization_stage(A) #first candidate
+ Ws = AggOps
+
+ #x[:] = 1 #TEST
+
+ self.candidates = [x]
+ #self.candidates = [1.0/D.data]
+
+ #create SA using x here
+ As,Is,Ps = sa_hierarchy(A,Ws,x)
+
+ for i in range(0):
+ x = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps)
+ #if i == 0:
+ # x = arange(20).repeat(20).astype(float)
+ #elif i == 1:
+ # x = arange(20).repeat(20).astype(float)
+ # x = numpy.ravel(transpose(x.reshape((20,20))))
+
+ As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)
+
+ self.candidates.append(x)
+
+ self.Ps = Ps
+ self.solver = multilevel_solver(As,Is)
+ self.AggOps = AggOps
+
+
+
+ def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps):
+ x = scipy.rand(A.shape[0])
+ b = zeros_like(x)
+
+
+ #x = arange(200).repeat(200).astype(float)
+ #x[:] = 1 #TEST
+
+ mu = 5
+
+ solver = multilevel_solver(As,Is)
+
+ for n in range(mu):
+ x = solver.solve(b, x0=x, tol=1e-8, maxiter=1)
+ #TEST FOR CONVERGENCE HERE
+
+ 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])
+
+ 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])
+
+ 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
+
+ return x
+
+
+ def __augment_cycle(self,A,As,Ps,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
+
+ new_As,new_Is,new_Ps,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])
+
+ I_l_new = smoothed_prolongator(P_l_new,A_l)
+ A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
+ W_m_new = vstack_csr(Ws[i+1],W_m_new)
+
+ new_As.append(A_m_new)
+ new_Ws.append(W_m_new)
+ new_Is.append(I_l_new)
+ new_Ps.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_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)
+
+ return new_As,new_Is,new_Ps,new_Ws
+
+
+ def __initialization_stage(self,A):
+ max_levels = 10
+ max_coarse = 50
+
+ AggOps = []
+ Is = []
+
+ # aSA parameters
+ mu = 5 # number of test relaxation iterations
+ epsilon = 0.1 # minimum acceptable relaxation convergence factor
+
+ scipy.random.seed(0)
+
+ #step 1
+ A_l = A
+ x = scipy.rand(A_l.shape[0])
+ skip_f_to_i = False
+
+ #step 2
+ b = zeros_like(x)
+ gauss_seidel(A_l,x,b,iterations=mu)
+ #step 3
+ #test convergence rate here
+
+ As = [A]
+
+ while len(AggOps) + 1 < max_levels and A_l.shape[0] > max_coarse:
+ W_l = sa_constant_interpolation(A_l,epsilon=0.08*0.5**(len(AggOps)-1)) #step 4b #TEST
+ #W_l = sa_constant_interpolation(A_l,epsilon=0) #step 4b #TEST
+ P_l,x = fit_candidate(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
+
+ AggOps.append(W_l)
+ Is.append(I_l)
+ As.append(A_l)
+
+ 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) #step 4h
+ x_A_x = inner(x,A_l*x)
+ if (x_A_x/inner(x_hat,A_l*x_hat))**(1.0/mu) < epsilon: #step 4i
+ print "sufficient convergence, skipping"
+ skip_f_to_i = True
+ if x_A_x == 0:
+ x = x_hat #need to restore x
+
+ #update fine-level candidate
+ for A_l,I in reversed(zip(As[1:],Is)):
+ gauss_seidel(A_l,x,zeros_like(x),iterations=mu) #TEST
+ x = I * x
+
+ gauss_seidel(A,x,b,iterations=mu) #TEST
+
+ return x,AggOps #first candidate,aggregation
+
+
+
+
+from scipy import *
+from utils import diag_sparse
+from multilevel import poisson_problem1D,poisson_problem2D
+#A = poisson_problem2D(100)
+A = io.mmread("tests/sample_data/laplacian_40_3dcube.mtx").tocsr()
+
+#A = A*A
+#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+#A = D * A * D
+#A = io.mmread("nos2.mtx").tocsr()
+asa = adaptive_sa_solver(A)
+x = rand(A.shape[0])
+b = zeros_like(x)
+
+
+print "solving"
+x_sol,residuals = asa.solver.solve(b,x,tol=1e-12,maxiter=30,return_residuals=True)
+residuals = array(residuals)/residuals[0]
+print "residuals ",residuals
+
+print asa.solver
+
+print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
+
+for c in asa.candidates:
+ 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)))
+
+def plot2d(x):
+ from pylab import pcolor
+ pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
+
Modified: trunk/scipy/sandbox/multigrid/coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/coarsen.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/coarsen.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -1,10 +1,9 @@
-from scipy import *
import multigridtools
import scipy
import numpy
-from utils import diag_sparse,inf_norm
+from utils import diag_sparse,infinity_norm
def rs_strong_connections(A,theta):
@@ -33,37 +32,45 @@
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+
return scipy.sparse.csr_matrix((Sx,Sj,Sp),A.shape)
-def sa_constant_interpolation(A,epsilon=None):
+def sa_constant_interpolation(A,epsilon):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
- if epsilon is not None:
- S = sa_strong_connections(A,epsilon)
- else:
- S = A
-
+ S = sa_strong_connections(A,epsilon)
+
+ #S.ensure_sorted_indices()
+
#tentative (non-smooth) interpolation operator I
- Ij = multigridtools.sa_get_aggregates(A.shape[0],S.indptr,S.indices)
- Ip = numpy.arange(len(Ij)+1)
- Ix = numpy.ones(len(Ij))
+ Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+ Pp = numpy.arange(len(Pj)+1)
+ Px = numpy.ones(len(Pj))
- return scipy.sparse.csr_matrix((Ix,Ij,Ip))
+ return scipy.sparse.csr_matrix((Px,Pj,Pp))
-
+##def sa_smoother(A,S,omega):
+## Bp,Bj,Bx = multigridtools.sa_smoother(A.shape[0],omega,A.indptr,A.indices,A.data,S.indptr,S.indices,S.data)
+##
+## return csr_matrix((Bx,Bj,Bp),dims=A.shape)
+
def sa_interpolation(A,epsilon,omega=4.0/3.0):
if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
+
+ P = sa_constant_interpolation(A,epsilon)
+
+## As = sa_strong_connections(A,epsilon)
+## S = sa_smoother(A,S,omega)
- I = sa_constant_interpolation(A,epsilon)
D_inv = diag_sparse(1.0/diag_sparse(A))
D_inv_A = D_inv * A
- D_inv_A *= -omega/inf_norm(D_inv_A)
+ D_inv_A *= omega/infinity_norm(D_inv_A)
- P = I + (D_inv_A*I) #same as P=S*I, (faster?)
-
- return P
+ I = P - (D_inv_A*P) #same as I=S*P, (faster?)
+
+ return I
Modified: trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h 2007-09-07 19:55:57 UTC (rev 3307)
@@ -5,8 +5,8 @@
#include <vector>
#include <iterator>
#include <assert.h>
+#include <cmath>
-
//#define DEBUG
@@ -20,24 +20,22 @@
Sp->push_back(0);
//compute diagonal values
- std::vector<T> diags(n_row);
+ std::vector<T> diags(n_row,T(0));
for(int i = 0; i < n_row; i++){
int row_start = Ap[i];
int row_end = Ap[i+1];
for(int jj = row_start; jj < row_end; jj++){
- if(Aj[jj] == i){
- diags[i] = Ax[jj];
- break;
- }
+ if(Aj[jj] == i){
+ diags[i] = Ax[jj];
+ break;
+ }
}
}
#ifdef DEBUG
for(int i = 0; i < n_row; i++){ assert(diags[i] > 0); }
#endif
-
-
for(int i = 0; i < n_row; i++){
int row_start = Ap[i];
int row_end = Ap[i+1];
@@ -45,14 +43,15 @@
T eps_Aii = epsilon*epsilon*diags[i];
for(int jj = row_start; jj < row_end; jj++){
- const int& j = Aj[jj];
- const T& Aij = Ax[jj];
+ const int j = Aj[jj];
+ const T Aij = Ax[jj];
if(i == j){continue;}
- if(Aij*Aij >= eps_Aii * diags[j]){
- Sj->push_back(j);
- Sx->push_back(Aij);
+ // |A(i,j)| < epsilon * sqrt(|A(i,i)|*|A(j,j)|)
+ if(Aij*Aij >= std::abs(eps_Aii * diags[j])){
+ Sj->push_back(j);
+ Sx->push_back(Aij);
}
}
Sp->push_back(Sj->size());
@@ -61,9 +60,9 @@
void sa_get_aggregates(const int n_row,
- const int Ap[], const int Aj[],
- std::vector<int> * Bj){
-
+ const int Ap[], const int Aj[],
+ std::vector<int> * Bj)
+{
std::vector<int> aggregates(n_row,-1);
int num_aggregates = 0;
@@ -72,21 +71,19 @@
for(int i = 0; i < n_row; i++){
if(aggregates[i] >= 0){ continue; } //already marked
- const int& row_start = Ap[i];
- const int& row_end = Ap[i+1];
+ const int row_start = Ap[i];
+ const int row_end = Ap[i+1];
-
//Determine whether all neighbors of this node are free (not already aggregates)
bool free_neighborhood = true;
for(int jj = row_start; jj < row_end; jj++){
if(aggregates[Aj[jj]] >= 0){
- free_neighborhood = false;
- break;
+ free_neighborhood = false;
+ break;
}
}
if(!free_neighborhood){ continue; } //bail out
-
//Make an aggregate out of this node and its strong neigbors
aggregates[i] = num_aggregates;
for(int jj = row_start; jj < row_end; jj++){
@@ -96,52 +93,49 @@
}
-
//Pass #2
std::vector<int> aggregates_copy(aggregates);
for(int i = 0; i < n_row; i++){
if(aggregates[i] >= 0){ continue; } //already marked
- const int& row_start = Ap[i];
- const int& row_end = Ap[i+1];
+ const int row_start = Ap[i];
+ const int row_end = Ap[i+1];
for(int jj = row_start; jj < row_end; jj++){
- const int& j = Aj[jj];
+ const int j = Aj[jj];
if(aggregates_copy[j] >= 0){
- aggregates[i] = aggregates_copy[j];
- break;
+ aggregates[i] = aggregates_copy[j];
+ break;
}
}
}
-
//Pass #3
for(int i = 0; i < n_row; i++){
if(aggregates[i] >= 0){ continue; } //already marked
-
- const int& row_start = Ap[i];
- const int& row_end = Ap[i+1];
+
+ const int row_start = Ap[i];
+ const int row_end = Ap[i+1];
aggregates[i] = num_aggregates;
for(int jj = row_start; jj < row_end; jj++){
- const int& j = Aj[jj];
+ const int j = Aj[jj];
if(aggregates[j] < 0){ //unmarked neighbors
- aggregates[j] = num_aggregates;
+ aggregates[j] = num_aggregates;
}
}
num_aggregates++;
}
-
#ifdef DEBUG
for(int i = 0; i < n_row; i++){ assert(aggregates[i] >= 0 && aggregates[i] < num_aggregates); }
#endif
- *Bj = aggregates;
+ aggregates.swap(*Bj);
}
Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/multilevel.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -10,9 +10,9 @@
from coarsen import sa_interpolation,rs_interpolation
from relaxation import gauss_seidel,jacobi
+from utils import infinity_norm
-
def poisson_problem1D(N):
"""
Return a sparse CSR matrix for the 1d poisson problem
@@ -21,7 +21,7 @@
"""
D = 2*numpy.ones(N)
O = -numpy.ones(N)
- return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocsr()
+ return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate zeros
def poisson_problem2D(N):
"""
@@ -33,7 +33,8 @@
T = -numpy.ones(N*N)
O = -numpy.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).tocsr()
+ return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate zeros
+
def ruge_stuben_solver(A,max_levels=10,max_coarse=500):
"""
@@ -41,8 +42,8 @@
References:
"Multigrid"
- Trottenberg, U., C. W. Oosterlee, and Anton Schuller. San Diego: Academic Press, 2001.
- See Appendix A
+ Trottenberg, U., C. W. Oosterlee, and Anton Schuller. San Diego: Academic Press, 2001.
+ See Appendix A
"""
As = [A]
@@ -58,7 +59,7 @@
return multilevel_solver(As,Ps)
-def smoothed_aggregation_solver(A,max_levels=10,max_coarse=500):
+def smoothed_aggregation_solver(A,max_levels=10,max_coarse=500,epsilon=0.08):
"""
Create a multilevel solver using Smoothed Aggregation (SA)
@@ -72,8 +73,9 @@
Ps = []
while len(As) < max_levels and A.shape[0] > max_coarse:
- P = sa_interpolation(A,epsilon=0.08*0.5**(len(As)-1))
-
+ P = sa_interpolation(A,epsilon=epsilon*0.5**(len(As)-1))
+ #P = sa_interpolation(A,epsilon=0.0)
+
A = (P.T.tocsr() * A) * P #galerkin operator
As.append(A)
@@ -108,7 +110,10 @@
"""number of unknowns on all levels / number of unknowns on the finest level"""
return sum([A.shape[0] for A in self.As])/float(self.As[0].shape[0])
-
+
+ def psolve(self, b):
+ return self.solve(b,maxiter=1)
+
def solve(self, b, x0=None, tol=1e-5, maxiter=100, callback=None, return_residuals=False):
"""
TODO
@@ -122,12 +127,12 @@
#TODO change use of tol (relative tolerance) to agree with other iterative solvers
A = self.As[0]
- residuals = [norm(b-A*x,2)]
+ residuals = [scipy.linalg.norm(b-A*x)]
while len(residuals) <= maxiter and residuals[-1]/residuals[0] > tol:
self.__solve(0,x,b)
- residuals.append(scipy.linalg.norm(b-A*x,2))
+ residuals.append(scipy.linalg.norm(b-A*x))
if callback is not None:
callback(x)
@@ -142,11 +147,11 @@
A = self.As[lvl]
if len(self.As) == 1:
- x[:] = scipy.linalg.solve(A.todense(),b)
- return x
+ x[:] = scipy.linsolve.spsolve(A,b)
+ return
self.presmoother(A,x,b)
-
+
residual = b - A*x
coarse_x = zeros((self.As[lvl+1].shape[0]))
@@ -154,7 +159,9 @@
if lvl == len(self.As) - 2:
#direct solver on coarsest level
- coarse_x[:] = scipy.linalg.solve(self.As[-1].todense(),coarse_b)
+ coarse_x[:] = scipy.linsolve.spsolve(self.As[-1],coarse_b)
+ #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:
self.__solve(lvl+1,coarse_x,coarse_b)
@@ -165,26 +172,32 @@
def presmoother(self,A,x,b):
gauss_seidel(A,x,b,iterations=1,sweep="forward")
+ #x += 4.0/(3.0*infinity_norm(A))*(b - A*x)
def postsmoother(self,A,x,b):
- gauss_seidel(A,x,b,iterations=1,sweep="backward")
+ gauss_seidel(A,x,b,iterations=1,sweep="forward")
+ #gauss_seidel(A,x,b,iterations=1,sweep="backward")
+ #x += 4.0/(3.0*infinity_norm(A))*(b - A*x)
if __name__ == '__main__':
from scipy import *
A = poisson_problem2D(200)
+ #A = io.mmread("rocker_arm_surface.mtx").tocsr()
+
ml = smoothed_aggregation_solver(A)
#ml = ruge_stuben_solver(A)
+
x = rand(A.shape[0])
b = zeros_like(x)
+ #b = rand(A.shape[0])
- resid = []
-
- for n in range(10):
- x = ml.solve(b,x,maxiter=1)
- resid.append(linalg.norm(A*x))
+ x_sol,residuals = ml.solve(b,x0=x,maxiter=40,tol=1e-10,return_residuals=True)
+ residuals = array(residuals)/residuals[0]
+ print residuals
+
Modified: trunk/scipy/sandbox/multigrid/relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/relaxation.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/relaxation.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -13,6 +13,12 @@
iterations - number of iterations to perform (default: 1)
sweep - slice of unknowns to relax (default: all in forward direction)
"""
+ if A.shape[0] != A.shape[1]:
+ raise ValueError,'expected symmetric matrix'
+
+ if A.shape[1] != len(x) or len(x) != len(b):
+ raise ValueError,'unexpected number of unknowns'
+
if sweep == 'forward':
row_start,row_stop,row_step = 0,len(x),1
elif sweep == 'backward':
Added: trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_coarsen.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/tests/test_coarsen.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -0,0 +1,160 @@
+from numpy.testing import *
+
+from numpy import sqrt,empty,ones,arange,array_split
+from scipy import rand
+from scipy.sparse import spdiags,csr_matrix,lil_matrix
+import numpy
+
+set_package_path()
+import scipy.multigrid
+from scipy.multigrid.coarsen import sa_strong_connections,sa_constant_interpolation
+from scipy.multigrid.multilevel import poisson_problem1D,poisson_problem2D
+restore_path()
+
+
+def reference_sa_strong_connections(A,epsilon):
+ A_coo = A.tocoo()
+ S = lil_matrix(A.shape)
+
+ for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
+ if i == j: continue #skip diagonal
+
+ if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
+ S[i,j] = v
+
+ return S.tocsr()
+
+
+# note that this method only tests the current implementation, not
+# all possible implementations
+def reference_sa_constant_interpolation(A,epsilon):
+ S = sa_strong_connections(A,epsilon)
+ S = array_split(S.indices,S.indptr[1:-1])
+
+ n = A.shape[0]
+
+ R = set(range(n))
+ j = 0
+
+ aggregates = empty(n,dtype=A.indices.dtype)
+ aggregates[:] = -1
+
+
+ # Pass #1
+ for i,row in enumerate(S):
+ Ni = set(row) | set([i])
+
+ if Ni.issubset(R):
+ R -= Ni
+ for x in Ni:
+ aggregates[x] = j
+ j += 1
+
+ # Pass #2
+ Old_R = R.copy()
+ for i,row in enumerate(S):
+ if i not in R: continue
+
+ for x in row:
+ if x not in Old_R:
+ aggregates[i] = aggregates[x]
+ R.remove(i)
+ break
+
+
+ # Pass #3
+ for i,row in enumerate(S):
+ if i not in R: continue
+ Ni = set(row) | set([i])
+
+ for x in Ni:
+ if x in R:
+ aggregates[x] = j
+ j += 1
+
+ assert(len(R) == 0)
+
+ Pj = aggregates
+ Pp = arange(n+1)
+ Px = ones(n)
+
+ return csr_matrix((Px,Pj,Pp))
+
+class test_sa_strong_connections(NumpyTestCase):
+ def check_simple(self):
+ N = 4
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+ S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+ assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense()) #all connections are strong
+ assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+
+ N = 100
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+ S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+ assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense()) #all connections are strong
+ assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+
+ def check_random(self):
+ numpy.random.seed(0)
+
+ for N in [2,3,5,10]:
+ A = csr_matrix(rand(N,N))
+ for epsilon in [0.0,0.1,0.5,0.8,1.0,10.0]:
+ S_result = sa_strong_connections(A,epsilon)
+ S_expected = reference_sa_strong_connections(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson1D(self):
+ for N in [2,3,5,7,10,11,19]:
+ A = poisson_problem1D(N)
+ for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+ S_result = sa_strong_connections(A,epsilon)
+ S_expected = reference_sa_strong_connections(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson2D(self):
+ for N in [2,3,5,7,10,11,19]:
+ A = poisson_problem2D(N)
+ for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+ S_result = sa_strong_connections(A,epsilon)
+ S_expected = reference_sa_strong_connections(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+## def check_sample_data(self):
+## for filename in all_matrices:
+## A = open_matrix(filename)
+
+
+S_result = None
+S_expected = None
+class test_sa_constant_interpolation(NumpyTestCase):
+ def check_random(self):
+ numpy.random.seed(0)
+
+ for N in [2,3,5,10]:
+ A = csr_matrix(rand(N,N))
+ for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+ S_result = sa_constant_interpolation(A,epsilon)
+ S_expected = reference_sa_constant_interpolation(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson1D(self):
+ for N in [2,3,5,7,10,11,20,21,29,30]:
+ A = poisson_problem1D(N)
+ for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+ S_result = sa_constant_interpolation(A,epsilon)
+ S_expected = reference_sa_constant_interpolation(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+ def check_poisson2D(self):
+ for N in [2,3,5,7,10,11,20,21,29,30]:
+ A = poisson_problem2D(N)
+ for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+ S_result = sa_constant_interpolation(A,epsilon)
+ S_expected = reference_sa_constant_interpolation(A,epsilon)
+ assert_array_equal(S_result.todense(),S_expected.todense())
+
+
+if __name__ == '__main__':
+ NumpyTest().run()
+
Modified: trunk/scipy/sandbox/multigrid/tests/test_relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_relaxation.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/tests/test_relaxation.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -7,6 +7,7 @@
set_package_path()
+import scipy.multigrid
from scipy.multigrid.relaxation import polynomial_smoother,gauss_seidel,jacobi
restore_path()
@@ -36,7 +37,51 @@
polynomial_smoother(A,x,b,[-0.14285714, 1., -2.])
assert_almost_equal(x,x0 - 0.14285714*A*A*r + A*r - 2*r)
+ def check_jacobi(self):
+ N = 1
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+ x = arange(N).astype(numpy.float64)
+ b = zeros(N)
+ jacobi(A,x,b)
+ assert_almost_equal(x,array([0]))
+ N = 3
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+ x = zeros(N)
+ b = arange(N).astype(numpy.float64)
+ jacobi(A,x,b)
+ assert_almost_equal(x,array([0.0,0.5,1.0]))
+
+ N = 3
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+ x = arange(N).astype(numpy.float64)
+ b = zeros(N)
+ jacobi(A,x,b)
+ assert_almost_equal(x,array([0.5,1.0,0.5]))
+
+ N = 1
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+ x = arange(N).astype(numpy.float64)
+ b = array([10])
+ jacobi(A,x,b)
+ assert_almost_equal(x,array([5]))
+
+ N = 3
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+ x = arange(N).astype(numpy.float64)
+ b = array([10,20,30])
+ jacobi(A,x,b)
+ assert_almost_equal(x,array([5.5,11.0,15.5]))
+
+ N = 3
+ A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+ x = arange(N).astype(numpy.float64)
+ x_copy = x.copy()
+ b = array([10,20,30])
+ jacobi(A,x,b,omega=1.0/3.0)
+ assert_almost_equal(x,2.0/3.0*x_copy + 1.0/3.0*array([5.5,11.0,15.5]))
+
+
def check_gauss_seidel(self):
N = 1
A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
Added: trunk/scipy/sandbox/multigrid/tests/test_utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_utils.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/tests/test_utils.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -0,0 +1,60 @@
+from numpy.testing import *
+
+import numpy
+import scipy
+from scipy import matrix,array,diag
+from scipy.sparse import csr_matrix
+
+
+set_package_path()
+from scipy.multigrid.utils import infinity_norm,diag_sparse
+restore_path()
+
+
+class test_utils(NumpyTestCase):
+ def check_infinity_norm(self):
+ A = matrix([[-4]])
+ assert_equal(infinity_norm(csr_matrix(A)),4)
+
+ A = matrix([[1,0,-5],[-2,5,0]])
+ assert_equal(infinity_norm(csr_matrix(A)),7)
+
+ A = matrix([[0,1],[0,-5]])
+ assert_equal(infinity_norm(csr_matrix(A)),5)
+
+ A = matrix([[1.3,-4.7,0],[-2.23,5.5,0],[9,0,-2]])
+ assert_equal(infinity_norm(csr_matrix(A)),11)
+
+ def check_diag_sparse(self):
+ #check sparse -> array
+ A = matrix([[-4]])
+ assert_equal(diag_sparse(csr_matrix(A)),[-4])
+
+ A = matrix([[1,0,-5],[-2,5,0]])
+ assert_equal(diag_sparse(csr_matrix(A)),[1,5])
+
+ A = matrix([[0,1],[0,-5]])
+ assert_equal(diag_sparse(csr_matrix(A)),[0,-5])
+
+ A = matrix([[1.3,-4.7,0],[-2.23,5.5,0],[9,0,-2]])
+ assert_equal(diag_sparse(csr_matrix(A)),[1.3,5.5,-2])
+
+ #check array -> sparse
+ A = matrix([[-4]])
+ assert_equal(diag_sparse(array([-4])).todense(),csr_matrix(A).todense())
+
+ A = matrix([[1,0],[0,5]])
+ assert_equal(diag_sparse(array([1,5])).todense(),csr_matrix(A).todense())
+
+ A = matrix([[0,0],[0,-5]])
+ assert_equal(diag_sparse(array([0,-5])).todense(),csr_matrix(A).todense())
+
+ A = matrix([[1.3,0,0],[0,5.5,0],[0,0,-2]])
+ assert_equal(diag_sparse(array([1.3,5.5,-2])).todense(),csr_matrix(A).todense())
+
+
+
+
+if __name__ == '__main__':
+ NumpyTest().run()
+
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-09-07 19:55:57 UTC (rev 3307)
@@ -6,17 +6,18 @@
csr_matrix,csc_matrix,extract_diagonal
-def inf_norm(A):
+def infinity_norm(A):
"""
Infinity norm of a sparse matrix (maximum absolute row sum). This serves
as an upper bound on spectral radius.
"""
- if not isspmatrix_csr(A):
- return ValueError,'expected csr_matrix'
-
- abs_A = csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape,check=False)
- return (abs_A * numpy.ones(A.shape[1],dtype=A.dtype)).max()
+ if isspmatrix_csr(A) or isspmatrix_csc(A):
+ #avoid copying index and ptr arrays
+ abs_A = A.__class__((abs(A.data),A.indices,A.indptr),dims=A.shape,check=False)
+ return (abs_A * numpy.ones(A.shape[1],dtype=A.dtype)).max()
+ else:
+ return (abs(A) * numpy.ones(A.shape[1],dtype=A.dtype)).max()
def diag_sparse(A):
"""
More information about the Scipy-svn
mailing list