Solve the roots of a cubic polynomial using the polynomial companion matrix (depends on LAPACK routine SGEEV)

polynomial_companion_matrix.f90

! Compile with:
! $ gfortran -Wall polynomial_companion_matrix.f90 -o polynomial_companion_matrix -llapack
!
! To run the code:
! $ ./polynomial_companion_matrix
!
program main
    implicit none
    integer, parameter :: sp = kind(1.0)
    integer, parameter :: n = 3
    real(sp) :: a(n), M(n,n),wr(n),wi(n)
    real(sp) :: work(3*n), vl(n),vr(n)
    integer :: k, info
    character(len=20) :: num

    if (command_argument_count() == n) then
        do k = 1, n
            call get_command_argument(k,num)
            read(num,*) a(k)
        end do
        print *, "Coefficients: ",a
    else
        write(*,"(A,I0,A)") "Please enter coefficients for the ", n, "-th order polynomial"
        write(*,"(A)") "(the highest order coefficient is assumed to be 1)"
        read(*,*) a(:)
    end if
    print *, "Polynomial coefficients (a_1,a_2,...):", a

    ! Form the companion matrix
    M = 0
    do k = 1, n-1
        M(k,k+1) = 1
    end do
    do k = n,1,-1
        M(n,n-k+1) = -a(k)
    end do

    print *, "Polynomial companion matrix:"
    do k = 1, n
        print *, M(k,:)
    end do

    ! Call LAPACK routine xGEEV
    call sgeev('N','N',n,M,n,wr,wi,vl,1,vr,1,work,3*n,info)
    if (info /= 0) then
        print *, "[sgeev] info = ", info
        if (info < 0) then
            print *, "The ",abs(info),"-th argument had an illegal value."
        else
            print *, "The QR algorithm failed to compute all the eigenvalues."
        end if
        error stop 1
    end if

    print *, "Real and Imaginary parts of the roots:"
    do k = 1, n
        print *, wr(k), wi(k)
    end do

 end program main