topic Hello Irina, in Intel® oneAPI Math Kernel Library

extended eigensolver routines

Marios_G_1 — Fri, 24 Jun 2016 15:05:07 GMT

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.

Hi Marios,

Irina_S_Intel — Sun, 26 Jun 2016 09:17:00 GMT

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

Hello Irina,

Marios_G_1 — Sun, 26 Jun 2016 15:27:19 GMT

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