1! 2! Copyright (C) 2003 A. Smogunov 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8subroutine gep_x(n, amt, bmt, eigen, veigen) 9! 10! It solves GEP: A X = lambda B X using LAPACK routines 11! 12 USE kinds, only : DP 13 USE cond, only : delgep 14 implicit none 15 16 integer :: i, n, info, lwork 17 complex(DP) :: trywork 18 real(DP), allocatable :: rwork(:) 19 complex(DP) :: & 20 amt(n,n), & ! A 21 bmt(n,n), & ! B 22 eigen(n), & ! lambda 23 veigen(n,n) ! X 24 complex(DP), allocatable :: alpha(:), beta(:), work(:) 25 26 allocate( alpha( n ) ) 27 allocate( beta( n ) ) 28 allocate( rwork( 8*n ) ) 29! 30! for some reason the lapack routine does not work if the diagonal elements 31! of the matrices are exactly zero. We tested these routines on 32! pc_ifc, ibmsp and origin. 33! 34 do i=1,n 35 amt(i,i)=amt(i,i)+delgep 36 bmt(i,i)=bmt(i,i)+delgep 37 enddo 38 39 call ZGGEV('N', 'V', n, amt, n, bmt, n, alpha, beta, veigen, n, veigen, & 40 n, trywork, -1, rwork, info) 41 42 lwork=abs(trywork) 43 allocate( work( lwork ) ) 44 45 call ZGGEV('N', 'V', n, amt, n, bmt, n, alpha, beta, veigen, n, veigen, & 46 n, work, lwork, rwork, info) 47 48 if(info.ne.0) call errore('gep_x','error on zggev',info) 49! 50! lambda = alpha / beta 51! 52 do i=1, n 53 eigen(i)=alpha(i)/beta(i) 54! write(6,'(i5, 2f40.20)') i, DBLE(eigen(i)), AIMAG(eigen(i)) 55 enddo 56 57 deallocate(work) 58 deallocate(rwork) 59 deallocate(beta) 60 deallocate(alpha) 61 62 return 63end subroutine gep_x 64