2
$\begingroup$

exchangers,

I have run into a bit of a puzzling problem. To solve an complex eigenvalue-problem, I make use of the LAPACK library function ZHEEV. To test the implementation I used a real symmetric matrix:

\begin{align*} \left[ \begin{array}{c c c c} 2 & 0 & 2 & 0 \\ 0 & 2 & 0 & 2 \\ 2 & 0 & 2 & 0 \\ 0 & 2 & 0 & 2 \\ \end{array} \right] \end{align*}

which should have eigenvalues: 2x0, 2x4 Using ZHEEV I get however 2x -2 and 2x 2. Modifying the implementation to use DSYEV on the other hand gives the correct results.

I am missing something very trivial here, but I do not see what it is.

This is the subroutine:

subroutine solvearray(m)                                                                                         
complex(16), dimension(:,:),allocatable :: dm
complex(16), dimension(:), allocatable :: freq, work
complex(16) :: One
doubleprecision, dimension(:),allocatable :: w, rwork
integer,dimension(:,:),allocatable, intent(inout) :: m 
integer :: n, k, h, l, info

One = cmplx(1.0,0.0,kind=SELECTED_REAL_KIND(15,300))
write(*,'(A)',advance='no') 'order of matrix = '   
read(*,*) n  
l = 2*n - 1 
allocate(work(2*l), rwork(3*n-2), w(n))  
w = 0.0  
work = 0.0 *One 
rwork = 0.0  
info = 0
allocate(dm(1:n,1:n))
dm(1:n,1:n) = m(1:n,1:n)*One
write(*,*) 'Setup for lapack workspace has been completed.'
call ZHEEV('V','U', n, dm, n, w, work, 2*l, rwork, info)

write(*,*) 'This is the eigenvector matrix:'
do k=1,n  
    write(*,'(255F8.4)') dm(k,1:n)
end do 
write(*,*) 'These are the eigenvalues :'
write(*,'(255F8.4)') w  
deallocate(m)           
deallocate(W,Work,rwork)
end subroutine
Improve this question
$\endgroup$
7
  • $\begingroup$ Your expectation that zheev should give the correct results is reasonable. You likely have a bug in your code or possibly one in your lapack library. $\endgroup$
    – Bill Greene
    Mar 26, 2017 at 13:22
  • 1
    $\begingroup$ It is impossible to answer this question without seeing the exact code you have. Knowing that you use ZHEEV just isn't much to go on. $\endgroup$
    – Kirill
    Mar 26, 2017 at 17:25
  • $\begingroup$ I have added the subroutine (removed most of the things I included during debugging, to make it a bit shorter) $\endgroup$
    – DannyVanpoucke
    Mar 28, 2017 at 14:50
  • $\begingroup$ hmm, I seem to have found the problem. Apparently complex*16 is not complex(16) but complex(8) instead, which unfortunately doesn't make ZHEEV trip up. link $\endgroup$
    – DannyVanpoucke
    Mar 28, 2017 at 14:58
  • 2
    $\begingroup$ @DannyVanpoucke As the guy who provided the answer you linked to can I reiterate what it essentially says there. You should never, Never, NEVER use Real(8) , Complex(16) and similar in Fortran. Please read through, understand why, and learn the right way to do it. $\endgroup$
    – Ian Bush
    Mar 28, 2017 at 16:25

1 Answer 1

Reset to default
2
$\begingroup$

Following a similar problem answered by Ian Bush , the problem is solved by knowing that complex(16) is not the same as complex*16. In the specific case of the 64bit compiler used complex*16 needed to be translated to complex(8)...or even better use SELECT_REAL_KIND to define a compiler independent definition: C_double=SELECTED_REAL_KIND(15,307).

Thank you all for the help and useful remarks.

Danny

Improve this answer
$\endgroup$
1
  • $\begingroup$ Just to add complex*16 is not and has never been standard Fortran, and there is absolutely no advantage in using it compared to the standard and, if used properly, completely portable kind mechanism in the language. $\endgroup$
    – Ian Bush
    Apr 17, 2017 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.