Incorrect answer using sgeev compared to matlab

Question

I was wondering if anyone has any experience using sgeev to compute e-vale/e-vecs in fortran. I am currently having an issue diagonalising a matrix and am not sure why. the matrix is

   1    2   4   4   22  -3  22
   3    3   8   -3  -22 -2  14
   8  -2.3  16  2.5 22  1   7
  -6   17   22  -9  22  17  -6
   7   1    22  2.5 16  -2.3 8
  14  -2   -22  -3  8   3   3
  22    -3  22  4   4   2   1

It is definitely diagonalisable as it works fine in matlab but i can't get it to work in fortran and have no idea why,

my calling for using sgeev is correct as it has been tested out with other matrices and it works fine, returning scaled results etc

I know the matrix has the property that the first column is the inverse of the last column etc but i thought using the general matrix form in fortran would be fine. If anyone could shed any light on this it would be very much appreciated

    'program trial
    implicit none
    integer, parameter :: M=7
    integer, parameter :: N=6
    real :: qqq(7,7), ttt(7,7)
    character*1 :: jobvl, jobvr
    real :: wr(M)
    real :: diag(M,M)
    real :: wi(M)
    real :: vl(M,M)
    integer :: LDVL=M
    integer :: IHI, ILO
    real :: vr(M,M)
    integer :: LDVR=M
    real :: work(4*M)
    integer :: lwork=4*M
    integer :: info, infonow, check
    character (len=40) :: print_file
    integer :: filenumber=1
    integer :: r, rr, rrr

    qqq(1,1)=1
    qqq(1,2)=2
    qqq(1,3)=4
    qqq(1,4)=4
    qqq(1,5)=22
    qqq(1,6)=-3
    qqq(1,7)=22

    qqq(2,1)=3
    qqq(2,2)=3
    qqq(2,3)=8
    qqq(2,4)=-3
    qqq(2,5)=-22
    qqq(2,6)=-2
    qqq(2,7)=14

    qqq(3,1)=8
    qqq(3,2)=-2.3
    qqq(3,3)=16
    qqq(3,4)=2.5
    qqq(3,5)=22
    qqq(3,6)=1
    qqq(3,7)=7

    qqq(4,1)=-6
    qqq(4,2)=17
    qqq(4,3)=22
    qqq(4,4)=-9
    qqq(4,5)=22
    qqq(4,6)=17
    qqq(4,7)=-6

    qqq(5,1)=7
    qqq(5,2)=1
    qqq(5,3)=22
    qqq(5,4)=2.5
    qqq(5,5)=16
    qqq(5,6)=-2.3
    qqq(5,7)=8

    qqq(6,1)=14
    qqq(6,2)=-2
    qqq(6,3)=-22
    qqq(6,4)=-3
    qqq(6,5)=8
    qqq(6,6)=3
    qqq(6,7)=3

    qqq(7,1)=22
    qqq(7,2)=-3
    qqq(7,3)=22
    qqq(7,4)=4
    qqq(7,5)=4
    qqq(7,6)=2
    qqq(7,7)=1


    do rr=1,7
    do r=1,7
    ttt(r,rr)=qqq(r,rr)
    end do
    end do

    jobvl='V'
    jobvr='V'

    call SGEEV(jobvl,jobvr,M,qqq,M,wr,wi,vl,LDVL,vr&
    ,LDVR,work,lwork,info)

    do rr=1,N+1
    do r=1,N+1
    if (r==rr) then
    diag(r,rr)=wr(r)
    else
    diag(r,rr)=0
    end if
    end do
    end do


    write(*,*) wr




    vl=transpose(vl)`

You say you "can't get it to work in fortran", but what does that mean? It doesn't compile, gives wrong numbers, calls you names? — francescalus, Sep 09 '14 at 09:26
Everything compiles fine but it gives a different answer to the one in matlab, when i compute vr*diag(wr)*(vl') it does not equal the original matrix qqq while in matlab it does — weddle_32, Sep 09 '14 at 13:01
*Does not equal* in the sense that `2 != 3` or in the sense that `2.23456789 != 2.23456788` ? — High Performance Mark, Sep 09 '14 at 15:10

score 2 · Accepted Answer · answered Sep 09 '14 at 16:37

There are several problems, which have to be solved:

1) The last four eigenvalues of your matrix are complex. As you are ignoring them, you can't get the correct result.

2) The eigenvectors to conjugate eigenvalues are also complex and are a combination of the reported eigenvectors of the eigenvalues.

Both problems can be solved by switching to the complex variant CGEEV.

The last thing is, that the diagonal form qqq = vr' * diag(wr) * vr only holds in general, if your matrix qqq is hermitian or real symmetric.

In your case, you have to calculate the inverse of the matrix vr. Another possibility is the generation of an orthonormal system from your eigenvectors, which matlab might do for you as default.

Incorrect answer using sgeev compared to matlab

1 Answers1