1
$\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   
$\endgroup$
  • $\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 '17 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 '17 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 '17 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 '17 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 '17 at 16:25
1
$\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

$\endgroup$
  • $\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 '17 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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