[Scipy-svn] r4063 - in trunk/scipy/sparse/linalg/eigen/lobpcg: . tests
scipy-svn@scip...
scipy-svn@scip...
Tue Apr 1 11:10:34 CDT 2008
Author: wnbell
Date: 2008-04-01 11:10:24 -0500 (Tue, 01 Apr 2008)
New Revision: 4063
Modified:
trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
trunk/scipy/sparse/linalg/eigen/lobpcg/tests/large_scale.py
Log:
minor cleanup
updated test
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-04-01 15:46:48 UTC (rev 4062)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-04-01 16:10:24 UTC (rev 4063)
@@ -251,7 +251,7 @@
# gramYBY is a Cholesky factor from now on...
gramYBY = la.cho_factor( gramYBY )
except:
- raise ValueError('cannot handle linear dependent constraints')
+ raise ValueError('cannot handle linearly dependent constraints')
applyConstraints( blockVectorX, gramYBY, blockVectorBY, blockVectorY )
@@ -265,14 +265,15 @@
gramXAX = sc.dot( blockVectorX.T, blockVectorAX )
# gramXBX is X^T * X.
gramXBX = sc.dot( blockVectorX.T, blockVectorX )
+
_lambda, eigBlockVector = symeig( gramXAX )
ii = nm.argsort( _lambda )[:sizeX]
if largest:
ii = ii[::-1]
_lambda = _lambda[ii]
+
eigBlockVector = nm.asarray( eigBlockVector[:,ii] )
-# pause()
- blockVectorX = sc.dot( blockVectorX, eigBlockVector )
+ blockVectorX = sc.dot( blockVectorX, eigBlockVector )
blockVectorAX = sc.dot( blockVectorAX, eigBlockVector )
if B is not None:
blockVectorBX = sc.dot( blockVectorBX, eigBlockVector )
@@ -285,7 +286,7 @@
residualNormsHistory = []
previousBlockSize = sizeX
- ident = nm.eye( sizeX, dtype = A.dtype )
+ ident = nm.eye( sizeX, dtype = A.dtype )
ident0 = nm.eye( sizeX, dtype = A.dtype )
##
@@ -369,22 +370,19 @@
xbp = sc.dot( blockVectorX.T, activeBlockVectorBP )
wbp = sc.dot( activeBlockVectorR.T, activeBlockVectorBP )
- gramA = nm.bmat( [[nm.diag( _lambda ), xaw, xap],
- [xaw.T, waw, wap],
- [xap.T, wap.T, pap]] )
- try:
- gramB = nm.bmat( [[ident0, xbw, xbp],
- [ xbw.T, ident, wbp],
- [ xbp.T, wbp.T, ident]] )
- except:
- print ident
- print xbw
- raise
+ gramA = nm.bmat( [[nm.diag( _lambda ), xaw, xap],
+ [ xaw.T, waw, wap],
+ [ xap.T, wap.T, pap]] )
+
+ gramB = nm.bmat( [[ident0, xbw, xbp],
+ [ xbw.T, ident, wbp],
+ [ xbp.T, wbp.T, ident]] )
else:
- gramA = nm.bmat( [[nm.diag( _lambda ), xaw],
- [xaw.T, waw]] )
- gramB = nm.bmat( [[ident0, xbw],
- [xbw.T, ident0]] )
+ gramA = nm.bmat( [[nm.diag( _lambda ), xaw],
+ [ xaw.T, waw]] )
+ gramB = nm.bmat( [[ident0, xbw],
+ [ xbw.T, ident0]] )
+
try:
assert nm.allclose( gramA.T, gramA )
except:
@@ -427,6 +425,7 @@
if verbosityLevel > 10:
print eigBlockVector
pause()
+
##
# Compute Ritz vectors.
if iterationNumber > 0:
@@ -434,22 +433,20 @@
eigBlockVectorR = eigBlockVector[sizeX:sizeX+currentBlockSize]
eigBlockVectorP = eigBlockVector[sizeX+currentBlockSize:]
- pp = sc.dot( activeBlockVectorR, eigBlockVectorR )\
- + sc.dot( activeBlockVectorP, eigBlockVectorP )
+ pp = sc.dot( activeBlockVectorR, eigBlockVectorR )
+ pp += sc.dot( activeBlockVectorP, eigBlockVectorP )
- app = sc.dot( activeBlockVectorAR, eigBlockVectorR )\
- + sc.dot( activeBlockVectorAP, eigBlockVectorP )
+ app = sc.dot( activeBlockVectorAR, eigBlockVectorR )
+ app += sc.dot( activeBlockVectorAP, eigBlockVectorP )
- bpp = sc.dot( activeBlockVectorBR, eigBlockVectorR )\
- + sc.dot( activeBlockVectorBP, eigBlockVectorP )
+ bpp = sc.dot( activeBlockVectorBR, eigBlockVectorR )
+ bpp += sc.dot( activeBlockVectorBP, eigBlockVectorP )
else:
eigBlockVectorX = eigBlockVector[:sizeX]
eigBlockVectorR = eigBlockVector[sizeX:]
- pp = sc.dot( activeBlockVectorR, eigBlockVectorR )
-
+ pp = sc.dot( activeBlockVectorR, eigBlockVectorR )
app = sc.dot( activeBlockVectorAR, eigBlockVectorR )
-
bpp = sc.dot( activeBlockVectorBR, eigBlockVectorR )
if verbosityLevel > 10:
@@ -457,9 +454,8 @@
print app
print bpp
pause()
-# print pp.shape, app.shape, bpp.shape
- blockVectorX = sc.dot( blockVectorX, eigBlockVectorX ) + pp
+ blockVectorX = sc.dot( blockVectorX, eigBlockVectorX ) + pp
blockVectorAX = sc.dot( blockVectorAX, eigBlockVectorX ) + app
blockVectorBX = sc.dot( blockVectorBX, eigBlockVectorX ) + bpp
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/tests/large_scale.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/tests/large_scale.py 2008-04-01 15:46:48 UTC (rev 4062)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/tests/large_scale.py 2008-04-01 16:10:24 UTC (rev 4063)
@@ -1,7 +1,7 @@
from scipy import array, arange, ones, sort, cos, pi, rand, \
set_printoptions, r_
-from scipy.sandbox import lobpcg
-from scipy.sparse import spdiags, speye
+from scipy.sparse.linalg import lobpcg
+from scipy import sparse
from pylab import loglog, show, xlabel, ylabel, title
set_printoptions(precision=8,linewidth=90)
import time
@@ -12,11 +12,11 @@
A moment-based method for large-scale generalized eigenvalue problems
Appl. Num. Anal. Comp. Math. Vol. 1 No. 2 (2004) """
- A = speye( n, n )
+ A = sparse.eye( n, n )
d0 = array(r_[5,6*ones(n-2),5])
d1 = -4*ones(n)
d2 = ones(n)
- B = spdiags([d2,d1,d0,d1,d2],[-2,-1,0,1,2],n,n)
+ B = sparse.spdiags([d2,d1,d0,d1,d2],[-2,-1,0,1,2],n,n)
k = arange(1,n+1)
w_ex = sort(1./(16.*pow(cos(0.5*k*pi/(n+1)),4))) # exact eigenvalues
@@ -28,13 +28,12 @@
#
# Large scale
#
-n = 25000
+n = 2500
A,B, w_ex = sakurai(n) # Mikota pair
X = rand(n,m)
data=[]
tt = time.clock()
-eigs,vecs, resnh = lobpcg.lobpcg(X,A,B,
- residualTolerance = 1e-6, maxIterations =500, retResidualNormsHistory=1)
+eigs,vecs, resnh = lobpcg(X,A,B, residualTolerance = 1e-6, maxIterations =500, retResidualNormsHistory=1)
data.append(time.clock()-tt)
print 'Results by LOBPCG for n='+str(n)
print
More information about the Scipy-svn
mailing list