subroutine subspace_expansion(this, basis_index)
type(DavidsonCalcType), intent(inout) :: this
integer, intent(in) :: basis_index
integer :: i
real(dp) :: dot_prod, norm
! Create the new basis state from the residual. This step performs
! t = (D - EI)^(-1) r,
! where D is the diagonal of the Hamiltonian matrix, E is the eigenvalue previously
! calculated, I is the identity matrix and r is the residual.
do i = 1, this%super%space_size
this%super%basis_vectors(i, basis_index) = this%residual(i) / (hamil_diag(i) - this%davidson_eigenvalue)
end do
! This step then maskes the new basis vector orthogonal to all other basis vectors, by doing
! t <- t - (t,v)v
! for each basis vector v, where (t,v) denotes the dot product.
do i = 1, basis_index - 1
dot_prod = dot_product(this%super%basis_vectors(:, basis_index), this%super%basis_vectors(:, i))
this%super%basis_vectors(:, basis_index) = &
this%super%basis_vectors(:, basis_index) - dot_prod * this%super%basis_vectors(:, i)
end do
! Finally we calculate the norm of the new basis vector and then normalise it to have a norm of 1.
! The new basis vector is stored in the next available column in the basis_vectors array.
norm = dot_product(this%super%basis_vectors(:, basis_index), this%super%basis_vectors(:, basis_index))
norm = sqrt(norm)
this%super%basis_vectors(:, basis_index) = this%super%basis_vectors(:, basis_index) / norm
end subroutine subspace_expansion