#!/usr/bin/env python
""" Test functions for the sparse.linalg.eigen.lobpcg module
"""
import numpy
from numpy.testing import *
from scipy import arange,ones,rand,set_printoptions,r_,diag,linalg
from scipy.linalg import eig
from scipy.sparse.linalg.eigen.lobpcg import lobpcg
set_printoptions(precision=3,linewidth=90)
def ElasticRod(n):
# Fixed-free elastic rod
L = 1.0
le=L/n
rho = 7.85e3
S = 1.e-4
E = 2.1e11
mass = rho*S*le/6.
k = E*S/le
A = k*(diag(r_[2.*ones(n-1),1])-diag(ones(n-1),1)-diag(ones(n-1),-1))
B = mass*(diag(r_[4.*ones(n-1),2])+diag(ones(n-1),1)+diag(ones(n-1),-1))
return A,B
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)
z = arange(2*n-1,0,-2)
A = diag(z)-diag(y,-1)-diag(y,1)
return A,B
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(A, X, B=B, tol=1e-5, maxiter=30)
eigs.sort()
#w,v = symeig(A,B)
w,v = eig(A,b=B)
w.sort()
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_Small():
A,B = ElasticRod(10)
compare_solutions(A,B,10)
A,B = MikotaPair(10)
compare_solutions(A,B,10)
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__":
run_module_suite()
|