prg_subgraphloop_mod.F90

Go to the documentation of this file.

 
!
!
 
module prg_subgraphloop_mod
 
  use prg_graph_mod
  use prg_sp2_mod
  use prg_timer_mod
  use prg_parallel_mod
  use bml
  use omp_lib
 
  implicit none
 
  private  !Everything is private by default
 
  integer, parameter :: dp = kind(1.0d0)
 
  public :: prg_subgraphsp2loop
  public :: prg_balanceparts
  public :: prg_partordering
  public :: prg_getpartitionhalosfromgraph
  public :: prg_getgrouppartitionhalosfromgraph
  public :: prg_collectmatrixfromparts
 
contains
 
  !! Subgraph SP2 loop.
  !!
  !! \param h_bml Input Hamiltonian matrix
  !! \param g_bml Input Graph as matrix
  !! \param rho_bml Output density matrix
  !! \param gp Input/Output graph partitioning
  !! \param threshold threshold for matrix operations
  subroutine prg_subgraphsp2loop(h_bml, g_bml, rho_bml, gp, threshold)
 
    type (bml_matrix_t), intent(in) :: h_bml, g_bml
    type (bml_matrix_t), intent(inout) :: rho_bml
    type (graph_partitioning_t), intent(inout) :: gp
    real(dp), intent(in) :: threshold
 
    integer :: i, j, k
    integer, allocatable :: vsize(:)
    type (bml_matrix_t) :: x_bml
    real(dp) :: thresh0
 
    allocate(vsize(2))
 
    thresh0 = 0.0_dp
 
    call prg_timer_start(subind_timer)
 
    ! Determine elements for each subgraph
    !$omp parallel do default(none) &
    !$omp private(i, vsize) &
    !$omp shared(gp, h_bml, g_bml)
    do i = 1, gp%totalParts
 
      call bml_matrix2submatrix_index(g_bml, &
           gp%sgraph(i)%nodeInPart, gp%nnodesInPart(i), &
           gp%sgraph(i)%core_halo_index, &
           vsize, .true., h_bml)
 
      gp%sgraph(i)%lsize = vsize(1)
      gp%sgraph(i)%llsize = vsize(2)
 
    enddo
    !$omp end parallel do
 
    !    do i = 1, gp%totalParts
    !        write(*,*)"i = ", i, " core size = ", gp%sgraph(i)%llsize, &
    !                  " full size = ", gp%sgraph(i)%lsize
    !        write(*,*)"nnodes = ", gp%nnodesInPart(i), &
    !                  (gp%sgraph(i)%nodeInPart(k),k=1,gp%nnodesInPart(i))
    !    enddo
 
    call prg_timer_stop(subind_timer)
 
    deallocate(vsize)
 
    ! Balance parts by size of subgraph
    if (getnranks() > 1) then
      call prg_balanceparts(gp)
      call prg_partordering(gp)
    endif
 
    ! Process each part one at a time
    !do i = 1, gp%nparts
    do i = gp%localPartMin(gp%myRank+1), gp%localPartMax(gp%myRank+1)
 
      call bml_zero_matrix(bml_matrix_dense, bml_element_real, dp, &
           gp%sgraph(i)%lsize, gp%sgraph(i)%lsize, x_bml);
 
      ! Extract subgraph and prg_normalize
      call prg_timer_start(subext_timer)
      call bml_matrix2submatrix(h_bml, x_bml, &
           gp%sgraph(i)%core_halo_index, gp%sgraph(i)%lsize)
      call prg_timer_stop(subext_timer)
 
      ! Run SP2 on subgraph/submatrix
      call prg_timer_start(subsp2_timer)
      !call prg_sp2_submatrix_inplace(x_bml, threshold, gp%pp, &
      call prg_sp2_submatrix_inplace(x_bml, thresh0, gp%pp, &
           gp%maxIter, gp%vv, gp%mineval, gp%maxeval, &
           gp%sgraph(i)%llsize)
      call prg_timer_stop(subsp2_timer)
 
      ! Reassemble into density matrix
      call prg_timer_start(suball_timer)
      call bml_submatrix2matrix(x_bml, rho_bml, &
           gp%sgraph(i)%core_halo_index, gp%sgraph(i)%lsize, &
           gp%sgraph(i)%llsize, threshold)
      call prg_timer_stop(suball_timer)
 
      call bml_deallocate(x_bml)
 
    enddo
 
    ! Fnorm
    call prg_fnormgraph(gp)
 
    ! Collect density matrix over local and distributed parts
    call prg_collectmatrixfromparts(gp, rho_bml)
 
  end subroutine prg_subgraphsp2loop
 
  subroutine prg_collectmatrixfromparts(gp, rho_bml)
 
    implicit none
 
    type (graph_partitioning_t), intent(inout) :: gp
    type (bml_matrix_t), intent(inout) :: rho_bml
 
    if (getnranks() > 1) then
 
      ! Save original domain
      call bml_save_domain(rho_bml)
 
      ! Update density matrix domain based on partitions of orbitals
      ! for gather
      call bml_update_domain(rho_bml, gp%localPartMin, gp%localPartMax, &
           gp%nnodesInPart)
      call bml_reorder(rho_bml, gp%reorder);
 
      ! Exchange/gather density matrix across ranks
      call bml_allgathervparallel(rho_bml)
 
      ! Reorder density matrix to match original
      call bml_reorder(rho_bml, gp%order)
 
      ! Restore to original domain
      call bml_restore_domain(rho_bml)
 
    endif
 
  end subroutine prg_collectmatrixfromparts
 
  !! Balance subgraphs by size.
  subroutine prg_balanceparts(gp)
 
    type (graph_partitioning_t), intent(inout) :: gp
 
    type (subgraph_t), allocatable :: temp_sgraph(:)
    type (subgraph_t) :: temp
    integer :: nranks, myrank, nparts, avgparts
    integer :: i, j, rid, sid, onum, pnode, ip
 
    integer :: tempp
    integer, allocatable :: tsize(:), torder(:)
 
    nranks = getnranks()
    myrank = getmyrank()
 
    nparts = gp%totalParts
    avgparts = nparts / nranks
 
    ! Sort parts by core+halo size
    allocate(temp_sgraph(nparts))
    allocate(tsize(nparts))
    allocate(torder(nparts))
 
    do i = 1, nparts
      temp_sgraph(i) = gp%sgraph(i)
      torder(i) = i
      tsize(i) = temp_sgraph(i)%llsize
    enddo
 
    !! Sort parts by core+halo size
    do i = 1, nparts-1
      do j = i, nparts
        if (tsize(i) .lt. tsize(j)) then
          tempp = tsize(i)
          tsize(i) = tsize(j)
          tsize(j) = tempp
 
          tempp = torder(i)
          torder(i) = torder(j)
          torder(j) = tempp
        endif
      enddo
    enddo
 
    !! Print ordered subgraph sizes after sorting
    !   if (printRank() .eq. 1) then
    !      write(*,*) "after sort:"
    !      do i = 1, nParts
    !        write(*,*)"order ", i, " part = ", temp_sgraph(i)%part, " size = ", &
    !                  temp_sgraph(i)%lsize
    !      enddo
    !    endif
 
    rid = 1
    do i = 1, nranks
      sid = sid + 1
      do j = i, nparts, nranks
        gp%sgraph(rid) = temp_sgraph(torder(j))
        !write(*,*) "rank = ", i, " part = "), rid, " from = ", j
        rid = rid + 1
        sid = sid + 1
      enddo
      if (sid .le. gp%localPartExtent(i)) then
        gp%sgraph(rid) = temp_sgraph(torder(nranks*avgparts+i))
        !write(*,*) "rank = ", i, " part = ", rid, " from = ", nranks*avgparts+i
        rid = rid + 1
      endif
    enddo
 
    !! Print ordered subgraph sizes after renumbering
    !    if (printRank() .eq. 1) then
    !      write(*,*) "after renumber:"
    !      do i = 1, nParts
    !        write(*,*)"order ", i, " part = ", gp%sgraph(i)%part, &
    !                  " core = ", gp%sgraph(i)%llsize, &
    !                  " size = ", gp%sgraph(i)%lsize, &
    !                  " min node = ", minval(gp%sgraph(i)%nodeInPart) + 1, &
    !                  " max node = ", maxval(gp%sgraph(i)%nodeInPart) + 1
    !      enddo
    !    endif
 
    !! Reset number of nodes per part
    do i = 1, nparts
      gp%nnodesInPart(i) = gp%sgraph(i)%llsize
      gp%nnodesInPartAll(i) = gp%sgraph(i)%llsize
    enddo
 
    deallocate(temp_sgraph)
    deallocate(tsize)
    deallocate(torder)
 
  end subroutine prg_balanceparts
 
  subroutine prg_partordering(gp)
 
    type (graph_partitioning_t), intent(inout) :: gp
 
    integer :: i, j, onum, pnode, nparts
 
    nparts = gp%totalParts
 
    !! Set ordering and reordering based on partitions
    onum = 0
    do i = 1, nparts
      !!write(*,*) i, " Part ", ip, " nnodes = ", gp%sgraph(i)%llsize
      do j = 1, gp%sgraph(i)%llsize
        pnode = gp%sgraph(i)%nodeInPart(j)
        !!write(*,*) "  ", j, "  Node ", pnode, " to ", onum
        gp%reorder(pnode+1) = onum
        gp%order(onum+1) = pnode
        onum = onum + 1
      enddo
    enddo
 
  end subroutine prg_partordering
 
  subroutine prg_getgrouppartitionhalosfromgraph(gp, g_bml, hnode, djflag)
 
    implicit none
 
    type (graph_partitioning_t), intent(inout) :: gp
    type (bml_matrix_t), intent(in) :: g_bml
    integer, intent(in) :: hnode(*)
    logical, intent(in) :: djflag
 
    integer :: i, j, inode
    integer, allocatable :: vsize(:)
 
    allocate(vsize(2))
 
    !$omp parallel do default(none) &
    !$omp private(i, j, inode, vsize) &
    !$omp shared(gp, g_bml, hnode, djflag)
    do i = 1, gp%totalParts
 
      do j = 1, gp%nnodesInPart(i)
        inode = gp%sgraph(i)%nodeInPart(j)
        gp%sgraph(i)%nodeInPart(j) = hnode(inode+1) - 1
      enddo
 
      call bml_matrix2submatrix_index(g_bml, &
           gp%sgraph(i)%nodeInPart, gp%nnodesInPart(i), &
           gp%sgraph(i)%core_halo_index, &
           vsize, djflag)
 
      gp%sgraph(i)%lsize = vsize(1)
      gp%sgraph(i)%llsize = vsize(2)
 
    enddo
    !$omp end parallel do
 
    deallocate(vsize)
 
  end subroutine prg_getgrouppartitionhalosfromgraph
 
  subroutine prg_getpartitionhalosfromgraph(gp, g_bml, djflag)
 
    implicit none
 
    type (graph_partitioning_t), intent(inout) :: gp
    type (bml_matrix_t), intent(in) :: g_bml
    logical, intent(in) :: djflag
 
    integer :: i
    integer, allocatable :: vsize(:)
 
    allocate(vsize(2))
 
    !$omp parallel do default(none) &
    !$omp private(i, vsize) &
    !$omp shared(gp, g_bml, djflag)
    do i = 1, gp%totalParts
 
      call bml_matrix2submatrix_index(g_bml, &
           gp%sgraph(i)%nodeInPart, gp%nnodesInPart(i), &
           gp%sgraph(i)%core_halo_index, &
           vsize, djflag)
 
      gp%sgraph(i)%lsize = vsize(1)
      gp%sgraph(i)%llsize = vsize(2)
 
    enddo
    !$omp end parallel do
 
    deallocate(vsize)
 
  end subroutine prg_getpartitionhalosfromgraph
 
end module prg_subgraphloop_mod