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
$ ifort dhsein.f -lmkl_core -lmkl_intel_thread -liomp5 -lmkl_intel_lp64 && ./a.out < dhsein.d
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.