Eigenvectors using MKL

Marcos_D_ · ‎04-11-2014

Hi, I need help to calculate Eigenvectors using MKL routines.

I have some codes that were compiled using Compaq Visual Fortran and they had IMSL libraries within. Now I'm using Visual Studio 2010+Visual Fortran XE 2013. And I'm trying to convert those old codes without IMSL, but using MKL (Lapack) functions. In this case is the function DEVCRG ( http://www.roguewave.com/Portals/0/products/imsl-numerical-libraries/fortran-library/docs/6.0/math/default.htm?turl=evcrg.htm ). As you can check they have an example in that page. That's the same matrix i'm trying in the code bellow. I've checked with MATLAB and I've found the same results.

So, here we go with MKL, following the decision tree ( http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-96CDD8D9-37D7-477A-8EEE-7E6F4081E3DF.htm ) in the manual I'm trying in this sequence those routines: ?gebal , ?gehrd , ?unghr , ?hseqr , ?trevc , ?gebak

I've arrived to find the eigenvalues, but I can't find the eigenvectors!!

0.2000E+01 0.4000E+01
0.2000E+01 -.4000E+01
0.1000E+01 -.2297E-15

    program eigenvalvec

    use MKL95_LAPACK
    
    implicit none
    double precision, allocatable :: x(:,:),scalev(:)
    double complex, allocatable :: xc(:,:),tauc(:),w(:),uph(:,:),q(:,:),vr(:,:),vl(:,:),qz(:,:)
    logical, allocatable :: select(:)
    integer :: n = 3.0
    integer :: nouput,i,ilo,ihi,info
    character*14 file10
    allocate(x(n,n),scalev(n),xc(n,n),tauc(n-1),w(n),uph(n,n),q(n,n),select(n),vr(n,n),vl(n,n),qz(n,n))
    
    nouput=1
    
    x(1,1)=8.0 ; x(1,2)=-1.0 ; x(1,3)=-5.0
    x(2,1)=-4.0 ; x(2,2)=4.0 ; x(2,3)=-2.0
    x(3,1)=18.0 ; x(3,2)=-5.0 ; x(3,3)=-7.0
    
    xc = CMPLX(x)
   
    call gebal(xc,scalev,ilo,ihi,'P',info) !Balances a general matrix to improve the accuracy of computed eigenvalues and eigenvectors
	
    uph = xc
    call gehrd(uph,tauc,ilo,ihi,info) !Reduces a general matrix to upper Hessenberg form

    q = uph
    call unghr(q,tauc,ilo,ihi,info) !Generates the complex unitary matrix Q determined by ?gehrd.
    
    call hseqr(uph,w,ilo,ihi,q,'S','I',info) !Computes all eigenvalues and (optionally) the Schur factorization of a matrix reduced to Hessenberg form.
    !w: EIGENVALUES!!!
    
    vl = q
    vr = q
    call trevc(uph,'B',select,vl,vr,n,info) !Computes selected eigenvectors of an upper (quasi-)triangular matrix computed by ?hseqr.
    
    call gebak(vr,scalev,ilo,ihi,'B','R',info) !Transforms eigenvectors of a balanced matrix to those of the original nonsymmetric matrix.

    end program eigenvalvec

Please, I really need some light on that... specially with "?trevc". Thank you very much!!!

Marcos_D_ · ‎04-13-2014

Someone told me to skip all that trouble and use the call GEEVX.

I feel I'm almost there to find the eigenvectors exactly like DEVCRG (from IMSL) and MATLAB. In that case:

0.3162 + 0.3162i 0.3162 - 0.3162i 0.4082

-0.0000 + 0.6325i -0.0000 - 0.6325i 0.8165

0.6325 0.6325 0.4082

Well, using geevx and setting "balanc" as 'N' or 'P' I find:

0.3162 0.3162 0.3162 0.3162 0.4082 -0.0000

-0.0000 0.6325 0.6325 0.0000 0.8165 0.0000

0.6325 0.0000 -0.0000 0.6325 0.4082 -0.0000

Now setting "balanc" as 'S' or 'B' I get:

0.3162 -0.3162 0.3162 -0.3162 0.4082 0.0000

0.6325 0.0000 -0.0000 -0.6325 0.8165 0.0000

-0.0000 -0.6325 0.6325 0.0000 0.4082 0.0000

I couldn't set the variable "abnrm". I can't set the right type of it. But anyway, I don't think it will affect the result, since it is only an output. Here is the new code I have:

program eigsreal

use MKL95_LAPACK

implicit none
double precision, allocatable :: x(:,:),scalev(:),rconde(:),rcondv(:)
double complex, allocatable :: xc(:,:),w(:),vr(:,:),vl(:,:)
integer :: n = 3.0
integer :: nouput,iochek,i,ilo,ihi,info
character*14 file10
allocate(x(n,n),scalev(n),xc(n,n),w(n),vr(n,n),vl(n,n),rconde(n),rcondv(n))

nouput=1

x(1,1)=8.0 ; x(1,2)=-1.0 ; x(1,3)=-5.0
x(2,1)=-4.0 ; x(2,2)=4.0 ; x(2,3)=-2.0
x(3,1)=18.0 ; x(3,2)=-5.0 ; x(3,3)=-7.0

xc = CMPLX(x)

!pg 1202
!?geevx - Computes the eigenvalues and left and right eigenvectors of a general matrix, with preliminary
!         matrix balancing, and computes reciprocal condition numbers for the eigenvalues and right eigenvectors.

!call geevx(a, w[,vl] [,vr] [,balanc] [,ilo] [,ihi] [,scale] [,abnrm] [,rconde] [,rcondv] [,info])

call geevx(xc,w,vl,vr,'P',ilo,ihi,scalev)

file10='geevx_vr.out'
open(unit=nouput,file=file10,status='unknown',iostat=iochek,form='formatted')
do i = 1,n
    write(nouput,'(36(F10.4,1x))') vr(i,:)
end do

end program eigsreal