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

Hello,

I am trying to diagonalize a matrix using the feast algorithm routines. I have written this small program to test

zfeast_heev. The following code tries to diagonalize a 4x4 complex Hermitian matrix:

program test implicit none integer, parameter :: n=4 ! h0 nxn matrix complex(16) :: h0(n,n), zero character(len=1), parameter :: UPLO='F' integer :: fpm(128), loop, M0,M,info real(8) :: epsout, E(n), res(n), Pi, alfa, Mk integer, parameter :: L=4 complex(16) :: X(n,n) integer :: i,j,k,i1,j1 real(8), parameter :: Emin=-5.0d0, Emax=5.0d0 complex(16), parameter :: ii=(0.0d0,1.0d0) open (1, file='eigenvalues.dat') open (2, file='check.dat') zero=dcmplx(0.d0, 0.d0) h0=zero h0(1,2)= 2.0d0 + 2.0d0*ii h0(1,3)= 3.0d0 - 2.0d0*ii h0(2,1)= 2.0d0 - 2.0d0*ii h0(2,4)= 3.0d0 - 2.0d0*ii h0(3,1)= 3.0d0 + 2.0d0*ii h0(3,4)= -2.0d0 - 2.0d0*ii h0(4,2)= 3.0d0 + 2.0d0*ii h0(4,3)= -2.0d0 - 2.0d0*ii M0=L M=M0 print *,'Search interval ', Emin,' ', Emax call feastinit(fpm) fpm(1)=1 print *, ' Testing zfeast_hcsrev ' call zfeast_heev(UPLO,n,h0,n, fpm, epsout, loop,Emin,Emax,M0,E,X,M,res,info) print *,' FEAST OUTPUT INFO ',info if(info/=0) stop 1 print *, 'Number of eigenvalues found ', M do i1=1,M write(1,*) E(i1) end do end program test

The eigenvalues of the matrix h0 are 4.58258, 4.58258, -4.58258, -4.58258 but this program gives me different results.

I hope someone can explain me what I am missing in order to make this code work properly.

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

Hi Marios,

The first problem I see is that matrix h0 is not Hermitian (`h0(3,4) is not equal to `

h0(4,3) conjugate). Intel MKL Extended Eigensolver supports only Symmetric or Hermitian matrices. Set h0(4,3)= -2.0d0 + 2.0d0*ii and that will give you matrix with described eigenvalues.

The second thing is that you use complex`(16) data type. By specifying parameter in a complex(16) type declaration, you specify both real and imaginary part 16-bytes long that gives overall size of 32 bytes. However, `

COMPLEX*16 type `required for `

zfeast_heev is of size 16, ` i.e. imaginary and real part both of size 8. So correct declaration of complex variables in your program would be: `

complex`(8) or `

complex*16.

Best regards,

Irina

` `

Link Copied

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

Hi Marios,

The first problem I see is that matrix h0 is not Hermitian (`h0(3,4) is not equal to `

h0(4,3) conjugate). Intel MKL Extended Eigensolver supports only Symmetric or Hermitian matrices. Set h0(4,3)= -2.0d0 + 2.0d0*ii and that will give you matrix with described eigenvalues.

The second thing is that you use complex`(16) data type. By specifying parameter in a complex(16) type declaration, you specify both real and imaginary part 16-bytes long that gives overall size of 32 bytes. However, `

COMPLEX*16 type `required for `

zfeast_heev is of size 16, ` i.e. imaginary and real part both of size 8. So correct declaration of complex variables in your program would be: `

complex`(8) or `

complex*16.

Best regards,

Irina

` `

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

Hello Irina,

Regarding the elements h0(3,4), h0(4,3) you are right. It was a misprint. I corrected it.

And yes you are correct about the complex declaration!

I changed it to complex(8) and everything works fine!

Thank you,

Marios

P.S: I can't believe that was the problem, I completely missed it...

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