diff --git a/lib/blacsgrid.fpp b/lib/blacsgrid.fpp
index 5b48e8f..4a99429 100644
--- a/lib/blacsgrid.fpp
+++ b/lib/blacsgrid.fpp
@@ -18,6 +18,7 @@ module blacsgrid_module
     integer :: nproc = -1  !< Nr. of processes in the grid.
     integer :: myrow = -1  !< Row of the current process.
     integer :: mycol = -1  !< Column of the current process.
+    integer :: leadproc = -1  !< Id of the lead process.
     integer :: leadrow = -1   !< Row of the lead process.
     integer :: leadcol = -1   !< Column of the lead process.
     logical :: lead = .false. !< Whether the current process is the lead.
@@ -137,8 +138,8 @@ contains
       end if
       self%leadrow = leadrow0
       self%leadcol = leadcol0
-      self%lead = (self%myrow == self%leadrow &
-          & .and. self%mycol == self%leadcol)
+      self%leadproc = blacs_pnum(self%ctxt, self%leadrow, self%leadcol)
+      self%lead = (self%myrow == self%leadrow .and. self%mycol == self%leadcol)
     end if
 
   end subroutine initgrid
diff --git a/lib/scalapackfx.fpp b/lib/scalapackfx.fpp
index e826304..66cf524 100644
--- a/lib/scalapackfx.fpp
+++ b/lib/scalapackfx.fpp
@@ -30,6 +30,7 @@ module scalapackfx_module
   public :: scalafx_ptrsm
   public :: scalafx_getdescriptor
   public :: scalafx_getlocalshape
+  public :: scalafx_getremoteshape
   public :: scalafx_infog2l
   public :: scalafx_indxl2g
   public :: scalafx_localindices
@@ -37,6 +38,7 @@ module scalapackfx_module
   public :: scalafx_pgesvd
   public :: scalafx_numroc
   public :: scalafx_indxg2p
+  public :: scalafx_infog2p
 
   !> Cholesky factorization of a symmetric/Hermitian positive definite matrix.
   interface scalafx_ppotrf
@@ -1961,6 +1963,45 @@ contains
 
   end subroutine scalafx_getlocalshape
 
+
+  !> Returns the shape of a remote part of a distributed array, given the processor number
+  !!
+  !! \param mygrid  BLACS grid descriptor.
+  !! \param desc  Global matrix descriptor.
+  !! \param iproc  Processor number
+  !! \param nrowloc  Nr. of local rows, returns 0 if iproc outside of grid.
+  !! \param ncolloc  Nr. of local columns, returns 0 if iproc outside of grid.
+  !!
+  subroutine scalafx_getremoteshape(mygrid, desc, iproc, nrowloc, ncolloc, storage)
+    type(blacsgrid), intent(in) :: mygrid
+    integer, intent(in) :: desc(DLEN_), iproc
+    integer, intent(out) :: nrowloc, ncolloc
+    logical, intent(in), optional :: storage
+
+    integer :: iprow, ipcol
+    logical :: isStorage
+
+    isStorage = .false.
+    if (present(storage)) then
+      isStorage = storage
+    end if
+
+    if (iproc < 0 .or. iproc >= mygrid%nproc) then
+      nrowloc = 0
+      ncolloc = 0
+    else
+      call blacsfx_pcoord(myGrid, iproc, iprow, ipcol)
+      if (isStorage) then
+        nrowloc = max(1,numroc(desc(M_), desc(MB_), iprow, desc(RSRC_), mygrid%nrow))
+      else
+        nrowloc = numroc(desc(M_), desc(MB_), iprow, desc(RSRC_), mygrid%nrow)
+      end if
+      ncolloc = numroc(desc(N_), desc(NB_), ipcol, desc(CSRC_), mygrid%ncol)
+    end if
+
+  end subroutine scalafx_getremoteshape
+
+
   !> Maps global position in a distributed matrix to local one.
   !!
   !! \param mygrid  BLACS descriptor.
@@ -2105,6 +2146,40 @@ contains
 
   end function scalafx_indxg2p
 
+
+  !> Processor grid location holding a global matrix element. If the element is outside of the
+  !! matrix, returns location (-1,-1)
+  subroutine scalafx_infog2p(mygrid, desc, grow, gcol, prow, pcol)
+
+    !> BLACS descriptor.
+    type(blacsgrid), intent(in) :: mygrid
+
+    !> Descriptor of the distributed matrix.
+    integer, intent(in) :: desc(DLEN_)
+
+    !> Global row index.
+    integer, intent(in) :: grow
+
+    !> Global column index
+    integer, intent(in) :: gcol
+
+    !> Row index in the BLACS grid
+    integer, intent(out) :: prow
+
+    !> Column index in BLACS grid
+    integer, intent(out) :: pcol
+
+    if (grow < 0 .or. grow > desc(M_) .or. gcol < 0 .or. gcol > desc(N_)) then
+      prow = -1
+      pcol = -1
+    else
+      prow = scalafx_indxg2p(grow, desc(MB_), desc(RSRC_), myGrid%nrow)
+      pcol = scalafx_indxg2p(gcol, desc(NB_), desc(CSRC_), myGrid%ncol)
+    end if
+
+  end subroutine scalafx_infog2p
+
+
   !> Maps a global position in a distributed matrix to local one.
   !!
   subroutine scalafx_localindices(mygrid, desc, grow, gcol, local, lrow, lcol)
@@ -2121,25 +2196,29 @@ contains
     !> Global column index
     integer, intent(in) :: gcol
 
-    !> Indicates whether given global index is local for the process.
+    !> Indicates whether given global index is local for the calling process.
     logical, intent(out) :: local
 
-    !> Row index in the local matrix (or 0 if global index is not local)
+    !> Row index in the local matrix on the process holding that matrix, if outside the matrix
+    !> returns 0
     integer, intent(out) :: lrow
 
-    !> Column index in the local matrix (or 0 if global index is not local)
+    !> Column index in the local matrix on the process holding that matrix, if outside the matrix
+    !> returns 0
     integer, intent(out) :: lcol
 
     !------------------------------------------------------------------------
 
     integer :: rsrc, csrc
 
-    call infog2l(grow, gcol, desc, mygrid%nrow, mygrid%ncol, mygrid%myrow,&
-        & mygrid%mycol, lrow, lcol, rsrc, csrc)
-    local = (rsrc == mygrid%myrow .and. csrc == mygrid%mycol)
-    if (.not. local) then
+    if (grow < 0 .or. grow > desc(M_) .or. gcol < 0 .or. gcol > desc(N_)) then
       lrow = 0
       lcol = 0
+      local = .false.
+    else
+      call infog2l(grow, gcol, desc, mygrid%nrow, mygrid%ncol, mygrid%myrow,&
+          & mygrid%mycol, lrow, lcol, rsrc, csrc)
+      local = (rsrc == mygrid%myrow .and. csrc == mygrid%mycol)
     end if
 
   end subroutine scalafx_localindices
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 38d204d..0f14ee7 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -9,6 +9,7 @@ set(common-dep-targets
   test_gemr2d
   test_linecomm
   test_psyr_pher
+  test_remoteelements
   test_svd
   test_tran)
 
diff --git a/test/test_remoteelements.f90 b/test/test_remoteelements.f90
new file mode 100644
index 0000000..7b1c79e
--- /dev/null
+++ b/test/test_remoteelements.f90
@@ -0,0 +1,145 @@
+!> Testing copy between matrix patterns
+program test_remoteelements
+  use, intrinsic :: iso_fortran_env, stdout => output_unit, stdin => input_unit
+  use test_common_module
+  use libscalapackfx_module
+  use blacsfx_module
+  implicit none
+
+  call main()
+
+contains
+
+  subroutine main()
+    type(blacsgrid) :: myGrid
+    integer, allocatable :: AA(:,:), AAreceive(:,:), All(:,:)
+    integer :: desc(DLEN_), nprow, npcol, mm, nn, mb, nb, iproc, nproc, ii, jj, kk, iGlob, jGlob
+    integer :: iprow, ipcol, iLoc, jLoc, locrow, loccol
+    logical :: isLocal
+
+    ! grid information
+    call blacsfx_pinfo(iproc, nproc)
+
+    ! approximately square grid for group
+    do npcol = floor(sqrt(real(nproc, dp))), nproc
+      if (mod(nproc, npcol) == 0) then
+        exit
+      end if
+    end do
+    nprow = nproc / npcol
+    call myGrid%initgrid(nprow, npcol)
+
+    if (myGrid%lead) then
+      write(stdout, "(A,2(1X,I0))") "# source processor grid:", nprow, npcol
+      write(stdOut, "(A)")'# Matrix size to re-distrbute?'
+      read(stdin, *) mm, nn
+      write(stdout,"(A)")'# Block sizes of matrix?'
+      read(stdin,*) mb, nb
+      call blacsfx_gebs(myGrid, mm)
+      call blacsfx_gebs(myGrid, nn)
+      call blacsfx_gebs(myGrid, mb)
+      call blacsfx_gebs(myGrid, nb)
+    else
+      call blacsfx_gebr(myGrid, mm)
+      call blacsfx_gebr(myGrid, nn)
+      call blacsfx_gebr(myGrid, mb)
+      call blacsfx_gebr(myGrid, nb)
+    end if
+
+    call scalafx_creatematrix(myGrid, mm, nn, mb, nb, AA, desc)
+
+    if (myGrid%lead) then
+      do ii = 0, nproc-1
+        call scalafx_getremoteshape(myGrid, desc, ii, locrow, loccol)
+        write(stdout, "(A,I0,A,1X,I0,1X,I0)") "# A block (data:",ii,"):", locrow, loccol
+        call scalafx_getremoteshape(myGrid, desc, ii, locrow, loccol, .true.)
+        write(stdout, "(T8,A,I0,A,1X,I0,1X,I0)") "(storage:",ii,"):", locrow, loccol
+      end do
+
+      do ii = 1, nn
+        do jj = 1, mm
+          call scalafx_infog2p(mygrid, desc, jj, ii, iprow, ipcol)
+          write(stdout, "(A,I0,1X,I0,A,I0,1X,I0,A,I0)") "# Matrix elements (",jj, ii,&
+              & ') stored on BLACS grid at ', iprow, ipcol, " proc. ",&
+              & blacsfx_pnum(myGrid, iprow, ipcol)
+        end do
+      end do
+
+    end if
+
+    call blacsfx_barrier(myGrid)
+
+    AA(:,:) = -1
+    call scalafx_getremoteshape(myGrid, desc, iProc, locrow, loccol, .true.)
+    call blacsfx_pcoord(myGrid, iProc, iprow, ipcol)
+    do ii = 1, loccol
+      iGlob = scalafx_indxl2g(ii, desc(NB_), ipcol, desc(CSRC_), myGrid%ncol)
+      do jj = 1, locrow
+        jGlob = scalafx_indxl2g(jj, desc(MB_), iprow, desc(RSRC_), myGrid%nrow)
+        AA(jj,ii) = jGlob + (iGlob-1) * mm
+      end do
+    end do
+    write(stdout,*)iproc,':', AA(:locrow,:loccol)
+
+    call blacsfx_barrier(myGrid)
+
+    if (myGrid%lead) then
+
+      allocate(All(mm,nn), source = 0)
+
+      do kk = 0, nproc-1
+
+        if (kk /=  myGrid%leadproc) then
+          write(stdout,*)'Follow', kk
+
+          call scalafx_getremoteshape(myGrid, desc, kk, locrow, loccol, .true.)
+          allocate(AAreceive(locrow, loccol), source=-1)
+
+          call blacsfx_pcoord(myGrid, kk, iprow, ipcol)
+          write(*,*)kk, iprow, ipcol,'Expecting', shape(AAreceive)
+          call blacsfx_gerv(myGrid, AAreceive, iprow, ipcol)
+
+        else
+
+          write(stdout,*)'Lead', kk
+          call scalafx_getlocalshape(mygrid, desc, locrow, loccol)
+          AAreceive = AA(:locrow,:loccol)
+
+        end if
+
+        call scalafx_getremoteshape(myGrid, desc, kk, locrow, loccol, .true.)
+        call blacsfx_pcoord(myGrid, kk, iprow, ipcol)
+        do ii = 1, loccol
+          iGlob = scalafx_indxl2g(ii, desc(NB_), ipcol, desc(CSRC_), myGrid%ncol)
+          do jj = 1, locrow
+            jGlob = scalafx_indxl2g(jj, desc(MB_), iprow, desc(RSRC_), myGrid%nrow)
+            write(stdout,*)'Proc',kk,',', iprow, ipcol,': (', jGlob,iGlob,')(',jj,ii,')'
+            All(jGlob,iGlob) = AAreceive(jj,ii)
+          end do
+        end do
+
+        deallocate(AAreceive)
+
+      end do
+
+    else
+
+      write(*,*)iproc, myGrid%myrow, myGrid%mycol,'sending', shape(AA)
+      call blacsfx_gesd(myGrid, AA, myGrid%leadrow, myGrid%leadcol)
+
+    end if
+
+    call blacsfx_barrier(myGrid)
+
+    if (myGrid%lead) then
+      do jj = 1, nn
+        write(stdout,*)All(:,jj)
+      end do
+    end if
+
+    ! Finish blacs.
+    call blacsfx_exit()
+
+  end subroutine main
+
+end program test_remoteelements