[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