subroutine subspace_expansion_lanczos(ivec, lanc_vecs, full_vec, proj_hamil, counts, disps)
integer, intent(in) :: ivec
real(dp), intent(inout) :: lanc_vecs(:, :)
real(dp), intent(inout) :: full_vec(:)
real(dp), intent(inout) :: proj_hamil(:, :)
integer(MPIArg), intent(in) :: counts(0:nProcessors - 1), disps(0:nProcessors - 1)
integer :: i, j
real(dp) :: overlap, tot_overlap, last_norm, norm, tot_norm
call MPIAllGatherV(lanc_vecs(:, ivec), full_vec, counts, disps)
! Multiply the last Lanczos vector by the Hamiltonian.
lanc_vecs(:, ivec + 1) = 0.0_dp
do i = 1, counts(iProcIndex)
do j = 1, sparse_ham(i)%num_elements
lanc_vecs(i, ivec + 1) = lanc_vecs(i, ivec + 1) + &
sparse_ham(i)%elements(j) * full_vec(sparse_ham(i)%positions(j))
end do
end do
overlap = dot_product(lanc_vecs(:, ivec + 1), lanc_vecs(:, ivec))
call MPISumAll(overlap, tot_overlap)
if (ivec == 1) then
lanc_vecs(:, ivec + 1) = lanc_vecs(:, ivec + 1) - tot_overlap * lanc_vecs(:, ivec)
else
last_norm = proj_hamil(ivec, ivec - 1)
lanc_vecs(:, ivec + 1) = lanc_vecs(:, ivec + 1) - tot_overlap * lanc_vecs(:, ivec) - last_norm * lanc_vecs(:, ivec - 1)
end if
norm = 0.0_dp
do i = 1, counts(iProcIndex)
norm = norm + lanc_vecs(i, ivec + 1) * lanc_vecs(i, ivec + 1)
end do
call MPISumAll(norm, tot_norm)
tot_norm = sqrt(tot_norm)
! The final new Lanczos vector.
lanc_vecs(:, ivec + 1) = lanc_vecs(:, ivec + 1) / tot_norm
! Update the subspace Hamiltonian. It can be shown that this has a tridiagonal form
! where the non-zero elements are equal to the following.
proj_hamil(ivec, ivec) = tot_overlap
proj_hamil(ivec, ivec + 1) = tot_norm
proj_hamil(ivec + 1, ivec) = tot_norm
end subroutine subspace_expansion_lanczos