from scipy import *
from scipy.sandbox import lobpcg
from scipy.sparse import spdiags, speye
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 = speye( 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)

    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 = 25000
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)
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 'Elapsed time',data[0]
loglog(arange(1,n+1),w_ex,'b.')
xlabel(r'Number $i$')
ylabel(r'$\lambda_i$')
title('Eigenvalue distribution')
show()