subroutine p_Zinv(A,Nblock)
complex(8),dimension(:,:),intent(inout) :: A
integer :: Nblock
integer :: Nb
integer :: Ns
integer :: Qrows,Qcols
integer :: i,j,lda,info
!
complex(8) :: guess_lwork(1)
integer :: guess_liwork(1)
complex(8),dimension(:),allocatable :: work
integer,dimension(:),allocatable :: Iwork
integer :: lwork
integer :: liwork
integer, allocatable :: ipiv(:)
!
integer,external :: numroc,indxG2L,indxG2P,indxL2G
!
complex(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"
!
!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)
!
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)
!
!< 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)call cpu_time(t_start)
allocate(Ipiv(Qrows*Nb))
call PZGETRF(Ns, Ns, A_loc, 1, 1, descAloc, IPIV, INFO)
if(info /= 0) then
print*, "PZGETRF returned info =", info
stop
end if
!
!
call PZGETRI( 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 PZGETRI( Ns, A_loc, 1, 1, descAloc, IPIV, Work, Lwork, Iwork, LIwork, INFO )
if(info /= 0) then
print*, "PZGETRI returned info =", info
stop
end if
!
if(master)call cpu_time(t_stop)
if(master)write(unit,"(A20,F21.12)")"Time inv A_loc:",t_stop-t_start
!
A=0d0
call Gather_BLACS(A_loc,A,descA,unit)
if(master)close(unit)
200 continue
return
!
end subroutine p_Zinv