Nonsymmetric eigenvalue problem

[asd__asdqwe]

Hello,

I'm trying to compute the right eigenvectors of a Hessenberg matrix. It seems I get the correct results when using the driver routine ?geev, but I can't obtain the same results with the computational routines ?hseqr and ?hsein (which should in the end be faster because there is no need to reduce A into Hessenberg form). Do you see what is wrong in the file dhsein.f ?

$ ifort dgeev.f  -lmkl_core -lmkl_intel_thread -liomp5 -lmkl_intel_lp64 && ./a.out < dhsein.d

 Eigenvector( 3)
  9.8107E-01
  1.2917E-01
  1.2998E-01
  6.2617E-02

$ ifort dhsein.f  -lmkl_core -lmkl_intel_thread -liomp5 -lmkl_intel_lp64 && ./a.out < dhsein.d

 Eigenvector( 3)
 -8.7754E-01
  5.4626E-01
  1.0000E+00
  0.0000E+00

Thank you.

[mecej4]

The MKL documentation for ?HSEQR says: 

      If info = 0 and job = 'E', the contents of h are unspecified on exit.

The implication is that, before calling HSEQR, you should save the matrix H, and use the saved copy in the subsequent call to HSEIN.

Another issue, a minor one, is that different eigenvector routines in Lapack may use different normalization for the eigenvectors, so you should always re-scale the vectors to some consistent convention of your choosing before you make comparisons.

It took me a bit of effort to work this puzzle out, mainly because I had not used these Lapack routines before -- it is not often the case in real world problems that an upper Hessenberg matrix lands on your lap. 

For your convenience, here is (a Lapack-95 version of) the test program that I wrote to track the problem.

program xhess
use lapack95
! read an upper Hessenberg matrix, use HSEQR and HSEIN to find eigenvalues and eigenvectors 
implicit none
integer i,j,n,m,mm,info
real(8), allocatable :: h(:,:),hsav(:,:),wr(:),wi(:),vr(:,:)
logical, allocatable :: select(:)
!
read(*,*)          ! title line
read(*,*) n
allocate (h(n,n),hsav(n,n),wr(n),wi(n),select(n),vr(n,2))
read(*,*)((h(i,j),j=1,n),i=1,n)
hsav=h                                  ! save for use with hsein
call hseqr(h,wr,wi,info=info)
write(*,*)'Info from hseqr = ',info
if(info /= 0)stop
write(*,10)(i,wr(i),wi(i),i=1,n)
10 format('Eigenvalues:',/,(I2,2ES17.9))

select([1,2,4]) = .false.
select(3) = .true.
call hsein(hsav,wr,wi,select,vr=vr,eigsrc='n',info=info)
write(*,*)'Info from hsein = ',info
if(info /= 0)stop
vr(:,1) = vr(:,1)/norm2(vr(:,1))
write(*,20)(i,vr(i,1),i=1,n)
20 format('Eigenvector 3',/,(i3,2x,ES17.9))
end program

What your code dhsein.f did, in effect, was to leave out line-13 and use h in place of hsav on line 22.