Merge branch 'stage-1.2' into release

CHANGELOG.rst

@@ -4,6 +4,15 @@ Change Log
 
 Notable project changes in various releases.
 
+1.2
+===
+
+Added
+-----
+
+* infog2l accepts also an array of indices
+
+
 1.1
 ===
 

CMakeLists.txt

@@ -8,7 +8,7 @@ include(ScalapackFxUtils)
 
 include(${CMAKE_CURRENT_SOURCE_DIR}/config.cmake)
 
-project(ScalapackFx VERSION 1.1.0 LANGUAGES Fortran)
+project(ScalapackFx VERSION 1.2.0 LANGUAGES Fortran)
 
 setup_build_type()
 

doc/doxygen/Doxyfile

@@ -32,7 +32,7 @@ PROJECT_NAME           = "ScaLAPACKFX"
 # This could be handy for archiving the generated documentation or
 # if some version control system is used.
 
-PROJECT_NUMBER = "1.1.0"
+PROJECT_NUMBER = "1.2.0"
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer

doc/sphinx/conf.py

@@ -45,10 +45,10 @@
 # built documents.
 #
 # The short X.Y version.
-version = '1.1'
+version = '1.2'
 
 # The full version, including alpha/beta/rc tags.
-release = '1.1.0'
+release = '1.2.0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

lib/scalapackfx.fpp

@@ -153,6 +153,11 @@ module scalapackfx_module
     module procedure scalafx_creatematrix_complex, scalafx_creatematrix_dcomplex
   end interface scalafx_creatematrix
 
+  !> Maps global position in a distributed matrix to local one.
+  interface scalafx_infog2l
+    module procedure scalafx_infog2l_single, scalafx_infog2l_array
+  end interface scalafx_infog2l
+
 !************************************************************************
 !*** ppotrf
 !************************************************************************
@@ -1965,7 +1970,8 @@ contains
   !! \param rsrc  Row of the process owning the local matrix.
   !! \param csrc  Column of the process owning the local matrix.
   !!
-  subroutine scalafx_infog2l(mygrid, desc, grow, gcol, lrow, lcol, rsrc, csrc)
+  subroutine scalafx_infog2l_single(mygrid, desc, grow, gcol,&
+      & lrow, lcol, rsrc, csrc)
     type(blacsgrid), intent(in) :: mygrid
     integer, intent(in) :: desc(DLEN_)
     integer, intent(in) :: grow, gcol
@@ -1975,8 +1981,92 @@ contains
     call infog2l(grow, gcol, desc, mygrid%nrow, mygrid%ncol, mygrid%myrow,&
         & mygrid%mycol, lrow, lcol, rsrc, csrc)
 
-  end subroutine scalafx_infog2l
+  end subroutine scalafx_infog2l_single
+
+  !> Maps global positions in a distributed matrix to local one.
+  !!
+  !! \param mygrid  BLACS descriptor.
+  !! \param desc  Descriptor of the distributed matrix.
+  !! \param grow  Global row indices.
+  !! \param gcol  Global column indices.
+  !! \param lrow  Local row indices on output.
+  !! \param lcol  Local column indices on output.
+  !! \param rsrc  Rows of the process owning the local matrix.
+  !! \param csrc  Columns of the process owning the local matrix.
+  !! \param calcAllIndices  Whether to calculate all lrow and lcol,
+  !!     even if the current process does not own them. (default: true)
+  !!
+  subroutine scalafx_infog2l_array(mygrid, desc, grow, gcol,&
+      & lrow, lcol, rsrc, csrc, calcAllIndices)
+    type(blacsgrid), intent(in) :: mygrid
+    integer, intent(in) :: desc(DLEN_)
+    integer, intent(in) :: grow(:), gcol(:)
+    integer, intent(out) :: lrow(:), rsrc(:)
+    integer, intent(out) :: lcol(:), csrc(:)
+    logical, intent(in), optional :: calcAllIndices
+
+    call scalapackfx_infog2l_helper(grow, desc(MB_), desc(RSRC_),&
+        & mygrid%myrow, mygrid%nrow, lrow, rsrc, calcAllIndices)
+
+    call scalapackfx_infog2l_helper(gcol, desc(NB_), desc(CSRC_),&
+        & mygrid%mycol, mygrid%ncol, lcol, csrc, calcAllIndices)
+
+  end subroutine scalafx_infog2l_array
+
+  !> Helper routine for scalafx_infog2l_array.
+  !!
+  !! \param globalInd  Global row/column indices.
+  !! \param descB  Either desc(MB_) or desc(NB_).
+  !! \param descSRC  Either desc(RSRC_) or desc(CSRC_).
+  !! \param myPos  Row/column of the current process.
+  !! \param nPos  Number of rows/columns.
+  !! \param localInd  Local row/column indices on output.
+  !! \param localPos  Rows/columns of the process owning the local matrix.
+  !! \param calcAllIndices  Whether to calculate all local indices,
+  !!     even if the current process does not own them. (default: true)
+  !!
+  subroutine scalapackfx_infog2l_helper(globalInd, descB, descSRC,&
+      & myPos, nPos, localInd, localPos, calcAllIndices)
+    integer, intent(in) :: globalInd(:)
+    integer, intent(in) :: descB, descSRC
+    integer, intent(in) :: myPos, nPos
+    integer, intent(out) :: localInd(:), localPos(:)
+    logical, intent(in), optional :: calcAllIndices
+
+    real(dp) :: inv
+    integer, dimension(size(globalInd)) :: blk
+    integer :: check, i
+    logical :: calcAllIndices_
+
+    ! Note that we explicitly multiply with a double here instead of
+    ! dividing by an integer to enhance performance.
+    inv = 1.0_dp / real(descB, kind=dp)
+    blk = (globalInd - 1) * inv
+
+    check = modulo(myPos - descSRC, nPos)
+
+    localPos = mod(blk + descSRC, nPos)
+
+    calcAllIndices_ = .true.
+    if (present(calcAllIndices)) then
+      calcAllIndices_ = calcAllIndices
+    end if
+
+    do i = 1, size(globalInd)
+      if (calcAllIndices_ .or. myPos == localPos(i)) then
+        localInd(i) = (blk(i) / nPos + 1) * descB + 1
+        if (check >= mod(blk(i), nPos)) then
+          if (myPos == localPos(i)) then
+            localInd(i) = localInd(i) + mod(globalInd(i) - 1, descB)
+          end if
+          localInd(i) = localInd(i) - descB
+        end if
+      else
+        localInd(i) = -1
+      end if
+    end do
 
+  end subroutine scalapackfx_infog2l_helper
 
   !> Maps local row or column index onto global matrix position.
   !!

lib/scalapackfx_tools.fpp

@@ -68,15 +68,33 @@
     integer, intent(in) :: ii, jj
     ${FTYPE}$, intent(inout) :: glob(:,:)
 
-    integer :: i2, j2, iloc, jloc, prow, pcol
+    integer :: i2, j2, nr, nc
+    integer, dimension(size(loc, dim=1)) :: irows, iloc, prow
+    integer, dimension(size(loc, dim=2)) :: icols, jloc, pcol
+
+    nr = size(loc, dim=1)
+    nc = size(loc, dim=2)
+
+    do i2 = 1, nr
+      irows(i2) = i2 + ii - 1
+    end do
+
+    do j2 = 1, nc
+      icols(j2) = j2 + jj - 1
+    end do
+
+    call scalafx_infog2l(mygrid, desc, irows, icols, iloc, jloc,&
+        & prow, pcol, .false.)
 
-    do j2 = 1, size(loc, dim=2)
-      do i2 = 1, size(loc, dim=1)
-        call scalafx_infog2l(mygrid, desc, i2 + ii - 1, j2 + jj - 1, &
-            & iloc, jloc, prow, pcol)
-        if (prow == mygrid%myrow .and. pcol == mygrid%mycol) then
-          glob(iloc, jloc) = glob(iloc, jloc) + loc(i2, j2)
+    do j2 = 1, nc
+      if (pcol(j2) /= mygrid%mycol) then
+        cycle
+      end if
+      do i2 = 1, nr
+        if (prow(i2) /= mygrid%myrow) then
+          cycle
         end if
+        glob(iloc(i2), jloc(j2)) = glob(iloc(i2), jloc(j2)) + loc(i2, j2)
       end do
     end do
 
@@ -109,17 +127,36 @@
     ${FTYPE}$, intent(in) :: glob(:,:)
     ${FTYPE}$, intent(out) :: loc(:,:)
 
-    integer :: i2, j2, iloc, jloc, prow, pcol
+    integer :: i2, j2, nr, nc
+    integer, dimension(size(loc, dim=1)) :: irows, iloc, prow
+    integer, dimension(size(loc, dim=2)) :: icols, jloc, pcol
+
+    nr = size(loc, dim=1)
+    nc = size(loc, dim=2)
+
+    do i2 = 1, nr
+      irows(i2) = i2 + ii - 1
+    end do
+
+    do j2 = 1, nc
+      icols(j2) = j2 + jj - 1
+    end do
+
+    call scalafx_infog2l(mygrid, desc, irows, icols, iloc, jloc,&
+        & prow, pcol, .false.)
 
-    loc(:,:) = 0.0_dp
-    do j2 = 1, size(loc, dim=2)
-      do i2 = 1, size(loc, dim=1)
-        call scalafx_infog2l(mygrid, desc, i2 + ii - 1, j2 + jj - 1, &
-            & iloc, jloc, prow, pcol)
-        if (prow == mygrid%myrow .and. pcol == mygrid%mycol) then
-          loc(i2, j2) = glob(iloc, jloc)
-        end if
-      end do
+    do j2 = 1, nc
+      if (pcol(j2) == mygrid%mycol) then
+        do i2 = 1, nr
+          if (prow(i2) == mygrid%myrow) then
+            loc(i2, j2) = glob(iloc(i2), jloc(j2))
+          else
+            loc(i2, j2) = 0.0_dp
+          end if
+        end do
+      else
+        loc(:, j2) = 0.0_dp
+      end if
     end do
 
   end subroutine cpg2l_${SUFFIX}$

test/CMakeLists.txt

@@ -3,6 +3,7 @@ set(targets
   test_subgrids)
 
 set(common-dep-targets
+  test_cpg2l
   test_det
   test_diag
   test_gemr2d

test/test_cpg2l.f90

@@ -0,0 +1,106 @@
+!> Testing rank one updates.
+program test_cpg2l
+  use, intrinsic :: iso_fortran_env, stdout => output_unit
+  use test_common_module
+  use libscalapackfx_module
+  implicit none
+
+
+  ! Block size (using an extremely small value for test purposes)
+  integer, parameter :: bsize = 2
+
+  call main()
+
+contains
+
+  subroutine main()
+    type(blacsgrid) :: grid1, grid2
+
+    integer :: nprow, npcol, iproc, nproc
+
+    ! Initialize blas and create a square processor grid
+    call blacsfx_pinfo(iproc, nproc)
+    do nprow = int(sqrt(real(nproc, dp))), nproc
+      if (mod(nproc, nprow) == 0) then
+        exit
+      end if
+    end do
+    npcol = nproc / nprow
+
+    call grid1%initgrid(nprow, npcol)
+    if (grid1%lead) then
+      write(stdout, "(A,2(1X,I0))") "# processor grid:", nprow, npcol
+    end if
+
+    call grid2%initgrid(1, nproc)
+    if (grid2%lead) then
+      write(stdout, "(A,2(1X,I0))") "# processor grid:", 1, nproc
+    end if
+
+    if (.not. readMatrixAndTest(grid1, 2, 2, 1, 2)) then
+      write(stdout, "(A)") "Test 1 failed"
+    end if
+    if (.not. readMatrixAndTest(grid1, 5, 5, 1, 1)) then
+      write(stdout, "(A)") "Test 2 failed"
+    end if
+    if (.not. readMatrixAndTest(grid2, 2, 2, 1, 2)) then
+      write(stdout, "(A)") "Test 3 failed"
+    end if
+    if (.not. readMatrixAndTest(grid2, 3, 5, 3, 1)) then
+      write(stdout, "(A)") "Test 4 failed"
+    end if
+
+    call grid1%destruct()
+    call grid1%initgrid(nproc, 1)
+    if (.not. readMatrixAndTest(grid1, 2, 2, 1, 2)) then
+      write(stdout, "(A)") "Test 5 failed"
+    end if
+
+    ! Finish blacs.
+    call blacsfx_exit()
+
+  end subroutine main
+
+  function readMatrixAndTest(mygrid, iSize, jSize, i0, j0) result(success)
+    type(blacsgrid), intent(inout) :: mygrid
+    integer, intent(in) :: iSize, jSize, i0, j0
+    logical :: success
+
+    real(dp), allocatable :: glob(:,:), localTest(:,:), localRef(:,:)
+    integer :: desc(DLEN_)
+    integer :: mm, nn, i, j
+    integer :: iloc, jloc, prow, pcol
+
+    ! Read in matrix from disc.
+    call readfromfile(mygrid, "hamsqr1.dat", bsize, bsize, glob, desc)
+    mm = desc(M_)
+    nn = desc(N_)
+    if (mygrid%lead) then
+      write(stdout, "(A,2(1X,I0))") "# global matrix size:", mm, nn
+      write(stdout, "(A,2(1X,I0))") "# local matrix size on leader:",&
+          & size(glob, dim=1), size(glob, dim=2)
+    end if
+
+    allocate(localRef(iSize,jSize), localTest(iSize,jSize))
+
+    localRef(:,:) = 0.0_dp
+    do j = 1, jSize
+      do i = 1, iSize
+        call scalafx_infog2l(mygrid, desc, i + i0 - 1, j + j0 - 1, iloc, jloc,&
+            & prow, pcol)
+        if (prow == mygrid%myrow .and. pcol == mygrid%mycol) then
+          localRef(i, j) = glob(iloc, jloc)
+        end if
+      end do
+    end do
+
+    localTest(:,:) = 0.0_dp
+    call scalafx_cpg2l(mygrid, desc, i0, j0, glob, localTest)
+
+    success = all(localTest == localRef)
+
+  end function readMatrixAndTest
+
+
+end program test_cpg2l
+