subspace_extraction_ss Subroutine

    subroutine subspace_extraction_ss(this, basis_index)

        type(davidson_ss), intent(inout) :: this
        integer, intent(in) :: basis_index
        integer :: lwork, info
        real(dp), allocatable, dimension(:) :: work
        real(dp) :: eigenvalue_list(basis_index)
        ! these are not stack-allocated since they are not always used
        real(dp), allocatable :: eigenvalue_list_imag(:)
        real(dp), allocatable :: left_eigenvectors(:, :)
        real(dp), allocatable :: right_eigenvectors(:, :)
        integer :: minInd, tmp(1)

        if (t_non_hermitian_2_body) then

            lwork = max(1, 4 * basis_index)
            allocate(work(lwork))
            allocate(eigenvalue_list_imag(basis_index))
            ! we are only interested in the right eigenvectors, so left are not referenced
            allocate(left_eigenvectors(1, basis_index))
            allocate(right_eigenvectors(basis_index, basis_index))
            ! call the general lapack routine
            call dgeev( &
                'N', &
                'V', &
                basis_index, &
                this%projected_hamil_work(1:basis_index, 1:basis_index), &
                basis_index, &
                eigenvalue_list, &
                eigenvalue_list_imag, &
                left_eigenvectors, &
                1, &
                right_eigenvectors, &
                basis_index, &
                work, &
                lwork, &
                info &
                )

            ! the eigenvalues do not come out sorted the way we would like, so get the
            ! target's index
            minInd = sum(minloc(eigenvalue_list))

            ! store the davidson vector
            this%davidson_eigenvalue = eigenvalue_list(minInd)
            this%eigenvector_proj(1:basis_index) = right_eigenvectors(1:basis_index, minInd)

            deallocate(left_eigenvectors)
            deallocate(eigenvalue_list_imag)
            deallocate(right_eigenvectors)
        else
            ! 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%projected_hamil_work(1:basis_index, 1:basis_index), &
                basis_index, &
                eigenvalue_list, &
                work, &
                lwork, &
                info &
                )

            this%davidson_eigenvalue = eigenvalue_list(1)
            ! The first column stores the ground state.
            this%eigenvector_proj(1:basis_index) = this%projected_hamil_work(1:basis_index, 1)
        end if
        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_this_proc 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_this_proc 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.
        if (this%space_size_this_proc > 0) then
            call dgemv('N', &
                       this%space_size_this_proc, &
                       basis_index, &
                       1.0_dp, &
                       this%basis_vectors(:, 1:basis_index), &
                       this%space_size_this_proc, &
                       this%eigenvector_proj(1:basis_index), &
                       1, &
                       0.0_dp, &
                       this%davidson_eigenvector, &
                       1)
        end if

    end subroutine subspace_extraction_ss