[Scipy-svn] r4414 - trunk/scipy/sparse/linalg/eigen/lobpcg
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Jun 6 02:29:52 EDT 2008
Author: wnbell
Date: 2008-06-06 01:29:48 -0500 (Fri, 06 Jun 2008)
New Revision: 4414
Modified:
trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
Log:
minor cleanup of lobpcg
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-06-04 16:21:10 UTC (rev 4413)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-06-06 06:29:48 UTC (rev 4414)
@@ -10,13 +10,11 @@
Examples in tests directory contributed by Nils Wagner.
"""
-import types
from warnings import warn
-import numpy as nm
-import scipy as sc
-import scipy.sparse as sp
-import scipy.io as io
+import numpy as np
+import scipy as sp
+
from scipy.sparse.linalg import aslinearoperator, LinearOperator
## try:
@@ -28,7 +26,7 @@
import scipy.linalg as sla
import scipy.lib.lapack as ll
if select is None:
- if nm.iscomplexobj( mtxA ):
+ if np.iscomplexobj( mtxA ):
if mtxB is None:
fun = ll.get_lapack_funcs( ['heev'], arrays = (mtxA,) )[0]
else:
@@ -51,7 +49,7 @@
else:
out = sla.eig( mtxA, mtxB, right = eigenvectors )
w = out[0]
- ii = nm.argsort( w )
+ ii = np.argsort( w )
w = w[slice( *select )]
if eigenvectors:
v = out[1][:,ii]
@@ -66,7 +64,8 @@
raw_input()
def save( ar, fileName ):
- io.write_array( fileName, ar, precision = 8 )
+ from scipy.io import write_array
+ write_array( fileName, ar, precision = 8 )
##
# 21.05.2007, c
@@ -78,7 +77,7 @@
if ar.ndim == 2:
return ar
else: # Assume 1!
- aux = nm.array( ar, copy = False )
+ aux = np.array( ar, copy = False )
aux.shape = (ar.shape[0], 1)
return aux
@@ -111,10 +110,10 @@
def applyConstraints( blockVectorV, factYBY, blockVectorBY, blockVectorY ):
"""Internal. Changes blockVectorV in place."""
- gramYBV = sc.dot( blockVectorBY.T, blockVectorV )
+ gramYBV = sp.dot( blockVectorBY.T, blockVectorV )
import scipy.linalg as sla
tmp = sla.cho_solve( factYBY, gramYBV )
- blockVectorV -= sc.dot( blockVectorY, tmp )
+ blockVectorV -= sp.dot( blockVectorY, tmp )
def b_orthonormalize( B, blockVectorV,
@@ -126,13 +125,13 @@
blockVectorBV = B( blockVectorV )
else:
blockVectorBV = blockVectorV # Shared data!!!
- gramVBV = sc.dot( blockVectorV.T, blockVectorBV )
+ gramVBV = sp.dot( blockVectorV.T, blockVectorBV )
gramVBV = sla.cholesky( gramVBV )
sla.inv( gramVBV, overwrite_a = True )
# gramVBV is now R^{-1}.
- blockVectorV = sc.dot( blockVectorV, gramVBV )
+ blockVectorV = sp.dot( blockVectorV, gramVBV )
if B is not None:
- blockVectorBV = sc.dot( blockVectorBV, gramVBV )
+ blockVectorBV = sp.dot( blockVectorBV, gramVBV )
if retInvR:
return blockVectorV, blockVectorBV, gramVBV
@@ -236,10 +235,10 @@
else:
lohi = (1, sizeX)
- A_dense = A(nm.eye(n))
+ A_dense = A(np.eye(n))
if B is not None:
- B_dense = B(nm.eye(n))
+ B_dense = B(np.eye(n))
_lambda, eigBlockVector = symeig(A_dense, B_dense, select=lohi )
else:
_lambda, eigBlockVector = symeig(A_dense, select=lohi )
@@ -248,7 +247,7 @@
if residualTolerance is None:
- residualTolerance = nm.sqrt( 1e-15 ) * n
+ residualTolerance = np.sqrt( 1e-15 ) * n
maxIterations = min( n, maxIterations )
@@ -283,7 +282,7 @@
blockVectorBY = blockVectorY
# gramYBY is a dense array.
- gramYBY = sc.dot( blockVectorY.T, blockVectorBY )
+ gramYBY = sp.dot( blockVectorY.T, blockVectorBY )
try:
# gramYBY is a Cholesky factor from now on...
gramYBY = sla.cho_factor( gramYBY )
@@ -299,32 +298,32 @@
##
# Compute the initial Ritz vectors: solve the eigenproblem.
blockVectorAX = A( blockVectorX )
- gramXAX = sc.dot( blockVectorX.T, blockVectorAX )
+ gramXAX = sp.dot( blockVectorX.T, blockVectorAX )
# gramXBX is X^T * X.
- gramXBX = sc.dot( blockVectorX.T, blockVectorX )
+ gramXBX = sp.dot( blockVectorX.T, blockVectorX )
_lambda, eigBlockVector = symeig( gramXAX )
- ii = nm.argsort( _lambda )[:sizeX]
+ ii = np.argsort( _lambda )[:sizeX]
if largest:
ii = ii[::-1]
_lambda = _lambda[ii]
- eigBlockVector = nm.asarray( eigBlockVector[:,ii] )
- blockVectorX = sc.dot( blockVectorX, eigBlockVector )
- blockVectorAX = sc.dot( blockVectorAX, eigBlockVector )
+ eigBlockVector = np.asarray( eigBlockVector[:,ii] )
+ blockVectorX = sp.dot( blockVectorX, eigBlockVector )
+ blockVectorAX = sp.dot( blockVectorAX, eigBlockVector )
if B is not None:
- blockVectorBX = sc.dot( blockVectorBX, eigBlockVector )
+ blockVectorBX = sp.dot( blockVectorBX, eigBlockVector )
##
# Active index set.
- activeMask = nm.ones( (sizeX,), dtype = nm.bool )
+ activeMask = np.ones( (sizeX,), dtype = np.bool )
lambdaHistory = [_lambda]
residualNormsHistory = []
previousBlockSize = sizeX
- ident = nm.eye( sizeX, dtype = A.dtype )
- ident0 = nm.eye( sizeX, dtype = A.dtype )
+ ident = np.eye( sizeX, dtype = A.dtype )
+ ident0 = np.eye( sizeX, dtype = A.dtype )
##
# Main iteration loop.
@@ -332,15 +331,15 @@
if verbosityLevel > 0:
print 'iteration %d' % iterationNumber
- aux = blockVectorBX * _lambda[nm.newaxis,:]
+ aux = blockVectorBX * _lambda[np.newaxis,:]
blockVectorR = blockVectorAX - aux
- aux = nm.sum( blockVectorR.conjugate() * blockVectorR, 0 )
- residualNorms = nm.sqrt( aux )
+ aux = np.sum( blockVectorR.conjugate() * blockVectorR, 0 )
+ residualNorms = np.sqrt( aux )
residualNormsHistory.append( residualNorms )
- ii = nm.where( residualNorms > residualTolerance, True, False )
+ ii = np.where( residualNorms > residualTolerance, True, False )
activeMask = activeMask & ii
if verbosityLevel > 2:
print activeMask
@@ -348,7 +347,7 @@
currentBlockSize = activeMask.sum()
if currentBlockSize != previousBlockSize:
previousBlockSize = currentBlockSize
- ident = nm.eye( currentBlockSize, dtype = A.dtype )
+ ident = np.eye( currentBlockSize, dtype = A.dtype )
if currentBlockSize == 0:
failureFlag = False # All eigenpairs converged.
@@ -390,44 +389,44 @@
aux = b_orthonormalize( B, activeBlockVectorP,
activeBlockVectorBP, retInvR = True )
activeBlockVectorP, activeBlockVectorBP, invR = aux
- activeBlockVectorAP = sc.dot( activeBlockVectorAP, invR )
+ activeBlockVectorAP = sp.dot( activeBlockVectorAP, invR )
##
# Perform the Rayleigh Ritz Procedure:
# Compute symmetric Gram matrices:
- xaw = sc.dot( blockVectorX.T, activeBlockVectorAR )
- waw = sc.dot( activeBlockVectorR.T, activeBlockVectorAR )
- xbw = sc.dot( blockVectorX.T, activeBlockVectorBR )
+ xaw = sp.dot( blockVectorX.T, activeBlockVectorAR )
+ waw = sp.dot( activeBlockVectorR.T, activeBlockVectorAR )
+ xbw = sp.dot( blockVectorX.T, activeBlockVectorBR )
if iterationNumber > 0:
- xap = sc.dot( blockVectorX.T, activeBlockVectorAP )
- wap = sc.dot( activeBlockVectorR.T, activeBlockVectorAP )
- pap = sc.dot( activeBlockVectorP.T, activeBlockVectorAP )
- xbp = sc.dot( blockVectorX.T, activeBlockVectorBP )
- wbp = sc.dot( activeBlockVectorR.T, activeBlockVectorBP )
+ xap = sp.dot( blockVectorX.T, activeBlockVectorAP )
+ wap = sp.dot( activeBlockVectorR.T, activeBlockVectorAP )
+ pap = sp.dot( activeBlockVectorP.T, activeBlockVectorAP )
+ xbp = sp.dot( blockVectorX.T, activeBlockVectorBP )
+ wbp = sp.dot( activeBlockVectorR.T, activeBlockVectorBP )
- gramA = nm.bmat( [[nm.diag( _lambda ), xaw, xap],
+ gramA = np.bmat( [[np.diag( _lambda ), xaw, xap],
[ xaw.T, waw, wap],
[ xap.T, wap.T, pap]] )
- gramB = nm.bmat( [[ident0, xbw, xbp],
+ gramB = np.bmat( [[ident0, xbw, xbp],
[ xbw.T, ident, wbp],
[ xbp.T, wbp.T, ident]] )
else:
- gramA = nm.bmat( [[nm.diag( _lambda ), xaw],
+ gramA = np.bmat( [[np.diag( _lambda ), xaw],
[ xaw.T, waw]] )
- gramB = nm.bmat( [[ident0, xbw],
+ gramB = np.bmat( [[ident0, xbw],
[ xbw.T, ident]] )
try:
- assert nm.allclose( gramA.T, gramA )
+ assert np.allclose( gramA.T, gramA )
except:
print gramA.T - gramA
raise
try:
- assert nm.allclose( gramB.T, gramB )
+ assert np.allclose( gramB.T, gramB )
except:
print gramB.T - gramB
raise
@@ -440,23 +439,23 @@
# Solve the generalized eigenvalue problem.
# _lambda, eigBlockVector = la.eig( gramA, gramB )
_lambda, eigBlockVector = symeig( gramA, gramB )
- ii = nm.argsort( _lambda )[:sizeX]
+ ii = np.argsort( _lambda )[:sizeX]
if largest:
ii = ii[::-1]
if verbosityLevel > 10:
print ii
- _lambda = _lambda[ii].astype( nm.float64 )
- eigBlockVector = nm.asarray( eigBlockVector[:,ii].astype( nm.float64 ) )
+ _lambda = _lambda[ii].astype( np.float64 )
+ eigBlockVector = np.asarray( eigBlockVector[:,ii].astype( np.float64 ) )
lambdaHistory.append( _lambda )
if verbosityLevel > 10:
print 'lambda:', _lambda
## # Normalize eigenvectors!
-## aux = nm.sum( eigBlockVector.conjugate() * eigBlockVector, 0 )
-## eigVecNorms = nm.sqrt( aux )
-## eigBlockVector = eigBlockVector / eigVecNorms[nm.newaxis,:]
+## aux = np.sum( eigBlockVector.conjugate() * eigBlockVector, 0 )
+## eigVecNorms = np.sqrt( aux )
+## eigBlockVector = eigBlockVector / eigVecNorms[np.newaxis,:]
# eigBlockVector, aux = b_orthonormalize( B, eigBlockVector )
if verbosityLevel > 10:
@@ -470,21 +469,21 @@
eigBlockVectorR = eigBlockVector[sizeX:sizeX+currentBlockSize]
eigBlockVectorP = eigBlockVector[sizeX+currentBlockSize:]
- pp = sc.dot( activeBlockVectorR, eigBlockVectorR )
- pp += sc.dot( activeBlockVectorP, eigBlockVectorP )
+ pp = sp.dot( activeBlockVectorR, eigBlockVectorR )
+ pp += sp.dot( activeBlockVectorP, eigBlockVectorP )
- app = sc.dot( activeBlockVectorAR, eigBlockVectorR )
- app += sc.dot( activeBlockVectorAP, eigBlockVectorP )
+ app = sp.dot( activeBlockVectorAR, eigBlockVectorR )
+ app += sp.dot( activeBlockVectorAP, eigBlockVectorP )
- bpp = sc.dot( activeBlockVectorBR, eigBlockVectorR )
- bpp += sc.dot( activeBlockVectorBP, eigBlockVectorP )
+ bpp = sp.dot( activeBlockVectorBR, eigBlockVectorR )
+ bpp += sp.dot( activeBlockVectorBP, eigBlockVectorP )
else:
eigBlockVectorX = eigBlockVector[:sizeX]
eigBlockVectorR = eigBlockVector[sizeX:]
- pp = sc.dot( activeBlockVectorR, eigBlockVectorR )
- app = sc.dot( activeBlockVectorAR, eigBlockVectorR )
- bpp = sc.dot( activeBlockVectorBR, eigBlockVectorR )
+ pp = sp.dot( activeBlockVectorR, eigBlockVectorR )
+ app = sp.dot( activeBlockVectorAR, eigBlockVectorR )
+ bpp = sp.dot( activeBlockVectorBR, eigBlockVectorR )
if verbosityLevel > 10:
print pp
@@ -492,17 +491,17 @@
print bpp
pause()
- blockVectorX = sc.dot( blockVectorX, eigBlockVectorX ) + pp
- blockVectorAX = sc.dot( blockVectorAX, eigBlockVectorX ) + app
- blockVectorBX = sc.dot( blockVectorBX, eigBlockVectorX ) + bpp
+ blockVectorX = sp.dot( blockVectorX, eigBlockVectorX ) + pp
+ blockVectorAX = sp.dot( blockVectorAX, eigBlockVectorX ) + app
+ blockVectorBX = sp.dot( blockVectorBX, eigBlockVectorX ) + bpp
blockVectorP, blockVectorAP, blockVectorBP = pp, app, bpp
- aux = blockVectorBX * _lambda[nm.newaxis,:]
+ aux = blockVectorBX * _lambda[np.newaxis,:]
blockVectorR = blockVectorAX - aux
- aux = nm.sum( blockVectorR.conjugate() * blockVectorR, 0 )
- residualNorms = nm.sqrt( aux )
+ aux = np.sum( blockVectorR.conjugate() * blockVectorR, 0 )
+ residualNorms = np.sqrt( aux )
if verbosityLevel > 0:
@@ -522,31 +521,31 @@
###########################################################################
if __name__ == '__main__':
- from scipy.sparse import spdiags, speye
+ from scipy.sparse import spdiags, speye, issparse
import time
## def B( vec ):
## return vec
n = 100
- vals = [nm.arange( n, dtype = nm.float64 ) + 1]
+ vals = [np.arange( n, dtype = np.float64 ) + 1]
A = spdiags( vals, 0, n, n )
B = speye( n, n )
# B[0,0] = 0
- B = nm.eye( n, n )
- Y = nm.eye( n, 3 )
+ B = np.eye( n, n )
+ Y = np.eye( n, 3 )
-# X = sc.rand( n, 3 )
+# X = sp.rand( n, 3 )
xfile = {100 : 'X.txt', 1000 : 'X2.txt', 10000 : 'X3.txt'}
- X = nm.fromfile( xfile[n], dtype = nm.float64, sep = ' ' )
+ X = np.fromfile( xfile[n], dtype = np.float64, sep = ' ' )
X.shape = (n, 3)
ivals = [1./vals[0]]
def precond( x ):
invA = spdiags( ivals, 0, n, n )
y = invA * x
- if sp.issparse( y ):
+ if issparse( y ):
y = y.toarray()
return as2d( y )
More information about the Scipy-svn
mailing list