[Scipy-svn] r4059 - in trunk/scipy/sparse/linalg/eigen: . lobpcg lobpcg/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Mar 31 16:51:46 EDT 2008
Author: wnbell
Date: 2008-03-31 15:51:41 -0500 (Mon, 31 Mar 2008)
New Revision: 4059
Modified:
trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
trunk/scipy/sparse/linalg/eigen/lobpcg/tests/test_lobpcg.py
trunk/scipy/sparse/linalg/eigen/setup.py
Log:
reworked LOBPCG unittests
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-03-31 17:35:30 UTC (rev 4058)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py 2008-03-31 20:51:41 UTC (rev 4059)
@@ -296,8 +296,7 @@
# gramYBY is a Cholesky factor from now on...
gramYBY = la.cho_factor( gramYBY )
except:
- print 'cannot handle linear dependent constraints'
- raise
+ raise ValueError('cannot handle linear dependent constraints')
applyConstraints( blockVectorX, gramYBY, blockVectorBY, blockVectorY )
Modified: trunk/scipy/sparse/linalg/eigen/lobpcg/tests/test_lobpcg.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/lobpcg/tests/test_lobpcg.py 2008-03-31 17:35:30 UTC (rev 4058)
+++ trunk/scipy/sparse/linalg/eigen/lobpcg/tests/test_lobpcg.py 2008-03-31 20:51:41 UTC (rev 4059)
@@ -1,11 +1,25 @@
+#!/usr/bin/env python
+""" Test functions for the sparse.linalg.eigen.lobpcg module
+"""
+
+import numpy
+from scipy.testing import *
+
from scipy import array, arange, ones, sort, cos, pi, rand, \
set_printoptions, r_, diag, linalg
-from scipy.sandbox import lobpcg
-from symeig import symeig
-from pylab import plot, show, legend, xlabel, ylabel
+from scipy.sparse.linalg.eigen import lobpcg
+
+try:
+ from symeig import symeig
+except:
+ raise ImportError('lobpcg requires symeig')
+
set_printoptions(precision=3,linewidth=90)
-def check1(n):
+
+
+def ElasticRod(n):
+ # Fixed-free elastic rod
L = 1.0
le=L/n
rho = 7.85e3
@@ -17,7 +31,9 @@
B = mass*(diag(r_[4.*ones(n-1),2])+diag(ones(n-1),1)+diag(ones(n-1),-1))
return A,B
-def check2(n):
+def MikotaPair(n):
+ # Mikota pair acts as a nice test since the eigenvalues
+ # are the squares of the integers n, n=1,2,...
x = arange(1,n+1)
B = diag(1./x)
y = arange(n-1,0,-1)
@@ -25,24 +41,39 @@
A = diag(z)-diag(y,-1)-diag(y,1)
return A,B
-n = 100 # Dimension
-def test_1and2():
- A,B = check1(n) # Fixed-free elastic rod
- A,B = check2(n) # Mikota pair acts as a nice test since the eigenvalues are the squares of the integers n, n=1,2,...
-
- m = 20
+def compare_solutions(A,B,m):
+ n = A.shape[0]
+
+ numpy.random.seed(0)
+
V = rand(n,m)
X = linalg.orth(V)
- eigs,vecs = lobpcg.lobpcg(X,A,B)
- eigs = sort(eigs)
+ eigs,vecs = lobpcg.lobpcg(X,A,B,residualTolerance=1e-5, maxIterations=30)
+ eigs.sort()
- w,v=symeig(A,B)
+ w,v = symeig(A,B)
- plot(arange(0,len(w[:m])),w[:m],'bx',label='Results by symeig')
- plot(arange(0,len(eigs)),eigs,'r+',label='Results by lobpcg')
- legend()
- xlabel(r'Eigenvalue $i$')
- ylabel(r'$\lambda_i$')
- show()
+ assert_almost_equal(w[:m/2],eigs[:m/2],decimal=2)
+
+ #from pylab import plot, show, legend, xlabel, ylabel
+ #plot(arange(0,len(w[:m])),w[:m],'bx',label='Results by symeig')
+ #plot(arange(0,len(eigs)),eigs,'r+',label='Results by lobpcg')
+ #legend()
+ #xlabel(r'Eigenvalue $i$')
+ #ylabel(r'$\lambda_i$')
+ #show()
+
+
+def test_ElasticRod():
+ A,B = ElasticRod(100)
+ compare_solutions(A,B,20)
+
+def test_MikotaPair():
+ A,B = MikotaPair(100)
+ compare_solutions(A,B,20)
+
+
+if __name__ == "__main__":
+ nose.run(argv=['', __file__])
Modified: trunk/scipy/sparse/linalg/eigen/setup.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/setup.py 2008-03-31 17:35:30 UTC (rev 4058)
+++ trunk/scipy/sparse/linalg/eigen/setup.py 2008-03-31 20:51:41 UTC (rev 4059)
@@ -6,6 +6,7 @@
config = Configuration('eigen',parent_package,top_path)
config.add_subpackage(('arpack'))
+ #config.add_subpackage(('lobpcg'))
return config
More information about the Scipy-svn
mailing list