- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- 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
- 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
- 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
- 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