# subspace_extraction Subroutine

## public subroutine subspace_extraction(this, basis_index)

### Arguments

Type IntentOptional Attributes Name
type(DavidsonCalcType), intent(inout) :: this
integer, intent(in) :: basis_index

## Source Code

    subroutine subspace_extraction(this, basis_index)

type(DavidsonCalcType), intent(inout) :: this
integer, intent(in) :: basis_index
integer :: lwork, info
real(dp), allocatable, dimension(:) :: work
real(dp) :: eigenvalue_list(basis_index)
! Scrap space for the diagonaliser.
lwork = max(1, 3 * basis_index - 1)
allocate(work(lwork))

! This routine diagonalises a symmetric matrix, A.
! V tells the routine to calculate eigenvalues *and* eigenvectors.
! U tells the routine to get the upper half of A (it is symmetric).
! basis_index is the number of rows and columns in A.
! A = projected_hamil_work. This matrix stores the eigenvectors in its columns on output.
! basis_index is the leading dimension of A.
! eigenvalue_list stores the eigenvalues on output.
! work is scrap space.
! lwork is the length of the work array.
! info = 0 on output is diagonalisation is successful.
call dsyev( &
'V', &
'U', &
basis_index, &
this%super%projected_hamil_work(1:basis_index, 1:basis_index), &
basis_index, &
eigenvalue_list, &
work, &
lwork, &
info &
)

! The first column stores the ground state.
this%davidson_eigenvalue = eigenvalue_list(1)
this%eigenvector_proj(1:basis_index) = this%super%projected_hamil_work(1:basis_index, 1)
deallocate(work)

! eigenvector_proj stores the eigenstate in the basis of vectors stored in the array
! basis_vectors. We now want it in terms of the original basis. To get this, multiply
! eigenvector_proj by basis_vectors(:, 1:basis_index):

! This function performs y := alpha*a*x + beta*y
! N specifies not to use the transpose of a.
! space_size is the number of rows in a.
! basis_index is the number of columns of a.
! alpha = 1.0_dp.
! a = basis_vectors(:,1:basis_index).
! space_size is the first dimension of a.
! input x = eigenvector_proj(1:basis_index).
! 1 is the increment of the elements of x.
! beta = 0.0_dp.
! output y = davidson_eigenvector.
! 1 is the incremenet of the elements of y.
call dgemv('N', &
this%super%space_size, &
basis_index, &
1.0_dp, &
this%super%basis_vectors(:, 1:basis_index), &
this%super%space_size, &
this%eigenvector_proj(1:basis_index), &
1, &
0.0_dp, &
this%davidson_eigenvector, &
1)

end subroutine subspace_extraction