- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page