create eigenvalues with Fortran

[Ingeborg]

Hey,

i have an noob problem. I have the most expierences in c and Matlab and i try to create an own roots function in matlab and started to learn FORTRAn. For this ive written a program which creates the funktion fun(x) = 3x^3 + 6x^2 + 9x + 12 and now i want to get from the companion matrix the eigenvalues which representate the roots from the function.

Code: Select all

program bsp

integer                 :: k = 1, lf = 3,SDIM,info
real                    :: fun(4)
real                    :: compagnion_mtx(3,3)
real                    :: real_eig(3), imag_eig(3)
double precision      :: work(9)
fun = (/3,6,9,12/)
compagnion_mtx = 0;
compagnion_mtx(1,:) = fun(2:4)/fun(1);
compagnion_mtx(2,1) = 1;
compagnion_mtx(3,2) = 1;

do
    if( k >= 4) exit
    write(*,*) compagnion_mtx(k,:)
    k = k+1;
end do
k = 1;

write(*,*), "test"

call DGEES ( 'N', 'N', 'N', lf, compagnion_mtx, lf, SDIM, real_eig, imag_eig, 'N', work, 9, 'N', info)

write(*,*), "test"

do
    if( k >= 4) exit
    write(*,*) real_eig(k),imag_eig(k)
    k = k+1;
end do
k = 1;

end


I have followed this documentation and filled parts of the input which wouldnt refered simply with 'N'. I am not sure what is wrong. It compiles something and then it crashes and finishes the execution.

I hope you can help me, thanks.

Inge

[admin]

Hi below is your code with a few modifications:
1) SGEES has to be use if you use real precision (DGEES for double precision)
2) you cannot use the 'N' if you want to skip the parameters. Even if they are not reference, you need to pass something.
3) It is recommended to compute the workspace size first with LWORK=-1 then call the routine with the required workspace. You can of course use allocatable array if you wish.

Julie

Code: Select all

program bsp

integer                 :: k = 1, lf = 3,SDIM,info, lwork
real                    :: fun(4)
real                    :: compagnion_mtx(3,3)
real                    :: real_eig(3), imag_eig(3)
real                    :: work(109), VS(1)
logical                 :: SLT, BWORK(1)

fun = (/3,6,9,12/)
compagnion_mtx = 0;
compagnion_mtx(1,:) = fun(2:4)/fun(1);
compagnion_mtx(2,1) = 1;
compagnion_mtx(3,2) = 1;

SDIM=0
write(*,*), "[Data] Compagnion Matrix ",lwork
do
    if( k >= 4) exit
    write(*,*) compagnion_mtx(k,:)
    k = k+1;
end do
k = 1;

call SGEES ( 'N', 'N', SLT, lf, compagnion_mtx, lf, SDIM, real_eig, imag_eig, VS ,1, work, -1, BWORK, info)
lwork=work(1)
write(*,*), "[Computing Workspace] Workspace Size=",lwork

write(*,*), "[Computing Eigenvalues] Calling SGEES"
call SGEES ( 'N', 'N', SLT, lf, compagnion_mtx, lf, SDIM, real_eig, imag_eig, VS ,1, work, lwork, BWORK, info)

write(*,*), "[Solution] Eigenvalues:"

do
    if( k >= 4) exit
    write(*,*) real_eig(k),imag_eig(k)
    k = k+1;
end do
k = 1;

end

[Ingeborg]

Hey Julie,

I just want to say thanks to you.

Ingeborg