p_Dinv Subroutine

subroutine p_Dinv(A,Nblock)
  real(8),dimension(:,:),intent(inout)       :: A
  integer                                    :: Nblock
  integer                                    :: Nb
  integer                                    :: Ns
  integer                                    :: Qrows,Qcols
  integer                                    :: i,j,lda,info
  !
  real(8)                                    :: guess_lwork(1)
  integer                                    :: guess_liwork(1)
  real(8),dimension(:),allocatable           :: work
  integer,dimension(:),allocatable           :: Iwork 
  integer                                    :: lwork
  integer                                    :: liwork
  integer, allocatable                       :: ipiv(:)
  !
  integer,external                           :: numroc,indxG2L,indxG2P,indxL2G
  real(8),external                           :: dlamch,Pilaenvx
  !
  real(8),dimension(:,:),allocatable         :: A_loc,Z_loc
  integer                                    :: rankX,rankY
  integer                                    :: sendR,sendC,Nipiv
  integer                                    :: Nrow,Ncol
  integer                                    :: myi,myj,unit
  integer,dimension(9)                       :: descA,descAloc,descZloc
  real(8)                                    :: t_stop,t_start
  logical                                    :: master
  !
  !
  Ns    = max(1,size(A,1))
  if(any(shape(A)/=[Ns,Ns]))stop "my_eighD error: A has illegal shape"
  !
  !< 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_inv.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)
  !
  if(master)call cpu_time(t_start)
  allocate(Ipiv(Qrows*Nb))
  call PDGETRF(Ns, Ns, A_loc, 1, 1, descAloc, IPIV, INFO)
  if(info /= 0) then
     print*, "PDGETRF returned info =", info
     stop
  end if
  !
  !
  call PDGETRI( Ns, A_loc, 1, 1, descAloc, IPIV, guess_lWORK, -1, guess_LIWORK, -1, INFO )
  lwork = guess_lwork(1)
  liwork= guess_liwork(1)
  allocate(work(lwork))
  allocate(iwork(liwork))
  call PDGETRI( Ns, A_loc, 1, 1, descAloc, IPIV, Work, Lwork, Iwork, LIwork, INFO )
  if(info /= 0) then
     print*, "PDGETRI ERROR. returned info =", info
     stop
  end if
  !
  if(master)call cpu_time(t_stop)
  if(master)write(unit,"(A20,F21.12)")"Time Invert :",t_stop-t_start
  !
  A=0d0
  call Gather_BLACS(A_loc,A,descA,unit)
  !
  if(master)close(unit)
  !
100 continue
  return
  !
end subroutine p_Dinv