subspace_extraction Subroutine

    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