[Scipy-svn] r3843 - in trunk/scipy/sandbox/multigrid: . examples tests
scipy-svn@scip...
scipy-svn@scip...
Wed Jan 16 15:57:57 CST 2008
Author: wnbell
Date: 2008-01-16 15:57:52 -0600 (Wed, 16 Jan 2008)
New Revision: 3843
Modified:
trunk/scipy/sandbox/multigrid/examples/adaptive.py
trunk/scipy/sandbox/multigrid/tests/test_sa.py
trunk/scipy/sandbox/multigrid/tests/test_utils.py
trunk/scipy/sandbox/multigrid/utils.py
Log:
fixed error in approximate_spectral_radius
Modified: trunk/scipy/sandbox/multigrid/examples/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/examples/adaptive.py 2008-01-16 13:15:30 UTC (rev 3842)
+++ trunk/scipy/sandbox/multigrid/examples/adaptive.py 2008-01-16 21:57:52 UTC (rev 3843)
@@ -10,7 +10,7 @@
#A = poisson( (200,200), spacing=(1,1e-2) ) #anisotropic
D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
-A = D * A * D
+#A = D * A * D
Modified: trunk/scipy/sandbox/multigrid/tests/test_sa.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_sa.py 2008-01-16 13:15:30 UTC (rev 3842)
+++ trunk/scipy/sandbox/multigrid/tests/test_sa.py 2008-01-16 21:57:52 UTC (rev 3843)
@@ -202,9 +202,9 @@
x_sol,residuals = ml.solve(b,x0=x,maxiter=20,tol=1e-12,return_residuals=True)
avg_convergence_ratio = (residuals[-1]/residuals[0])**(1.0/len(residuals))
+
+ assert(avg_convergence_ratio < 0.25)
- assert(avg_convergence_ratio < 0.5)
-
def test_DAD(self):
A = poisson( (100,100), format='csr' )
@@ -218,20 +218,17 @@
B = ones((A.shape[0],1))
- Dinv_B = D_inv * B
-
#TODO force 2 level method and check that result is the same
- sa1 = smoothed_aggregation_solver(A, B, max_levels=2, rescale=False)
- sa2 = smoothed_aggregation_solver(D*A*D, D_inv * B, max_levels=2, rescale=False)
+ #sa1 = smoothed_aggregation_solver(A, B, max_levels=2, rescale=False)
+ sa2 = smoothed_aggregation_solver(D*A*D, D_inv * B, max_levels=2, rescale=True)
#assert_almost_equal( sa2.Ps[0], sa1.Ps[0]
x_sol,residuals = sa2.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
avg_convergence_ratio = (residuals[-1]/residuals[0])**(1.0/len(residuals))
- print avg_convergence_ratio
-
- assert(avg_convergence_ratio < 0.2)
+
+ assert(avg_convergence_ratio < 0.25)
Modified: trunk/scipy/sandbox/multigrid/tests/test_utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_utils.py 2008-01-16 13:15:30 UTC (rev 3842)
+++ trunk/scipy/sandbox/multigrid/tests/test_utils.py 2008-01-16 21:57:52 UTC (rev 3843)
@@ -1,39 +1,38 @@
from scipy.testing import *
-import numpy
-import scipy
-from numpy import matrix,array,diag,zeros,sqrt
+from numpy import matrix, array, diag, zeros, sqrt
from scipy import rand
from scipy.sparse import csr_matrix
-from scipy.linalg import norm
+from scipy.linalg import eigvals, norm
-
from scipy.sandbox.multigrid.utils import *
from scipy.sandbox.multigrid.utils import symmetric_rescaling
class TestUtils(TestCase):
-# def test_approximate_spectral_radius(self):
-# cases = []
-#
-# cases.append( matrix([[-4]]) )
-# cases.append( array([[-4]]) )
-#
-# cases.append( array([[2,0],[0,1]]) )
-# cases.append( array([[-2,0],[0,1]]) )
-#
-# cases.append( array([[100,0,0],[0,101,0],[0,0,99]]) )
-#
-# for i in range(1,5):
-# cases.append( rand(i,i) )
-#
-# # method should be almost exact for small matrices
-# for A in cases:
-# A = A.astype(float)
-# Asp = csr_matrix(A)
-# assert_almost_equal( approximate_spectral_radius(A,tol=1e-2), norm(A,2), decimal=1 )
-# assert_almost_equal( approximate_spectral_radius(Asp,tol=1e-2), norm(A,2), decimal=1 )
-#
-# #TODO test larger matrices
+ def test_approximate_spectral_radius(self):
+ cases = []
+
+ cases.append( matrix([[-4]]) )
+ cases.append( array([[-4]]) )
+
+ cases.append( array([[2,0],[0,1]]) )
+ cases.append( array([[-2,0],[0,1]]) )
+
+ cases.append( array([[100,0,0],[0,101,0],[0,0,99]]) )
+
+ for i in range(1,5):
+ cases.append( rand(i,i) )
+
+ # method should be almost exact for small matrices
+ for A in cases:
+ A = A.astype(float)
+ Asp = csr_matrix(A)
+
+ expected = max([norm(x) for x in eigvals(A)])
+ assert_almost_equal( approximate_spectral_radius(A), expected )
+ assert_almost_equal( approximate_spectral_radius(Asp), expected )
+
+ #TODO test larger matrices
def test_infinity_norm(self):
A = matrix([[-4]])
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2008-01-16 13:15:30 UTC (rev 3842)
+++ trunk/scipy/sandbox/multigrid/utils.py 2008-01-16 21:57:52 UTC (rev 3843)
@@ -3,16 +3,16 @@
import numpy
import scipy
-from scipy import ravel,arange,concatenate,tile,asarray,sqrt,diff, \
- rand,zeros,empty,asmatrix,dot
-from scipy.linalg import norm,eigvals
+from scipy import ravel, arange, concatenate, tile, asarray, sqrt, diff, \
+ rand, zeros, ones, empty, asmatrix, dot
+from scipy.linalg import norm, eigvals
from scipy.sparse import isspmatrix, isspmatrix_csr, isspmatrix_csc, \
isspmatrix_bsr, csr_matrix, csc_matrix, bsr_matrix, coo_matrix, \
extract_diagonal
from scipy.sparse.sputils import upcast
-def approximate_spectral_radius(A,tol=0.1,maxiter=20,symmetric=None):
+def approximate_spectral_radius(A,tol=0.1,maxiter=10,symmetric=None):
"""approximate the spectral radius of a matrix
*Parameters*:
@@ -36,8 +36,8 @@
An approximation to the spectral radius of A (scalar value)
"""
- from scipy.sandbox.arpack import eigen
- return norm(eigen(A, k=1, ncv=min(10,A.shape[0]), which='LM', tol=tol, return_eigenvectors=False))
+ #from scipy.sandbox.arpack import eigen
+ #return norm(eigen(A, k=1, ncv=min(10,A.shape[0]), which='LM', tol=tol, return_eigenvectors=False))
if not isspmatrix(A):
A = asmatrix(A) #convert dense arrays to matrix type
@@ -98,8 +98,8 @@
V = V[1:]
H[1,0] = H[0,1]
beta = H[2,1]
-
- return norm(H[:j+1,:j+1],2)
+
+ return max([norm(x) for x in eigvals(H[:j+1,:j+1])])
@@ -112,9 +112,9 @@
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),shape=A.shape)
- return (abs_A * numpy.ones(A.shape[1],dtype=A.dtype)).max()
+ return (abs_A * ones(A.shape[1],dtype=A.dtype)).max()
else:
- return (abs(A) * numpy.ones(A.shape[1],dtype=A.dtype)).max()
+ return (abs(A) * ones(A.shape[1],dtype=A.dtype)).max()
def diag_sparse(A):
"""
More information about the Scipy-svn
mailing list