polyroots Subroutine

subroutine polyroots(n, p, wr, wi, info)

    implicit none

    integer, intent(in) :: n !! polynomial degree
    real(wp), dimension(n+1), intent(in) :: p !! polynomial coefficients array, in order of decreasing powers
    real(wp), dimension(n), intent(out) :: wr !! real parts of roots
    real(wp), dimension(n), intent(out) :: wi !! imaginary parts of roots
    integer, intent(out) :: info !! output from the lapack solver.
                                 !!
                                 !! * if `info=0` the routine converged.
                                 !! * if `info=-999`, then the leading coefficient is zero.

    integer :: i !! counter
    real(wp), allocatable, dimension(:,:) :: a !! companion matrix
    real(wp), allocatable, dimension(:) :: work !! work array
    real(wp), dimension(1) :: vl, vr !! not used here

#ifdef REAL32
    interface
        subroutine sgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
            implicit none
            character :: jobvl, jobvr
            integer :: info, lda, ldvl, ldvr, lwork, n
            real :: a(lda, *), vl(ldvl, *), vr(ldvr, *), wi(*), work(*), wr(*)
        end subroutine sgeev
    end interface
#elif REAL128
    ! do not have a quad solver in lapack
#else
    interface
        subroutine dgeev(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info)
            implicit none
            character :: jobvl, jobvr
            integer :: info, lda, ldvl, ldvr, lwork, n
            double precision :: a(lda, *), vl(ldvl, *), vr(ldvr, *), wi(*), work(*), wr(*)
        end subroutine dgeev
    end interface
#endif

    ! error check:
    if (abs(p(1)) == 0.0_wp) then
        info = -999
        return
    end if

    ! allocate the work array:
    allocate (work(3*n))

    ! create the companion matrix
    allocate (a(n, n))
    a = 0.0_wp
    do i = 1, n - 1
        a(i, i + 1) = 1.0_wp
    end do
    do i = n, 1, -1
        a(n, n - i + 1) = -p(i + 1)/p(1)
    end do

    ! call the lapack solver:
#ifdef REAL32
    ! single precision
    call sgeev('N', 'N', n, a, n, wr, wi, vl, 1, vr, 1, work, 3*n, info)
#elif REAL128
    error stop 'do not have a quad solver in lapack'
#else
    ! by default, use double precision:
    call dgeev('N', 'N', n, a, n, wr, wi, vl, 1, vr, 1, work, 3*n, info)
#endif

end subroutine polyroots