using LinearAlgebra ## HELPERS function orthonormalize(V::Matrix{Float64}) Q,R = qr(V) return Q*Matrix{Float64}(I,size(V)...) end function eachcolumn(V::AbstractMatrix) [V[:,i] for i in 1:size(V,2)] end ## THE ALGORITHM function Davidson(A::AbstractMatrix, NTargetEigvals::Int=1, subspaceincrement::Int=8, maxsubspacesize::Int=200, ϑ::Float64=1e-8) k = subspaceincrement M = min(maxsubspacesize, size(A,2)) V = Matrix{Float64}(I,size(A,1),k) eigvals = ones(NTargetEigvals) for m in k:k:M V = orthonormalize(V) AV = A * V T = Symmetric(V' * AV) eig = eigen(T) w = [AV*y - λ*(V*y) for (λ,y) in zip(eig.values[1:k],eachcolumn(eig.vectors)[1:k])] V = [V w...] eigvalsN = eig.values[1:NTargetEigvals] eigvals = norm(eigvalsN-eigvals)>=ϑ ? eigvalsN : return eig.values end error("Davidson did not converge in $maxsubspacesize-dimensional subspace") end ## THE PROBLEM testmatrix= Symmetric(diagm(0=>1:1200) + 0.0001*rand(1200,1200)) println("Exact eigenvalues: ",eigen(testmatrix).values[1:4]) @time eigen(testmatrix).values; println("Davidson eigenvls: ",Davidson(testmatrix,4)[1:4]) @time Davidson(testmatrix,4); nothing