lobpcg_sakurai.py

"""
LOBPCG: block-preconditionned solver
=======================================

This example demos the LOBPCG block-preconditionned solver.
"""

from scipy import array, arange, ones, sort, cos, pi, rand, \
     set_printoptions, r_
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

def sakurai(n):
    """ Example taken from
        T. Sakurai, H. Tadano, Y. Inadomi and U. Nagashima
        A moment-based method for large-scale generalized eigenvalue problems
        Appl. Num. Anal. Comp. Math. Vol. 1 No. 2 (2004) """

    A = sparse.eye( n, n )
    d0 = array(r_[5,6*ones(n-2),5])
    d1 = -4*ones(n)
    d2 =  ones(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

    return A,B, w_ex

m = 3  # Blocksize

#
# Large scale
#
n = 2500
A,B, w_ex = sakurai(n) # Mikota pair
import numpy
X,_ = numpy.linalg.qr(rand(n,3*m))
data=[]
tt = time.clock()
eigs,vecs, resnh = lobpcg(A,X,B, tol=1e-6, maxiter=500, retResidualNormsHistory=1,largest=True)
from scipy.sparse.linalg import eigs as eeigs
from scipy.sparse.linalg import inv, svds
#eigs2,vecs2 = eeigs(A,m,B)
ZZ=inv(B).dot(A)
ZZZ= ZZ.T.dot(ZZ)
Zs=svds(ZZZ,2499,return_singular_vectors=False)
print Zs #numpy.sqrt(Zs[-m:])
#eigs2,vecs2 = eeigs(ZZ.dot(A))
data.append(time.clock()-tt)
print 'Results by LOBPCG for n='+str(n)
print
print eigs
print
print 'Exact eigenvalues'
print
print w_ex[:m]
print
print 'Eig eigenvalues'
print
#print eigs2[:m]
print
print 'Elapsed time',data[0]
loglog(arange(1,n+1),w_ex,'b.')
xlabel(r'Number $i$')
ylabel(r'$\lambda_i$')
title('Eigenvalue distribution')
show()