Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## extended eigensolver routines Beginner
152 Views

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.

1 Solution Employee
152 Views

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

` `

2 Replies Employee
153 Views

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

` ` Beginner
152 Views

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... 