[Scipy-svn] r3617 - trunk/scipy/sandbox/multigrid
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Dec 5 14:12:46 EST 2007
Author: wnbell
Date: 2007-12-05 13:12:44 -0600 (Wed, 05 Dec 2007)
New Revision: 3617
Modified:
trunk/scipy/sandbox/multigrid/utils.py
Log:
removed use of ARPACK
changed diagonal scaling
Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py 2007-12-05 07:23:34 UTC (rev 3616)
+++ trunk/scipy/sandbox/multigrid/utils.py 2007-12-05 19:12:44 UTC (rev 3617)
@@ -3,21 +3,50 @@
import numpy
import scipy
-from numpy import ravel,arange,concatenate,tile,asarray,sqrt,diff
-from scipy.linalg import norm
+from scipy import ravel,arange,concatenate,tile,asarray,sqrt,diff, \
+ rand,zeros,empty,asmatrix,dot
+from scipy.linalg import norm,eig
from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
csr_matrix,csc_matrix,extract_diagonal, \
coo_matrix
-def approximate_spectral_radius(A,tol=0.1,maxiter=None):
+def approximate_spectral_radius(A,tol=0.1,maxiter=8):
"""
- Approximate the spectral radius of a symmetric matrix using ARPACK
+ Approximate the spectral radius of a matrix
"""
- from scipy.sandbox.arpack import eigen
- return norm(eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False))
+ #from scipy.sandbox.arpack import eigen
+ #return norm(eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False))
+
+ numpy.random.seed(0) #make results deterministic
+ #TODO profile vs V -> V.T
+ #TODO make algorithm adaptive
+ if not isspmatrix(A):
+ #convert dense arrays to matrix type
+ A = asmatrix(A)
+
+ v0 = rand(A.shape[0])
+ H = zeros((maxiter+1,maxiter))
+ V = zeros((A.shape[0],maxiter+1))
+
+ V[:,0]= v0/norm(v0)
+ for j in range(maxiter):
+ w = A * V[:,j]
+ for i in range(j+1):
+ H[i,j] = dot(w,V[:,i])
+ w -= H[i,j]*V[:,i]
+ H[j+1,j] = norm(w)
+ if (H[j+1,j] < 1e-12): break
+ V[:,j+1] = (1.0/H[j+1,j]) * w
+ # end
+ m=j
+ Heig,tmp = eig(H[0:m,0:m])
+ return max( [norm(x) for x in Heig] )
+
+
+
def infinity_norm(A):
"""
Infinity norm of a sparse matrix (maximum absolute row sum). This serves
@@ -55,10 +84,10 @@
raise ValueError,'expected square matrix'
D = diag_sparse(A)
- mask = D <= 0
+ mask = D == 0
- D[mask] = 0
- D_sqrt = sqrt(D)
+ #D[mask] = 0
+ D_sqrt = sqrt(abs(D))
D_sqrt_inv = 1.0/D_sqrt
D_sqrt_inv[mask] = 0
More information about the Scipy-svn
mailing list