Add indxg2p and expose numroc

Extend test_gemr2d to use these to describe the resulting matrices and
their storage locations.

lib/blacsfx.fpp

@@ -642,5 +642,4 @@ contains
   end subroutine blacsfx_exit
 
 
-
 end module blacsfx_module

lib/scalapack.fpp

@@ -9,7 +9,7 @@ module scalapack_module
 
   public :: psygst, phegst, psyev, pheev, psyevd, pheevd, psyevr, pheevr
   public :: ptrsm, ppotrf, ppotri, ptrtri, pgetrf, pgesvd
-  public :: sl_init, numroc, infog2l, indxl2g, descinit
+  public :: sl_init, numroc, infog2l, indxl2g, descinit, indxg2p
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!! ppotrf
@@ -516,6 +516,12 @@ module scalapack_module
       integer, intent(in) :: indxglob, nb, iproc, isrcproc, nprocs
     end function indxl2g
 
+    !> Finds processor id where global index is stored.
+    function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
+      integer :: indxg2p
+      integer, intent(in) :: indxglob, nb, iproc, isrcproc, nprocs
+    end function indxg2p
+
     !> Initializes a descriptor for a distributed array.
     subroutine descinit(desc, mm, nn, mb, nb, irsrc, icsrc, ictxt, lld, info)
       integer, intent(out) :: desc(*)

lib/scalapackfx.fpp

@@ -35,6 +35,8 @@ module scalapackfx_module
   public :: scalafx_localindices
   public :: scalafx_creatematrix
   public :: scalafx_pgesvd
+  public :: scalafx_numroc
+  public :: scalafx_indxg2p
 
   !> Cholesky factorization of a symmetric/Hermitian positive definite matrix.
   interface scalafx_ppotrf
@@ -2071,17 +2073,38 @@ contains
   !> Maps local row or column index onto global matrix position.
   !!
   !! \param indxloc  Local index on input.
-  !! \param mygrid  BLACS descriptor.
-  !! \param blocksize  Block size for direction (row or column)
+  !! \param nb  Block size for the column/row
+  !! \param iproc  Local array index
+  !! \param isrcproc coordinate of the process that possesses the first row/column of the
+  !! distributed matrix
+  !! \param nprocs  Total number of processes over which the matrix is distributed
   !!
-  function scalafx_indxl2g(indxloc, crB, mycr, crsrc, ncr)
+  function scalafx_indxl2g(indxloc, nb, iproc, isrcproc, nprocs)
     integer :: scalafx_indxl2g
-    integer, intent(in) :: indxloc, crB, mycr, crsrc, ncr
+    integer, intent(in) :: indxloc, nb, iproc, isrcproc, nprocs
 
-    scalafx_indxl2g = indxl2g(indxloc, crB, mycr, crsrc, ncr)
+    scalafx_indxl2g = indxl2g(indxloc, nb, iproc, isrcproc, nprocs)
 
   end function scalafx_indxl2g
 
+
+  !> Maps global matrix position onto processor.
+  !!
+  !! \param indxloc  Global index on input.
+  !! \param nb  Block size for the column/row
+  !! \param isrcproc coordinate of the process that possesses the first row/column of thexs
+  !! distributed matrix
+  !! \param nprocs  Total number of processes over which the matrix is distributed
+  !!
+  function scalafx_indxg2p(indxglob, nb, isrcproc, nprocs)
+    integer :: scalafx_indxg2p
+    integer, intent(in) :: indxglob, nb, isrcproc, nprocs
+    integer :: iDummy
+
+    scalafx_indxg2p = indxg2p(indxglob, nb, iDummy, isrcproc, nprocs)
+
+  end function scalafx_indxg2p
+
   !> Maps a global position in a distributed matrix to local one.
   !!
   subroutine scalafx_localindices(mygrid, desc, grow, gcol, local, lrow, lcol)
@@ -2122,4 +2145,18 @@ contains
   end subroutine scalafx_localindices
 
 
+  !> Number of rows or columns of distributed matrix owned by a specific process.
+  !!
+  !! \param nn Number of columns (or rows) of global matrix.
+  !! \param nb Column (or row) block size.
+  !! \param iproc The coordinate of the process whose local array row or column is to be determined.
+  !! \param nproc The total number processes over which the matrix is distributed.
+  function scalafx_numroc(nn, nb, iproc, isrcproc, nproc)
+    integer :: scalafx_numroc
+    integer, intent(in) :: nn, nb, iproc, isrcproc, nproc
+
+    scalafx_numroc = numroc(nn, nb, iproc, isrcproc, nproc)
+
+  end function scalafx_numroc
+
 end module scalapackfx_module

test/test_gemr2d.f90

@@ -14,18 +14,18 @@ program test_gemr2d
 
 contains
 
-  ! number elements in global matrix
+  ! number up the elements in the global matrix
   pure function numelem(iGlob, jGlob, mm)
     integer, intent(in) :: iGlob, jGlob, mm
     integer :: numelem
     numelem = jGlob + (iGlob-1) * mm
   end function numelem
 
   subroutine main()
-    type(blacsgrid) :: mygridA, mygridB
+    type(blacsgrid) :: myGridA, myGridB
     integer, allocatable :: AA(:,:), BB(:,:)
     integer :: descA(DLEN_), descB(DLEN_)
-    integer :: nprow, npcol, mm, nn, iproc, nproc, ii, jj, iGlob, jGlob
+    integer :: nprow, npcol, mm, nn, iproc, nproc, ii, jj, iGlob, jGlob, iprow, ipcol
     integer, allocatable :: iErr(:)
 
     ! Block size for matrix B
@@ -42,94 +42,161 @@ subroutine main()
       end if
     end do
     nprow = nproc / npcol
-    call mygridA%initgrid(nprow, npcol)
-    call mygridB%initgrid(1, nproc)
-    if (mygridA%lead) then
+    call myGridA%initgrid(nprow, npcol)
+    call myGridB%initgrid(1, nproc)
+    if (myGridA%lead) then
       write(stdout, "(A,2(1X,I0))") "# source processor grid:", nprow, npcol
     end if
-    if (mygridB%lead) then
+    if (myGridB%lead) then
       write(stdout, "(A,2(1X,I0))") "# destination processor grid:", 1, nproc
     end if
 
-    if (mygridA%lead) then
+    call blacsfx_barrier(myGridA)
+    call blacsfx_barrier(myGridB)
+
+    if (myGridA%lead) then
       write(stdOut, *)'# Matrix size to re-distrbute?'
       read(stdin, *) mm, nn
-      call blacsfx_gebs(mygridA, mm)
-      call blacsfx_gebs(mygridA, nn)
+      call blacsfx_gebs(myGridA, mm)
+      call blacsfx_gebs(myGridA, nn)
     else
-      call blacsfx_gebr(mygridA, mm)
-      call blacsfx_gebr(mygridA, nn)
+      call blacsfx_gebr(myGridA, mm)
+      call blacsfx_gebr(myGridA, nn)
     end if
 
     ! Block size for matrix B
     mbB = mm
     nbB = 1
 
-    if (mygridA%lead) then
+    if (myGridA%lead) then
       write(stdout, "(A,2(1X,I0))") "# A processor grid:", nprow, npcol
       write(stdout, "(A,1X,I0,1X,I0)") "# A block size:", mbA, nbA
     end if
-    if (mygridB%lead) then
+    if (myGridB%lead) then
       write(stdout, "(A,2(1X,I0))") "# B processor grid:", nprow, npcol
       write(stdout, "(A,1X,I0,1X,I0)") "# B block size:", mbB, nbB
     end if
 
-    call scalafx_creatematrix(mygridA, mm, nn, mbA, nbA, AA, descA)
-    call scalafx_creatematrix(mygridB, mm, nn, mbB, nbB, BB, descB)
+    call scalafx_creatematrix(myGridA, mm, nn, mbA, nbA, AA, descA)
+    call scalafx_creatematrix(myGridB, mm, nn, mbB, nbB, BB, descB)
+
+    if (myGridA%lead) then
+      do ii = 0, nproc-1
+        call blacsfx_pcoord(myGridA, ii, iprow, ipcol)
+        write(stdout, "(A,I0,A,1X,I0,1X,I0)") "# A block (proc:",ii,"):",&
+            & max(1, scalafx_numroc(descA(M_), descA(MB_), iprow, descA(RSRC_), myGridA%nrow)),&
+            & scalafx_numroc(descA(N_), descA(NB_), ipcol, descA(CSRC_), myGridA%ncol)
+      end do
+
+      write(stdout,*)'Processor locations in the A BLACS grid'
+      do ii = 0, nproc - 1
+        call blacsfx_pcoord(myGridA, ii, iprow, ipcol)
+        write(stdout, "(X,A,I0,A,I0,X,I0)")'Proc ', ii, ' is at BLACS grid loc. ', iprow, ipcol
+      end do
+
+      write(stdout,*)'Matrix element locations in the A BLACS grid'
+      do ii = 1, nn
+        do jj = 1, mm
+          iprow = scalafx_indxg2p(jj, descA(MB_), descA(RSRC_), myGridA%nrow)
+          ipcol = scalafx_indxg2p(ii, descA(NB_), descA(CSRC_), myGridA%ncol)
+          write(stdout, "(I0,X,I0,A,I0,X,I0,A,I0)")jj, ii, ' is at BLACS location ', ipcol, iprow,&
+              & ' proc:', blacsfx_pnum(myGridA, iprow, ipcol)
+        end do
+      end do
+    end if
+
+    call blacsfx_barrier(myGridA)
+    call blacsfx_barrier(myGridB)
+
+    if (myGridB%lead) then
+      do ii = 0, nproc - 1
+        call blacsfx_pcoord(myGridB, ii, iprow, ipcol)
+        write(stdout, "(A,I0,A,1X,I0,1X,I0)") "# B block (proc:",ii,"):",&
+            & max(1, scalafx_numroc(descB(M_), descB(MB_), iprow, descB(RSRC_), myGridB%nrow)),&
+            & scalafx_numroc(descB(N_), descB(NB_), ipcol, descB(CSRC_), myGridB%ncol)
+      end do
+
+      write(stdout,*)'Processor locations in the B BLACS grid'
+      do ii = 0, nproc - 1
+        call blacsfx_pcoord(myGridB, ii, iprow, ipcol)
+        write(stdout, "(X,A,I0,A,I0,X,I0)")'Proc ', ii, ' is at BLACS grid loc. ', iprow, ipcol
+      end do
+
+      write(stdout,*)'Matrix element locations in the B BLACS grid'
+      do ii = 1, nn
+        do jj = 1, mm
+          iprow = scalafx_indxg2p(jj, descB(MB_), descB(RSRC_), myGridB%nrow)
+          ipcol = scalafx_indxg2p(ii, descB(NB_), descB(CSRC_), myGridB%ncol)
+          write(stdout, "(I0,X,I0,A,I0,X,I0,A,I0)")jj, ii, ' is at BLACS location ', ipcol, iprow,&
+              & ' proc:', blacsfx_pnum(myGridB, iprow, ipcol)
+        end do
+      end do
+
+    end if
+
+    call blacsfx_barrier(myGridA)
+    call blacsfx_barrier(myGridB)
 
     AA(:,:) = 0
     BB(:,:) = 0
 
-    write(stdOut,"(A,I0,A,I0,',',I0,A)")'Processor ', mygridA%iproc,' holds A(', shape(AA),&
-        & ') elements'
+    write(stdOut,"(A,I0,A,I0,',',I0,A,I0,X,I0)")'Processor ', myGridA%iproc,' (grid A) holds A(',&
+        & shape(AA), ') elements vs ',&
+        & max(1, scalafx_numroc(descA(M_), descA(MB_), myGridA%myrow, descA(RSRC_), myGridA%nrow)),&
+        & scalafx_numroc(descA(N_), descA(NB_), myGridA%mycol, descA(CSRC_), myGridA%ncol)
 
-    write(stdOut,"(A,I0,A,I0,',',I0,A)")'Processor ', mygridB%iproc,' holds B(', shape(BB),&
-        & ') elements'
+    write(stdOut,"(A,I0,A,I0,',',I0,A,I0,X,I0)")'Processor ', myGridB%iproc,' (grid B) holds B(',&
+        & shape(BB), ') elements vs ',&
+        & max(1, scalafx_numroc(descB(M_), descB(MB_), myGridB%myrow, descB(RSRC_), myGridB%nrow)),&
+        & scalafx_numroc(descB(N_), descB(NB_), myGridB%mycol, descB(CSRC_), myGridB%ncol)
 
-    ! number off elements in A
+    ! number off the elements in A
     do ii = 1, size(AA, dim=2)
-      iGlob = scalafx_indxl2g(ii, descA(NB_), mygridA%mycol, descA(CSRC_), mygridA%ncol)
+      iGlob = scalafx_indxl2g(ii, descA(NB_), myGridA%mycol, descA(CSRC_), myGridA%ncol)
       do jj = 1, size(AA, dim=1)
-        jGlob = scalafx_indxl2g(jj, descA(MB_), mygridA%myrow, descA(RSRC_), mygridA%nrow)
+        jGlob = scalafx_indxl2g(jj, descA(MB_), myGridA%myrow, descA(RSRC_), myGridA%nrow)
         AA(jj,ii) = numelem(iGlob, jGlob, mm)
       end do
     end do
 
+    call blacsfx_barrier(myGridA)
+    call blacsfx_barrier(myGridB)
+
     ! A -> B
-    call blacsfx_gemr2d(mm, nn, aa, 1, 1, descA, bb, 1, 1, descB, mygridA%ctxt)
+    call blacsfx_gemr2d(mm, nn, aa, 1, 1, descA, bb, 1, 1, descB, myGridA%ctxt)
 
     aa(:,:) = 0
 
     ! B -> A
-    call blacsfx_gemr2d(mm, nn, bb, 1, 1, descB, aa, 1, 1, descA, mygridA%ctxt)
+    call blacsfx_gemr2d(mm, nn, bb, 1, 1, descB, aa, 1, 1, descA, myGridA%ctxt)
 
     allocate(iErr(0:nProc-1))
     iErr(:) = 0
     ! check it's what is expected
     if (all(shape(AA) > [0,0])) then
       lpCol: do ii = 1, size(AA, dim=2)
-        iGlob = scalafx_indxl2g(ii, descA(NB_), mygridA%mycol, descA(CSRC_), mygridA%ncol)
+        iGlob = scalafx_indxl2g(ii, descA(NB_), myGridA%mycol, descA(CSRC_), myGridA%ncol)
         if (iGlob > nn) then
           cycle lpCol
         end if
         lpRow: do jj = 1, size(AA, dim=1)
-          jGlob = scalafx_indxl2g(jj, descA(MB_), mygridA%myrow, descA(RSRC_), mygridA%nrow)
+          jGlob = scalafx_indxl2g(jj, descA(MB_), myGridA%myrow, descA(RSRC_), myGridA%nrow)
           if (jGlob > mm) then
             cycle lpRow
           end if
           if (AA(jj,ii) /= numelem(iGlob, jGlob, mm)) then
-            write(stdOut,*) "Error on processor ", mygridA%iproc, " Mismatch in element ", jGlob,&
+            write(stdOut,*) "Error on processor ", myGridA%iproc, " Mismatch in element ", jGlob,&
                 & iGlob
-            iErr(mygridA%iproc) = numelem(iGlob, jGlob, mm)
+            iErr(myGridA%iproc) = numelem(iGlob, jGlob, mm)
             exit lpcol
           end if
         end do lpRow
       end do lpCol
     end if
 
-    call blacsfx_gsum(mygridA, iErr)
+    call blacsfx_gsum(myGridA, iErr)
 
-    if (mygridB%lead) then
+    if (myGridB%lead) then
       if (any(iErr /= 0)) then
         write(stdOut,*)"Errors for matrix elements"
         write(stdOut,*) iErr