subspace_extraction_sl Subroutine

    subroutine subspace_extraction_sl()

        integer :: lwork, info, i
        real(dp), allocatable :: work(:), left_pert_overlaps(:), left_pert_overlaps_all(:)

        ! Scrap space for the diagonaliser.
        lwork = max(1, 3 * n_lanc_vecs_sl - 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).
        ! n_lanc_vecs_sl is the number of rows and columns in A.
        ! A = sl_hamil. This matrix stores the eigenvectors in its columns on output.
        ! n_lanc_vecs_sl is the leading dimension of A.
        ! sl_h_eigv stores the eigenvalues on output.
        ! work is scrap space.
        ! lwork is the length of the work array.
        ! info = 0 on output if diagonalisation is successful.
        call dsyev('V', 'U', n_lanc_vecs_sl, sl_hamil, n_lanc_vecs_sl, sl_h_eigv, work, lwork, info)

        deallocate(work)

        ! The first Lanczos vector differs from the perturbed ground state only
        ! by a normalisation factor, which is accounted for here.
        trans_amps_right = sl_hamil(1, :) * right_pert_norm / sqrt(pops_norm)

        ! Calculate the overlaps between the left-perturbed ground state and all
        ! of the Lanczos vectors.
        allocate(left_pert_overlaps(n_lanc_vecs_sl))
        allocate(left_pert_overlaps_all(n_lanc_vecs_sl))
        left_pert_overlaps = matmul(pert_ground_left, sl_vecs)
        ! Sum the contributions to the overlaps from all processes.
        call MPISumAll(left_pert_overlaps, left_pert_overlaps_all)
        ! And then use these overlaps to calculate the overlaps of the
        ! left-perturbed ground state with the final eigenvectors.
        trans_amps_left = matmul(left_pert_overlaps_all, sl_hamil)
        deallocate(left_pert_overlaps)
        deallocate(left_pert_overlaps_all)

    end subroutine subspace_extraction_sl