[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