subroutine p_zeigh_simple(A,W,Nblock,method,jobz,uplo,vl,vu,il,iu,tol)
complex(8),dimension(:,:),intent(inout) :: A ! M v = E v/v(i,j) = ith component of jth vec.
real(8),dimension(size(A,2)),intent(inout) :: W ! eigenvalues
integer :: Nblock
integer :: Nb
character(len=*),optional :: method
character(len=1),optional :: jobz,uplo
character(len=1) :: jobz_,uplo_,range
character(len=20) :: method_
real(8),optional :: vl,vu,tol
integer,optional :: il,iu
real(8) :: vL_,vU_,tol_
integer :: iL_,iU_
integer :: Ns
integer :: Qrows,Qcols
integer :: i,j,lda,info,ldz
integer :: lwork,liwork,lrwork,mW,mZ
complex(8),dimension(:),allocatable :: Work,RRwork
real(8),dimension(:),allocatable :: Rwork
integer,dimension(:),allocatable :: Iwork
complex(8) :: guess_lwork(1)
complex(8) :: guess_lrrwork(1)
real(8) :: guess_lrwork(1)
integer :: guess_liwork(1)
logical :: boolV,boolI
real(8),dimension(:,:),allocatable :: Z
integer,dimension(:),allocatable :: Isuppz
integer,dimension(:),allocatable :: Ifail
integer,dimension(:),allocatable :: Iclustr
real(8),dimension(:),allocatable :: Gap
!
integer,external :: numroc,indxG2L,indxL2G
real(8),external :: dlamch
!
complex(8),dimension(:,:),allocatable :: A_loc,Z_loc
integer :: rankX,rankY
integer :: Nrow,Ncol
integer :: sendR,sendC
integer :: myi,myj,unit
integer,dimension(9) :: descA,descAloc,descZloc
real(8) :: t_stop,t_start
logical :: master
!
method_='zheevr';if(present(method))method_=trim(method)
!
jobz_='V' ;if(present(jobz))jobz_=jobz
uplo_='L' ;if(present(uplo))uplo_=uplo
vl_ = 1d0 ;if(present(vL))vL_=vL
vu_ = 1d0 ;if(present(vU))vU_=vU
iL_ = 1 ;if(present(iL))iL_=iL
iU_ = 1 ;if(present(iU))iU_=iU
tol_ = dlamch('s') ;if(present(tol))tol_=tol
!
!>DEBUG
if(method_=="zheev" .OR. method_=="zheevd") then
write(0,"(A)")"WARNING p_eigh: pzheev & pzheevd are bugged. Fall on pzheevr."
method_='zheevr'
endif
if(uplo_=='U') then
write(0,"(A)")"WARNING p_eigh: pzheevr OR pzheevx are bugged. Only work with with uplo=L"
uplo_ ='L'
endif
!<DEBUG
!
W=0d0
!
range='A'
boolV=present(vL).AND.present(vU)
boolI=present(iL).OR.present(iU)
if(boolV.and.boolI)stop "vX and iX arguments both present. Can not set RANGE"
if(boolV)range='V'
if(boolI)range='I'
Ns = max(1,size(A,1))
if(any(shape(A)/=[Ns,Ns]))stop "my_eighD error: A has illegal shape"
!
!INIT SCALAPACK TREATMENT:
!< Get coordinate of the processes
call blacs_gridinfo( p_context, p_Nx, p_Ny, rankX, rankY)
master = (rankX==0).AND.(rankY==0)
if(rankX<0.AND.rankY<0)goto 200
!
Nb = Nblock
!
Qrows = numroc(Ns, Nb, rankX, 0, p_Nx)
Qcols = numroc(Ns, Nb, rankY, 0, p_Ny)
!
call descinit( descA, Ns, Ns, Nb, Nb, 0, 0, p_context, Qrows, info )
call descinit( descAloc, Ns, Ns, Nb, Nb, 0, 0, p_context, Qrows, info )
!
if(master)then
unit = free_unit()
open(unit,file="p_eigh.info")
write(unit,"(A20,I8,A5,I8)")"Grid=",p_Nx,"x",p_Ny
write(unit,"(A20,I8,A5,I8)")"Qrows x Qcols=",Qrows,"x",Qcols
endif
!
!< allocate local distributed A
allocate(A_loc(Qrows,Qcols));A_loc=zero
!
!< Distribute A
call Distribute_BLACS(A,A_loc,descAloc,unit)
!
!< Allocate distributed eigenvector matrix
allocate(Z_loc(Qrows,Qcols));Z_loc=zero
call descinit( descZloc, Ns, Ns, Nb, Nb, 0, 0, p_context, Qrows, info )
!
if(master)write(unit,"(A20,A21)")"Using Method:",method_
if(master)call cpu_time(t_start)
select case(method_)
case default
call PZHEEVR(jobz_,range,uplo_,&
Ns,A_loc,1,1,descAloc,vl_,vu_,il_,iu_,mW,mZ,W,Z_loc,1,1,descZloc,&
guess_lwork,-1,guess_lrwork,-1,guess_liwork,-1,info)
lwork = guess_lwork(1) ; lrwork= guess_lrwork(1) ; liwork= guess_liwork(1)
allocate(work(lwork))
allocate(rwork(lrwork))
allocate(iwork(liwork))
call PZHEEVR(jobz_,range,uplo_,&
Ns,A_loc,1,1,descAloc,vl_,vu_,il_,iu_,mW,mZ,W,Z_loc,1,1,descZloc,&
work, lwork, rwork, lrwork, iwork, liwork,info)
!
case ("zheevx")
allocate(Ifail(Ns))
allocate(Iclustr(2*p_Nx*p_Ny))
allocate(Gap(p_Nx*p_Ny))
call PZHEEVX(jobz_,range,uplo_,&
Ns,A_loc,1,1,descAloc,vl_,vu_,il_,iu_,0d0,mW, mZ,W,-1d0,Z_loc,1,1,descAloc,&
guess_lwork,-1,guess_lrwork,-1,guess_liwork,-1,ifail,iclustr,gap,info)
lwork = guess_lwork(1)
lrwork= guess_lrwork(1)
liwork= guess_liwork(1)
allocate(work(lwork))
allocate(rwork(lrwork))
allocate(iwork(liwork))
call PZHEEVX(jobz_,range,uplo_,&
Ns,A_loc,1,1,descAloc,vl_,vu_,il_,iu_,0d0,mW,mZ,W,-1d0,Z_loc,1,1,descZloc,&
work,lwork,rwork,lrwork,iwork,liwork,ifail,iclustr,gap,info)
!
!
! >>>ACTHUNG<< BUGGED DO NOT USE: --> zheevr, w UPLO='L'
case ("zheev") !Bugged
call PZHEEV(jobz_,uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
guess_lwork,-1,guess_lrwork,-1,info)
lwork = guess_lwork(1) ; lrwork = guess_lrwork(1)
allocate(work(lwork))
allocate(rwork(lrwork))
call PZHEEV(jobz_,uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
work,lwork,rwork,lrwork,info)
!
!
! >>>ACTHUNG<< BUGGED DO NOT USE: --> zheevr, w UPLO='L'
case ("zheevd")
call PZHEEVD('V',uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
guess_lwork,-1,guess_lrwork,-1,guess_liwork,-1,INFO)
lwork = guess_lwork(1) ; lrwork= guess_lrwork(1) ; liwork= guess_liwork(1)
allocate(work(lwork))
allocate(rwork(lrwork))
allocate(iwork(liwork))
call PZHEEVD('V',uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
work,lwork,rwork,lrwork,iwork,liwork,info)
end select
if(master)call cpu_time(t_stop)
if(master)write(unit,"(A20,F21.12)")"Time diag:",t_stop-t_start
!
if(jobz_=='V')then
A=zero
call Gather_BLACS(Z_loc,A,descA,unit)
endif
!
if(master)close(unit)
200 continue
return
end subroutine p_zeigh_simple