[Scipy-svn] r3440 - in trunk/scipy/sandbox/multigrid: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Oct 16 02:19:15 EDT 2007
Author: wnbell
Date: 2007-10-16 01:19:12 -0500 (Tue, 16 Oct 2007)
New Revision: 3440
Modified:
trunk/scipy/sandbox/multigrid/adaptive.py
trunk/scipy/sandbox/multigrid/sa.py
trunk/scipy/sandbox/multigrid/tests/test_sa.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
add block option to adaptive SA
Modified: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/adaptive.py 2007-10-16 06:19:12 UTC (rev 3440)
@@ -105,13 +105,13 @@
self.Rs = []
- if self.A.shape[0] <= self.opts['coarse: max size']:
+ if self.A.shape[0] <= max_coarse:
raise ValueError,'small matrices not handled yet'
#first candidate
- x,AggOps = self.__initialization_stage(A,blocks=blocks,\
- max_levels=max_levels,max_coarse=max_coarse,\
- mu=mu,epsilon=epsilon)
+ x,AggOps = self.__initialization_stage(A, blocks=blocks,\
+ max_levels=max_levels, max_coarse=max_coarse,\
+ mu=mu, epsilon=epsilon)
Ws = AggOps
@@ -134,7 +134,7 @@
self.AggOps = AggOps
- def __initialization_stage(self,A,max_levels,max_coarse,mu,epsilon):
+ def __initialization_stage(self,A,blocks,max_levels,max_coarse,mu,epsilon):
AggOps = []
Is = []
@@ -270,6 +270,9 @@
from scipy import *
from utils import diag_sparse
from multilevel import poisson_problem1D,poisson_problem2D
+
+blocks = None
+
A = poisson_problem2D(200)
#A = io.mmread("tests/sample_data/laplacian_41_3dcube.mtx").tocsr()
#A = io.mmread("laplacian_40_3dcube.mtx").tocsr()
@@ -280,16 +283,17 @@
#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()
+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=5)
+asa = adaptive_sa_solver(A,max_candidates=4,mu=12)
scipy.random.seed(0) #make tests repeatable
x = rand(A.shape[0])
b = A*rand(A.shape[0])
print "solving"
-if False:
+if True:
x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
else:
residuals = []
Modified: trunk/scipy/sandbox/multigrid/sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/sa.py 2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/sa.py 2007-10-16 06:19:12 UTC (rev 3440)
@@ -21,34 +21,37 @@
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:
- num_dofs = A.shape[0]
- num_blocks = blocks.max() + 1
-
- if num_dofs != len(blocks):
- raise ValueError,'improper block specification'
-
- # for non-scalar problems, use pre-defined blocks in aggregation
- # the strength of connection matrix is based on the Frobenius norms of the blocks
-
- B = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
- Bt = B.T.tocsr()
-
- #1-norms of blocks entries of A
- Block_A = Bt * csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape) * B
-
- S = sa_strong_connections(Block_A,epsilon)
- S.data[:] = 1
+ A_filtered = A #TODO subtract weak blocks from diagonal blocks?
- Mask = B * S * Bt
+## num_dofs = A.shape[0]
+## num_blocks = blocks.max() + 1
+##
+## if num_dofs != len(blocks):
+## raise ValueError,'improper block specification'
+##
+## # for non-scalar problems, use pre-defined blocks in aggregation
+## # the strength of connection matrix is based on the 1-norms of the blocks
+##
+## B = csr_matrix((ones(num_dofs),blocks,arange(num_dofs + 1)),dims=(num_dofs,num_blocks))
+## Bt = B.T.tocsr()
+##
+## #1-norms of blocks entries of A
+## Block_A = Bt * csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape) * B
+##
+## S = sa_strong_connections(Block_A,epsilon)
+## S.data[:] = 1
+##
+## Mask = B * S * Bt
+##
+## A_strong = A ** Mask
+## #A_weak = A - A_strong
+## A_filtered = A_strong
- A_filtered = A ** Mask
-
return A_filtered
-def sa_strong_connections(A,epsilon,blocks=None):
+def sa_strong_connections(A,epsilon):
if not isspmatrix_csr(A): raise TypeError('expected csr_matrix')
Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
@@ -66,7 +69,7 @@
if num_dofs != len(blocks):
raise ValueError,'improper block specification'
- print "SA has blocks"
+ #print "SA has blocks"
# for non-scalar problems, use pre-defined blocks in aggregation
# the strength of connection matrix is based on the Frobenius norms of the blocks
Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py 2007-10-16 06:19:12 UTC (rev 3440)
@@ -1,7 +1,7 @@
from numpy.testing import *
from numpy import sqrt,empty,ones,arange,array_split,eye,array, \
- zeros,diag,zeros_like,diff
+ zeros,diag,zeros_like,diff,matrix
from numpy.linalg import norm
from scipy import rand
from scipy.sparse import spdiags,csr_matrix,lil_matrix, \
@@ -51,19 +51,42 @@
def check_sa_strong_connections(self):
for A in self.cases:
for epsilon in [0.0,0.1,0.5,1.0,10.0]:
+ S_expected = reference_sa_strong_connections(A,epsilon)
S_result = sa_strong_connections(A,epsilon)
- S_expected = reference_sa_strong_connections(A,epsilon)
assert_almost_equal(S_result.todense(),S_expected.todense())
#assert_array_equal(sparsity(S_result).todense(),sparsity(S_expected).todense())
def check_sa_constant_interpolation(self):
for A in self.cases:
for epsilon in [0.0,0.1,0.5,1.0]:
- S_result = sa_constant_interpolation(A,epsilon)
S_expected = reference_sa_constant_interpolation(A,epsilon)
+
+ S_result = sa_constant_interpolation(A,epsilon,blocks=None)
assert_array_equal(S_result.todense(),S_expected.todense())
+
+ #blocks=1...N should be the same as blocks=None
+ S_result = sa_constant_interpolation(A,epsilon,blocks=arange(A.shape[0]))
+ assert_array_equal(S_result.todense(),S_expected.todense())
+ #check simple block examples
+ A = csr_matrix(arange(16).reshape(4,4))
+ A = A + A.T
+ blocks = array([0,0,1,1])
+ S_result = sa_constant_interpolation(A,epsilon=0.0,blocks=blocks)
+ S_expected = matrix([1,1,1,1]).T
+ assert_array_equal(S_result.todense(),S_expected)
+
+ S_result = sa_constant_interpolation(A,epsilon=0.5,blocks=blocks)
+ S_expected = matrix([1,1,1,1]).T
+ assert_array_equal(S_result.todense(),S_expected)
+
+ S_result = sa_constant_interpolation(A,epsilon=2.0,blocks=blocks)
+ S_expected = matrix([[1,0],[1,0],[0,1],[0,1]])
+ assert_array_equal(S_result.todense(),S_expected)
+
+
+
class TestFitCandidates(NumpyTestCase):
def setUp(self):
self.cases = []
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-10-15 22:46:49 UTC (rev 3439)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-10-16 06:19:12 UTC (rev 3440)
@@ -3,10 +3,11 @@
import numpy
import scipy
-from numpy import ravel,arange
+from numpy import ravel,arange,concatenate
from scipy.linalg import norm
from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
- csr_matrix,csc_matrix,extract_diagonal
+ csr_matrix,csc_matrix,extract_diagonal, \
+ coo_matrix
def approximate_spectral_radius(A,tol=0.1,maxiter=20):
More information about the Scipy-svn
mailing list