implicitly restart Lanczos method

Question

I want to write simple toy code for implicitly restart Lanczos method. Without implicit restarting, the code is perfectly working but when I turn on the restart, I cannot get proper solution

To my knowledge, newly constructed w should be orthogonal to all of the new Lanczos vectors. For the first restart, the orthogonality is well preserved but from the second restart, the orthogonality is significantly broken down and the program does not find proper eigenvalues.

I already spent several tens of hours to fix it. I almost gave up...... Here is my python code

"""
Author: Sunghwan Choi
Date Created: June 19, 2017
Python Version: 2.7 or 3.5

Reference for Lanczos algorithm
http://www.netlib.org/utk/people/JackDongarra/etemplates/node104.html
Reference for implicit restart
http://www.netlib.org/utk/people/JackDongarra/etemplates/node118.html
"""

import numpy as np
from scipy.sparse.linalg import eigsh
#from scipy.sparse import eye
from scipy.sparse import coo_matrix
from numpy import eye

def clustering(eigvals,eigvecs,tol=1e-2):
    ret_eigvals=[]
    ret_eigvecs=[]
    for i in range(len(eigvals)):
        for ret_eigval, ret_eigvec in zip (ret_eigvals,ret_eigvecs):
            if (np.abs(eigvals[i]/ret_eigval-1.0)<tol ):
                break
        else:
            ret_eigvals.append(eigvals[i])
            ret_eigvecs.append(eigvecs[:,i])
    ret_eigvals=np.array(ret_eigvals)
    ret_eigvecs=np.array(ret_eigvecs).T
    return ret_eigvals,ret_eigvecs

def check_conv(matrix, cal_eigval, cal_eigvec, tol):
    indices=[]
    for i in range(len(cal_eigval)):
        if(np.linalg.norm(matrix.dot(cal_eigvec[:,i]) - cal_eigval[i]*cal_eigvec[:,i])< tol):
            indices.append(i)
    return indices

################ input
size=1600
max_step=20000
which='SA'
#implicit=False
implicit=True
energy_range=[0.0,6.0]
tol = 1e-5 
n_eig=6
n_tol_check=40 # n_tol_check>n_eig ==0
######################

# generate 1D harmonic oscillator
h=0.1
matrix=-5/2*eye(size)
matrix+=4/3*(eye(size,k=1)+eye(size,k=-1))
matrix+=-1/12*(eye(size,k=2)+eye(size,k=-2))
matrix=-0.5*matrix/(h*h)
distance =lambda index: (index-size/2)*h
matrix+=np.diagflat( list(map( lambda i: 0.5*distance(i)**2, range(size))))

# solve eigenvalue problem to check validity
true_eigval,true_eigvec = eigsh(matrix,k=50,which=which)
indices = np.all([true_eigval>energy_range[0], true_eigval<energy_range[1]],axis=0)
true_eigval = true_eigval[indices]
true_eigvec = true_eigvec[:,indices]

#initialize variables
alpha=[]; beta=[]
index_v=0
restart_interval = n_tol_check+n_eig if implicit is not False else max_step
T = np.zeros((restart_interval,restart_interval))
v = np.zeros((size,restart_interval))
#Q=np.eye(restart_interval)

#generate initial vector

np.random.seed(1)
initial_vec = np.random.random(size)

#initial_vec = np.loadtxt("tmp")
w = v[:,index_v] = initial_vec/np.linalg.norm(initial_vec)
init_beta = np.linalg.norm(w)
# start Lanczos i_step
for i_step in range(max_step):
    if (i_step is 0):
        v[:,index_v] = w/init_beta
    else:
        v[:,index_v] = w/T[index_v,index_v-1]
    w=matrix.dot(v[:,index_v])

    if (i_step is 0):
        w=w-init_beta*v[:,index_v-1]
    else:
        w=w-T[index_v,index_v-1]*v[:,index_v-1]
    T[index_v,index_v]=np.dot(w,v[:,index_v])
    w -=T[index_v,index_v]*v[:,index_v]

    #check convergence
    if ((i_step+1)%n_tol_check==n_eig and i_step>n_eig):
        # calculate eigenval of T matrix
        cal_eigval, cal_eigvec_= np.linalg.eigh(T[:index_v+1,:index_v+1])
        cal_eigvec = np.dot(v[:,:index_v+1],cal_eigvec_)
        #check tolerance
        conv_indices = check_conv(matrix, cal_eigval, cal_eigvec,tol)
        #filter energy_range
        indices = np.all([cal_eigval[conv_indices]>energy_range[0], cal_eigval[conv_indices]<energy_range[1]],axis=0)
        #check clustering
        conv_cal_eigval,conv_cal_eigvec = clustering((cal_eigval[conv_indices])[indices], (cal_eigvec[conv_indices])[indices])
        if (len(conv_cal_eigval)>=n_eig):
            break
        # implicit restarting 
        if (implicit is True):
            Q=np.eye(restart_interval)
            # do shift & QR decomposition
            indices = np.argsort(np.abs(cal_eigval-np.mean(energy_range)))
            for index in indices[n_eig:]:
                new_Q,new_R = np.linalg.qr(T-cal_eigval[index]*np.eye(len(T)))
                T = np.dot(new_Q.T,np.dot(T,new_Q))
                v = np.dot(v,new_Q)
                Q = np.dot(Q,new_Q)

            w=v[:,n_eig]*T[n_eig,n_eig-1]+w*Q[-1,n_eig-1]
            v[:,n_eig:]=0.0

            T[:,n_eig:] = 0.0
            T[n_eig:,:] = 0.0
            #for debug
            #print(np.dot(w.T, v))

            # reset index 
            index_v=n_eig-1
    index_v+=1
    T[index_v,index_v-1]=np.linalg.norm(w)
    T[index_v-1,index_v]=np.linalg.norm(w)
else:
    print("not converged")
    exit(-1)

print ("energy window: (", energy_range[0],",",energy_range[1],")")
print ("true eigenvalue")
print(true_eigval)
print ("eigenvalue from Lanczos w/ implicit restart (",i_step+1,")")
print(conv_cal_eigval)