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