subroutine p_deigh_simple(A,W,Nblock,method,jobz,uplo,vl,vu,il,iu,tol)
real(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,mW,mZ
real(8),dimension(:),allocatable :: work,Rwork,Gap
integer :: guess_liwork(1)
real(8) :: guess_lwork(1)
logical :: boolV,boolI
real(8),dimension(:,:),allocatable :: Z
integer,dimension(:),allocatable :: Isuppz,Ifail,Iclustr
integer,dimension(:),allocatable :: Iwork
!
integer,external :: numroc,indxG2L,indxG2P,indxL2G
real(8),external :: dlamch
!
real(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_='dsyevr';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
!
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 100
!
Nb = Nblock
!
Qrows = numroc(Ns, Nb, rankX, 0, p_Nx)
Qcols = numroc(Ns, Nb, rankY, 0, p_Ny)
!
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))
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 )
!
!< Distribute A
call Distribute_BLACS(A,A_loc,descAloc,unit)
!
!< Allocate distributed eigenvector matrix
allocate(Z_loc(Qrows,Qcols));Z_loc=0d0
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 PDSYEVR(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_liwork,-1,info)
lwork = guess_lwork(1)
liwork= guess_liwork(1)
allocate(work(lwork))
allocate(iwork(liwork))
call PDSYEVR(jobz_,range,uplo_,&
Ns,A_loc,1,1,descAloc,vl_,vu_,il_,iu_,mW,mZ,W,Z_loc,1,1,descZloc,&
work, lwork, iwork, liwork,info)
!
case ("dsyev")
call PDSYEV(jobz_,uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
guess_lwork,-1,info)
lwork = guess_lwork(1)
allocate(work(lwork))
call PDSYEV( jobz_,uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
work,lwork,info)
!
!
case ("dsyevd")
call PDSYEVD(jobz_,uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
guess_lwork,-1,guess_liwork,-1,INFO)
lwork = guess_lwork(1)
liwork= guess_liwork(1)
allocate(work(lwork))
allocate(iwork(liwork))
call PDSYEVD(jobz_,uplo_,&
Ns,A_loc,1,1,descAloc,W,Z_loc,1,1,descZloc,&
work,lwork,iwork,liwork,info)
!
case ("dsyevx")
allocate(Ifail(Ns))
allocate(Iclustr(2*p_Nx*p_Ny))
allocate(Gap(p_Nx*p_Ny))
call PDSYEVX(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_liwork,-1,ifail,iclustr,gap,info)
lwork = guess_lwork(1)
liwork= guess_liwork(1)
allocate(work(lwork))
allocate(iwork(liwork))
call PDSYEVX(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,iwork,liwork,ifail,iclustr,gap,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=0d0
call Gather_BLACS(Z_loc,A,descA,unit)
endif
if(master)close(unit)
100 continue
return
!
end subroutine p_deigh_simple