prg_graph_mod.F90

Go to the documentation of this file.

 
!
!
 
module prg_graph_mod
 
  use bml
  use prg_parallel_mod
  use omp_lib
 
  implicit none
 
  private  !Everything is private by default
 
  integer, parameter :: dp = kind(1.0d0)
 
  public :: subgraph_t
  public :: graph_partitioning_t
  public :: prg_initsubgraph
  public :: prg_destroysubgraph
  public :: prg_initgraphpartitioning
  public :: prg_destroygraphpartitioning
  public :: prg_printgraphpartitioning
  public :: prg_equalpartition
  public :: prg_equalgrouppartition
  public :: prg_filepartition
  public :: prg_fnormgraph
  public :: prg_sedacspartition
  public :: prg_setpartition
 
  type subgraph_t
 
    integer :: part
 
    integer :: hsize
 
    integer :: lsize
 
    integer :: llsize
 
    integer, allocatable :: core_halo_index(:)
 
    integer, allocatable :: nodeinpart(:)
 
    !    real(dp) :: vvx(100)
 
  end type subgraph_t
 
  type graph_partitioning_t
 
    character(len=100) :: pname
 
    integer :: myrank
 
    integer :: totalprocs
 
    integer :: totalparts
 
    integer :: totalnodes
 
    integer :: totalnodes2
 
    integer :: globalpartmin
 
    integer :: globalpartmax
 
    integer :: globalpartextent
 
    integer, allocatable :: localpartmin(:)
 
    integer, allocatable :: localpartmax(:)
 
    integer, allocatable :: localpartextent(:)
 
    integer, allocatable :: order(:)
 
    integer, allocatable :: reorder(:)
 
    integer :: nparts
 
    integer :: nx, ny, nz
 
    integer, allocatable :: nnodesinpart(:)
 
    integer, allocatable :: nnodesinpartall(:)
 
    integer :: pp(100)
 
    integer :: maxiter
 
    real(dp) :: ehomo
 
    real(dp) :: elumo
 
    real(dp) :: mineval
 
    real(dp) :: maxeval
 
    real(dp) :: vv(100)
 
    type (subgraph_t), allocatable :: sgraph(:)
 
  end type graph_partitioning_t
 
contains
 
  subroutine prg_initsubgraph(sg, pnum, hsize)
 
    type (subgraph_t), intent(inout) :: sg
    integer, intent(in) :: pnum, hsize
 
    sg%part = pnum
    sg%hsize = hsize
    sg%lsize = 0
    sg%llsize = 0
 
    allocate(sg%core_halo_index(hsize))
 
  end subroutine prg_initsubgraph
 
  subroutine prg_destroysubgraph(sg)
 
    type (subgraph_t), intent(inout) :: sg
 
    if (allocated(sg%core_halo_index)) deallocate(sg%core_halo_index)
    if (allocated(sg%nodeInPart)) deallocate(sg%nodeInPart)
 
  end subroutine prg_destroysubgraph
 
  subroutine prg_initgraphpartitioning(gp, pname, np, nnodes, nnodes2)
 
    type (graph_partitioning_t), intent(inout) :: gp
    character(len=*), intent(in) :: pname
    integer, intent(in) :: np, nnodes, nnodes2
 
    integer :: nprocs, avgparts
    integer :: i, it, last, nleft
 
    nprocs = getnranks()
    gp%myRank = getmyrank()
 
    !! Global
    gp%pname = pname
    gp%totalProcs = nprocs
    gp%totalParts = np
    gp%totalNodes = nnodes
    gp%totalNodes2 = nnodes2
 
    !! Global bounds
    gp%globalPartMin = 1
    gp%globalPartMax = np
    gp%globalPartExtent = gp%globalPartMax - gp%globalPartMin + 1
 
    allocate(gp%localPartMin(nprocs))
    allocate(gp%localPartMax(nprocs))
    allocate(gp%localPartExtent(nprocs))
 
    ! Distribute parts evenly among ranks
    avgparts = int(gp%totalParts / nprocs)
    do i = 1, nprocs
      gp%localPartExtent(i) = avgparts
    enddo
 
    nleft = gp%totalParts - nprocs * avgparts
    if (nleft .gt. 0) then
      do i = 1, nleft
        gp%localPartExtent(i) = gp%localPartExtent(i) + 1
      enddo
    endif
 
    gp%localPartMin(1) = 1
    gp%localPartMax(1) = gp%localPartExtent(1)
    do i = 2, nprocs
      gp%localPartMin(i) = gp%localPartMax(i-1) + 1
      gp%localPartMax(i) = gp%localPartMin(i) + gp%localPartExtent(i) - 1
    enddo
 
    gp%nparts = gp%localPartExtent(gp%myRank+1)
 
    allocate(gp%nnodesInPart(np))
    allocate(gp%nnodesInPartAll(np))
    allocate(gp%sgraph(np))
 
    !! For reordering
    allocate(gp%order(nnodes))
    allocate(gp%reorder(nnodes))
 
    gp%maxIter = 0
    gp%mineval = 0.0_dp
    gp%maxeval = 0.0_dp
 
    if (printrank() .eq. 1) then
      write(*,*)
      write(*,*) "total procs = ", gp%totalProcs
      write(*,*) "total nodes = ", gp%totalNodes
      write(*,*) "total nodes2 = ", gp%totalNodes2
      write(*,*) "total parts = ", gp%totalParts
      write(*,*) "local parts = ", gp%nparts
      write(*,*)
      write(*,*) "globalPartMin = ", gp%globalPartMin, &
           " globalPartMax = ", gp%globalPartMax, &
           " globalPartExtent = ", gp%globalPartExtent
 
      do i = 1, nprocs
        write(*,*) "rank = ", i-1, &
             " localPartMin = ", gp%localPartMin(i), &
             " localPartMax = ", gp%localPartMax(i), &
             " localPartExtent = ", gp%localPartExtent(i)
      enddo
 
      write(*,*)
    endif
 
  end subroutine prg_initgraphpartitioning
 
  subroutine prg_destroygraphpartitioning(gp)
 
    type (graph_partitioning_t), intent(inout) :: gp
 
    integer :: i
 
    if (allocated(gp%localPartMin))deallocate(gp%localPartMin)
    if (allocated(gp%localPartMax))deallocate(gp%localPartMax)
    if (allocated(gp%localPartExtent))deallocate(gp%localPartExtent)
 
    if (allocated(gp%order))deallocate(gp%order)
    if (allocated(gp%reorder))deallocate(gp%reorder)
 
    if(allocated(gp%nnodesInPart)) deallocate(gp%nnodesInPart)
    if(allocated(gp%nnodesInPartAll))deallocate(gp%nnodesInPartAll)
 
    if (allocated(gp%sgraph)) then
      do i = 1, gp%totalParts
        call prg_destroysubgraph(gp%sgraph(i))
      enddo
      deallocate(gp%sgraph)
    endif
 
  end subroutine prg_destroygraphpartitioning
 
  subroutine prg_printgraphpartitioning(gp)
 
    type (graph_partitioning_t), intent(in) :: gp
 
    integer :: i, j
 
    if (gp%myRank .ne. 0) return
 
    ! Global data
    write(*,*) ""
    write(*,*) ""
    write(*,*) "Graph partitioning:"
    write(*,*) ""
 
    write(*,*) "name = ", gp%pname
    write(*,*) "totalProcs = ", gp%totalProcs
    write(*,*) "totalParts = ", gp%totalParts
    write(*,*) "totalNodes = ", gp%totalNodes
    write(*,*) "totalNodes2 = ", gp%totalNodes2
    write(*,*) ""
    write(*,*) "globalPartMin = ", gp%globalPartMin, &
         " globalPartMax = ", gp%globalPartMax, &
         " globalPartExtent = ", gp%globalPartExtent
    write(*,*) ""
 
    ! Local data
    write(*,*) "local parts = ", gp%nparts
    do i = 1, gp%totalProcs
      write(*,*) "rank = ", i-1, " localPartMin = ", gp%localPartMin(i), &
           " localPartMax = ", gp%localPartMax(i), &
           " localPartExtent = ", gp%localPartExtent(i)
    enddo
    write(*,*) ""
 
    ! SP2 data
    write(*,*) "Number of iterations = ", gp%maxIter
    write(*,*) "SP2 sequence = ", (gp%pp(i),i=1,gp%maxIter)
    write(*,*) "mineval = ", gp%mineval , "maxeval = ", gp%maxeval
    write(*,*) ""
 
    ! For each subgraph
    do i = 1, gp%nparts
      write(*,*) "part = ", i, " hsize = ", gp%sgraph(i)%hsize, &
           " lsize = ", gp%sgraph(i)%lsize, &
           " llsize = ", gp%sgraph(i)%llsize
      write(*,*) ""
      write(*,*) "Number of core nodes in part = ", gp%nnodesInPart(i)
      write(*,*) "nodeInPart = ", &
           (gp%sgraph(i)%nodeInPart(j), j=1,gp%nnodesInPart(i))
      write(*,*) ""
 
      write(*,*) "core_halo_index = ", &
           (gp%sgraph(i)%core_halo_index(j), j=1,gp%sgraph(i)%lsize)
      write(*,*) ""
 
    enddo
 
  end subroutine prg_printgraphpartitioning
 
 
  subroutine prg_equalpartition(gp, nodesPerPart, nnodes)
 
    type (graph_partitioning_t), intent(inout) :: gp
    integer, intent(in) :: nodesperpart, nnodes
 
    integer :: i, j, it, np, psize
    character(len=100) :: pname
 
    !! Init graph partitioning
    if (nnodes.le.nodesperpart)then
      np = 1
    else
      np = ceiling(real(nnodes) / real(nodesperpart))
    endif
    pname = '("equalParts")'
    call prg_destroygraphpartitioning(gp)
    call prg_initgraphpartitioning(gp, pname, np, nnodes, nnodes)
 
    ! Assign node ids (mapped to orbitals as rows) to each node in each
    ! partition
    !$omp parallel do default(none) private(i) &
    !$omp private(it,j,psize) &
    !$omp shared(gp,nnodes,nodesPerPart)
    do i = 1, gp%totalParts
      call prg_initsubgraph(gp%sgraph(i), i, nnodes)
      if ((i * nodesperpart) .le. nnodes) then
        psize = nodesperpart
      else
        psize = nnodes - (nodesperpart * (i-1))
      endif
      allocate(gp%sgraph(i)%nodeInPart(psize))
      do j = 1, psize
        it = (i-1) * nodesperpart + j-1;
        gp%sgraph(i)%nodeInPart(j) = it
      enddo
      gp%nnodesInPart(i) = psize
      gp%nnodesInPartAll(i) = psize
    enddo
    !$omp end parallel do
 
  end subroutine prg_equalpartition
 
  function time2milliseconds() result(mls)
 
    real(8) :: mls
    integer :: timevector(8)
 
    call date_and_time(values=timevector)
    mls = timevector(5)*60*60*1000 + timevector(6)*60*1000 + &
         timevector(7)*1000 + timevector(8)
 
  end function time2milliseconds
 
 
  subroutine prg_setpartition(gp,whichParts,nparts,nnodes,verb)
    implicit none
    type (graph_partitioning_t), intent(inout) :: gp
    integer, intent(in) :: nnodes, nparts
    integer, allocatable, intent(inout) :: whichparts(:)
    integer, optional, intent(in) :: verb
    integer :: i, cnt, part, j, it, psize, cut, cutold, ac, nac
    integer ::  mdim, nodesperpart,tobedistributed,sumsizes,rem
    character(len=100) :: pname
    logical :: writeout
 
    if(present(verb))then
      if(verb >= 1)then
        write(*,*)"Setting partitioning in prg_setPartition..."
        if(verb >= 2) writeout = .true.
      endif
    endif
 
    pname = '("setPartition")'
    call prg_destroygraphpartitioning(gp)
    call prg_initgraphpartitioning(gp, pname, nparts, nnodes, nnodes)
 
    !Assign node ids to each node in each part
    !$omp parallel do default(none) private(i) &
    !$omp private(cnt,it,j,psize,part) &
    !$omp shared(gp,nnodes,whichParts)
    do i = 1, gp%totalParts
      call prg_initsubgraph(gp%sgraph(i), i, nnodes)
      psize = int(sum(whichparts, mask=(whichparts==i))/i)
      allocate(gp%sgraph(i)%nodeInPart(psize))
      cnt = 0
      do j = 1,nnodes
        part = whichparts(j)
        if(part == i)then
          cnt = cnt + 1
          gp%sgraph(i)%nodeInPart(cnt) = j - 1
        endif
      enddo
      gp%nnodesInPart(i) = psize
      gp%nnodesInPartAll(i) = psize
    enddo
    !$omp end parallel do
 
end subroutine prg_setpartition
  subroutine prg_sedacspartition(gp,coords,whichParts_guess_saved,g_bml,nparts,nnodes,verb)
    implicit none
    type (graph_partitioning_t), intent(inout) :: gp
    integer, intent(in) :: nnodes, nparts
    integer, allocatable, intent(inout) :: whichparts_guess_saved(:)
    integer, optional, intent(in) :: verb
    integer :: i, cnt, part, j, it, psize, cut, cutold, ac, nac
    integer ::  mdim, nodesperpart,tobedistributed,sumsizes,rem
    integer, allocatable :: whichparts(:)
    integer, allocatable :: graph(:,:),degs(:)
    real(kind(1.0)), allocatable :: row(:)
    real(dp) :: bal,sumdegs, mlsi
    real(dp), allocatable, intent (in) :: coords(:,:)
    character(len=100) :: pname
    type (bml_matrix_t), intent(in) :: g_bml
    logical :: writeout
 
    if(present(verb))then
      if(verb >= 1)then
        write(*,*)"Doing SEDACS type of partitioning ..."
        if(verb >= 2) writeout = .true.
      endif
    endif
 
    !If there is no guess we will do a "chunk/block" partitioning
    !This initial pratitioning has all the nodes evenly distributed
    !if(allocated(whichParts_guess_saved))deallocate(whichParts_guess_saved)
    if(.not. allocated(whichparts_guess_saved))then
      nodesperpart = int(nnodes/nparts)
      rem = nnodes - nodesperpart*nparts
      tobedistributed = rem
 
      pname = '("equalParts")'
      call prg_destroygraphpartitioning(gp)
      call prg_initgraphpartitioning(gp, pname, nparts, nnodes, nnodes)
 
      sumsizes = 0
      do i = 1, gp%totalParts
        call prg_initsubgraph(gp%sgraph(i), i, nnodes)
        if (sumsizes .le. nnodes - nodesperpart) then
          psize = nodesperpart + min(tobedistributed,1)
          tobedistributed = max(tobedistributed - 1,0)
        else
          psize = nnodes - sumsizes
        endif
        write(*,*)i,psize,sumsizes
        allocate(gp%sgraph(i)%nodeInPart(psize))
        do j = 1, psize
          it = sumsizes + j-1;
          gp%sgraph(i)%nodeInPart(j) = it
        enddo
        sumsizes = sumsizes + psize
        gp%nnodesInPart(i) = psize
        gp%nnodesInPartAll(i) = psize
      enddo
 
      call prg_get_parts_indices(gp,nnodes,whichparts)
      allocate(whichparts_guess_saved(nnodes))
    else
      allocate(whichparts(nnodes))
      whichparts = whichparts_guess_saved
 
    endif
 
    !Get degrees and build 2D-sedacs-like graph
    !In this case the degreess are saved in a separate array
    mdim = bml_get_m(g_bml)
    allocate(graph(nnodes,mdim))
    allocate(degs(nnodes))
    allocate(row(nnodes))
    row = 0.0
    degs = 0
    do i = 1,nnodes
      call bml_get_row(g_bml,i,row)
      do j = 1,nnodes
        if(row(j) .ge. 0.5)then
          degs(i) = degs(i) + 1
          graph(i,degs(i)) = j
        endif
      enddo
    enddo
    deallocate(row)
 
    sumdegs = sum(degs)
    mlsi = time2milliseconds()
    call prg_get_graphcut(whichparts,degs,graph,cut) !Get the initial cut
    write(*,*)"Time for initial cut",time2milliseconds() - mlsi
    if(writeout) write(*,*)"Iter, Cut, RelCut, and Bal",0,cut, cut/sumdegs,bal
    cutold = cut
    do i = 1,200
      mlsi = time2milliseconds()
      call do_flips_precomp(whichparts,coords,graph,g_bml,degs,nnodes,nparts,bal)
      write(*,*)"Time for real flip",time2milliseconds() - mlsi
      mlsi = time2milliseconds()
      call prg_get_graphcut(whichparts,degs,graph,cut)
      write(*,*)"Time for real cut",time2milliseconds() - mlsi
      mlsi = time2milliseconds()
      call prg_get_balancing(whichparts,nparts,bal)
      write(*,*)"Time for real bal",time2milliseconds() - mlsi
      if(writeout) write(*,*)"Iter, Cut, RelCut, and Bal",i,cut, cut/sumdegs,bal
      if(cut >= cutold)then
        exit
      else
        cutold = cut
      endif
    enddo
 
    deallocate(graph)
    deallocate(degs)
 
    pname = '("sedacsPartition")'
    call prg_destroygraphpartitioning(gp)
    call prg_initgraphpartitioning(gp, pname, nparts, nnodes, nnodes)
 
    !Assign node ids to each node in each part
    !$omp parallel do default(none) private(i) &
    !$omp private(cnt,it,j,psize,part) &
    !$omp shared(gp,nnodes,whichParts)
    do i = 1, gp%totalParts
      call prg_initsubgraph(gp%sgraph(i), i, nnodes)
      psize = int(sum(whichparts, mask=(whichparts==i))/i)
      allocate(gp%sgraph(i)%nodeInPart(psize))
      cnt = 0
      do j = 1,nnodes
        part = whichparts(j)
        if(part == i)then
          cnt = cnt + 1
          gp%sgraph(i)%nodeInPart(cnt) = j - 1
        endif
      enddo
      gp%nnodesInPart(i) = psize
      gp%nnodesInPartAll(i) = psize
    enddo
    !$omp end parallel do
 
    whichparts_guess_saved = whichparts
    deallocate(whichparts)
 
  end subroutine prg_sedacspartition
 
 
  !! Do node partition flips with precomputed cuts.
  !! \brief This function is a special case of do_flips where
  !! the cuts around a node are precomputed for all possible
  !! part index that same node could have. This will differ from the do_flips
  !! since everytime there is a flip, there is no actualization of the cuts. The
  !! price to pay is the need of more iterations until convergence.
  !! \param whichPart partition indexing vector.
  !! \param graph Graph to be partition. graph[i,0] = degree of node i.
  !! graph[i,j>0] = the node conected to node i.
  !! \param nnodes Number of nodes.
  !! \param nparts Number of parts.
  !! \return whichPartNew New partition indexing verctor.
  !
  subroutine do_flips_precomp(whichParts,coords,graph,adj,degs,nnodes,nparts,bal)
    implicit none
    integer, allocatable :: cutsI(:,:)
    integer, intent(in) :: nnodes,nparts
    integer, allocatable, intent(inout) :: graph(:,:),degs(:)
    integer, allocatable, intent(inout) :: whichParts(:)
    real(dp) :: bal
    real(kind(1.0)), allocatable :: row(:)
    type(bml_matrix_t) :: adj
    integer :: i,j,deg,ind,ii,origCut,newCut, sumDeltaCut,ac,nac,newCut1
    integer :: origCutI,newCutI,origCutJ,newCutJ
    integer :: partIndexI, partIndexII,partIndexJ,k,cut,cutOld,newPartIndex,inDefect,discon
    integer, allocatable :: momentI(:),momentsI(:,:),partSizes(:)
    real(dp), allocatable, intent (in) :: coords(:,:)
    logical :: disconnected, conditionToFlip,actualize
 
 
    allocate(cutsi(nnodes,nparts))
    cutsi = 0
    sumdeltacut = 0
    ac = 0
    nac = 0
    indefect = 0
 
    !$omp parallel do default(none) private(i) &
    !$omp private(ii,deg,partIndexII,ind) &
    !$omp shared(degs,graph,nnodes,whichParts,cutsI)
    do i = 1,nnodes
      deg = degs(i)
      !Get the max cut a node could have
      cutsi(i,:) = deg
      !Lets look at every neighbor
      do ii = 1,deg
        ind = graph(i,ii)
        partindexii = whichparts(ind)
        !Everytime there is a neighbor in a certain part
        !it will decrese the cut of I if I would be on that
        !same part.
        cutsi(i,partindexii) = cutsi(i,partindexii) - 1
      enddo
    enddo
    !$omp end parallel do
 
    !Get sizes
    allocate(row(nnodes))
    allocate(partsizes(nparts))
    partsizes = 0
 
    do i = 1,nnodes
      partindexi = whichparts(i)
      partsizes(partindexi) = partsizes(partindexi) + 1
    enddo
 
    actualize = .true.
 
!    !$omp parallel do default(none) private(i) &
!    !$omp private(row) &
!    !$omp private(j,partIndexI,partIndexJ,origCut,origCutI,origCutJ) &
!    !$omp private(newPartIndex,newCut,newCutI,newCutJ,ind,newCut1,inDefect) &
!    !$omp shared(adj,whichParts,cutsI,graph,partSizes,actualize) &
!    !$omp shared(degs,nnodes)
    do i = 1,nnodes
      call bml_get_row(adj,i,row)
      if(cutsi(i,whichparts(i)) > 0)then
        do j = i+1, nnodes
        !do j = 1,nnodes
          !if(row(j) > 0.5)then
          if((whichparts(i) .ne. whichparts(j)))then
            !conditionToFlip = .false.
            !disconnected = .false.
 
            !Look at their neighbors and count the cuts
            partindexi = whichparts(i)
            partindexj = whichparts(j)
            origcuti = cutsi(i,partindexi)
            origcutj = cutsi(j,partindexj)
 
            newcuti = cutsi(i,partindexj)
 
            if(newcuti < degs(i))then
              newcutj = cutsi(j,partindexi)
              if(newcutj < degs(j))then
 
                !Now we know the cut when I is in partIndexI and J
                origcut = origcuti + origcutj
                newcut = newcuti + newcutj
 
                if((newcut <= origcut))then
 
          !        !$OMP CRITICAL      
                  !ac = ac + 1
                  whichparts(i) = partindexj
                  whichparts(j) = partindexi
 
                  !Actualizing possible cuts of neighbors of i and j
                  if(actualize)then
                  do ii = 1,degs(i)
                    ind = graph(i,ii)
                    !The neighs of i "if now they are in the "new color
                    !of i", their cut "at that color" will be decreased
                    !by one.
                    cutsi(ind,partindexj) = cutsi(ind,partindexj) - 1
                    cutsi(ind,partindexi) = cutsi(ind,partindexi) + 1
                  enddo
                  do ii = 1,degs(j)
                    ind = graph(j,ii)
                    !Same for the neighs of j
                    cutsi(ind,partindexi) = cutsi(ind,partindexi) - 1
                    cutsi(ind,partindexj) = cutsi(ind,partindexj) + 1
                  enddo
                  endif
          !        !$OMP END CRITICAL 
                endif
              endif
            else
              !nac = nac + 1
            endif
 
          elseif(whichparts(i) .eq. whichparts(j))then
            if(row(j) < 0.5)then
              origcut = 0
              newcut = 0
              partindexi = whichparts(i)
              partindexj = whichparts(j)
 
              origcut = cutsi(i,partindexi)
              origcut = origcut + cutsi(j,partindexj)
 
              ! Change color of I
              indefect = minloc(partsizes,dim=1)
 
              newpartindex = indefect
 
              newcut = cutsi(i,newpartindex)
              newcut = newcut + cutsi(j,partindexj)
 
              newcut1 = cutsi(i,partindexi)
              newcut1 = newcut1 + cutsi(j,newpartindex)
 
              if(newcut < origcut)then
                whichparts(i) = newpartindex
                whichparts(j) = partindexj
                if(newpartindex .ne.partindexj)then
                  partsizes(newpartindex) = partsizes(newpartindex) + 1
                  partsizes(partindexj) = partsizes(partindexj) - 1
                endif
                !sumDeltaCut = newCut - origCut + sumDeltaCut
              elseif(newcut1 < origcut)then
                whichparts(i) = partindexi
                whichparts(j) = newpartindex
                if(newpartindex .ne. partindexi)then
                  partsizes(newpartindex) = partsizes(newpartindex) + 1
                  partsizes(partindexj) = partsizes(partindexj) - 1
                endif
                !sumDeltaCut = newCut1 - origCut + sumDeltaCut
              else
                whichparts(i) = partindexi
                whichparts(j) = partindexj
              endif
            endif
          else
            cycle
          endif
          !endif
        enddo
      endif
    enddo
!    !$omp end parallel do
 
    deallocate(partsizes)
    deallocate(cutsi)
    deallocate(row)
 
  end subroutine do_flips_precomp
 
 
 
  subroutine prg_get_balancing(whichParts,np,bal)
    implicit none
    integer, allocatable, intent(in) :: whichParts(:)
    integer, allocatable :: partsSizes(:)
    integer :: np, i, nnodes
    real(dp), intent(out) :: bal
 
    allocate(partssizes(np))
    nnodes = size(whichparts)
    partssizes = 0
 
!    !$omp parallel do default(none) private(i) &
!    !$omp private(nnodes) &
!    !$omp shared(whichParts,np) &
!    !$omp reduction(+:partsSizes)
    do i = 1,nnodes
      !write(*,*)"whichParts",i,nnodes,whichParts(i),np
      if(whichparts(i) > np)then
        write(*,*)"!!!ERROR: Part index",whichparts(i),&
             "is larger than the total number of parts,",np
        stop
      endif
      partssizes(whichparts(i)) = partssizes(whichparts(i)) + 1
    enddo
!    !$omp end parallel do
 
    bal = real(maxval(partssizes))/real(minval(partssizes))
    deallocate(partssizes)
  end subroutine prg_get_balancing
 
  subroutine prg_get_parts_indices(gp,nnodes,whichParts)
    type (graph_partitioning_t), intent(in) :: gp
    integer, allocatable,intent(out) :: whichParts(:)
    integer, intent(in) :: nnodes
    integer :: i,j,k,kk,partIndex
 
    if(.not. allocated(whichparts)) allocate(whichparts(nnodes))
 
    partindex = 0
    whichparts = 0
 
    !$omp parallel do default(none) private(i) &
    !$omp private(k,kk) &
    !$omp shared(gp,nnodes,whichParts)
    do i = 1,gp%totalParts !Loop over parts
      do j = 1,gp%nnodesInPart(i)  !Loop over nodes in part
        k = gp%sgraph(i)%nodeInPart(j)
        kk = k + 1 !Indexig starts from 0 for all sgraphs (beacuse of historical reasons)
        whichparts(kk) = i
      enddo
    enddo
    !$omp end parallel do
 
  end subroutine prg_get_parts_indices
 
 
  subroutine prg_get_graphcut(whichParts,degs,graph,cut)
    integer, allocatable, intent(in) :: whichParts(:)
    integer, intent(out) :: cut
    integer :: nnodes,i,j,partIndexI,partIndexJ,k,kk
    integer :: ind
    integer, allocatable, intent(in) :: graph(:,:),degs(:)
 
    cut = 0
    nnodes = size(whichparts)
 
    do i = 1,nnodes
      !write(*,*)"graph",graph(i,1:degs(i))
      partindexi = whichparts(i)
      !Look at the neighbors and check if they are in different part
      do j = 1,degs(i)
        ind = graph(i,j)
        partindexj = whichparts(ind)
        if(int(partindexi - partindexj) .ne. 0)then
          cut = cut + 1
        endif
      enddo
    enddo
    cut = 0.5_dp*cut !Because of the double counting
 
  end subroutine prg_get_graphcut
 
 
  subroutine prg_equalgrouppartition(gp, hindex, ngroup, nodesPerPart, nnodes)
 
    type (graph_partitioning_t), intent(inout) :: gp
    integer, intent(in) :: nodesperpart, nnodes, ngroup
    integer, intent(in) :: hindex(2,ngroup)
 
    integer :: i, j, k, ll, np, psize, ind, ptotal
    character(len=100) :: pname
 
    !! Init graph partitioning
    np = ceiling(real(ngroup) / real(nodesperpart))
    write(pname, '("equalGroupParts")')
    call prg_destroygraphpartitioning(gp)
    call prg_initgraphpartitioning(gp, pname, np, ngroup, nnodes)
 
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    !! partition
    !$omp parallel do default(none) &
    !$omp private(i, j, k, ll, ind, psize, ptotal) &
    !$omp shared(gp, hindex, nnodes, ngroup, nodesPerPart)
    do i = 1, gp%totalParts
      call prg_initsubgraph(gp%sgraph(i), i, nnodes)
 
      !! Figure out number of groups in part
      if ((i * nodesperpart) .le. ngroup) then
        psize = nodesperpart
      else
        psize = ngroup - (nodesperpart * (i-1))
      endif
 
      !! Figure out total nodes/rows in part
      ind = (i-1)*nodesperpart
      ptotal = 0
      do j = 1, psize
        ptotal = ptotal + hindex(2, ind+j) - hindex(1, ind+j) + 1
      enddo
      gp%nnodesInPart(i) = ptotal
      gp%nnodesInPartAll(i) = ptotal
 
      !! Enumerate all nodes in part
      allocate(gp%sgraph(i)%nodeInPart(ptotal))
      ll = 1
      do j = 1, psize
        do k = hindex(1, ind+j), hindex(2, ind+j)
          gp%sgraph(i)%nodeInPart(ll) = k-1
          ll = ll + 1
        enddo
      enddo
    enddo
    !$omp end parallel do
 
    !    do i = 1, gp%totalParts
    !      write(*,*) "part ", i, ": ", gp%nnodesInPart(i), " nodes"
    !      write(*,*) "     ", &
    !        (gp%sgraph(i)%nodeInPart(ll),ll = 1, gp%nnodesInPart(i))
    !    enddo
 
  end subroutine prg_equalgrouppartition
 
  subroutine prg_filepartition(gp, partFile)
 
    type (graph_partitioning_t), intent(inout) :: gp
    character(len=*), intent(in) :: partfile
 
    call prg_readpart(gp, partfile)
 
  end subroutine prg_filepartition
 
  subroutine prg_readpart(gp, partFile)
 
    character(len=*), intent(in) :: partFile
    type (graph_partitioning_t), intent(inout) :: gp
 
    integer :: pfile
    integer :: totalNodes, totalParts
    integer :: i, j, ip, pnode
    character(len=100) :: pname
 
    pfile = 10
    open(unit=pfile, status='old', file=partfile)
 
    read(pfile, *) pname
    read(pfile, *) totalnodes, totalparts
 
    call prg_initgraphpartitioning(gp, pname, totalparts, totalnodes, totalnodes)
 
    !! Read in part sizes
    read(pfile, *) (gp%nnodesInPartAll(i), i=1,gp%totalParts)
    do i = 1, gp%totalParts
      gp%nnodesInPart(i) = gp%nnodesInPartAll(i)
    enddo
 
    !! Read in nodes for each part
    do i = 1, gp%totalParts
      read(pfile, *) ip
      call prg_initsubgraph(gp%sgraph(i), i, totalnodes)
      allocate(gp%sgraph(i)%nodeInPart(gp%nnodesInPart(i)))
      read(pfile, *) (gp%sgraph(i)%nodeInPart(j),j=1,gp%nnodesInPart(i))
      do j = 1, gp%nnodesInPart(i)
        gp%sgraph(i)%nodeInPart(j) = gp%sgraph(i)%nodeInPart(j) + 1
      enddo
    enddo
 
    close(pfile)
 
  end subroutine prg_readpart
 
  subroutine prg_fnormgraph(gp)
 
    type(graph_partitioning_t), intent(inout) :: gp
 
    integer :: i, j
#ifdef DO_MPI
    real(dp), allocatable :: slocal(:), sglobal(:)
#endif
 
#ifdef DO_MPI
    ! Sum traces from all parts on all ranks
    if (getnranks() .gt. 1) then
      allocate(slocal(gp%maxIter))
      allocate(sglobal(gp%maxIter))
      do i = 1, gp%maxIter
        slocal(i) = gp%vv(i)
      enddo
      call sumrealparallel(slocal, sglobal, gp%maxIter);
      do i = 1, gp%maxIter
        gp%vv(i) = sglobal(i)
      enddo
      deallocate(slocal, sglobal)
    endif
#endif
 
    !! Take sqrt for fnorm per iter
    do i = 1, gp%maxIter
      gp%vv(i) = sqrt(gp%vv(i))
    enddo
 
    if (printrank() .eq. 1) then
      write(*,*)
      write(*,*) "prg_fnormGraph:"
      do i = 1, gp%maxIter
        write(*,*) "iter = ", i, " fnorm = ", gp%vv(i)
      enddo
    endif
 
  end subroutine prg_fnormgraph
 
end module prg_graph_mod