Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Solve the roots of a cubic polynomial using the polynomial companion matrix (depends on LAPACK routine SGEEV)
! 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment