prg_partition_mod.F90

Go to the documentation of this file.

 
module prg_partition_mod
 
  use bml
  use prg_graph_mod
  use prg_parallel_mod
  use prg_extras_mod
  use, intrinsic :: iso_c_binding
 
  implicit none
 
  private  !Everything is private by default
 
  integer, parameter :: dp = kind(1.0d0)
 
  integer, parameter :: metis_index_kind = metis_index_kind
 
  integer, parameter :: metis_real_kind = kind(metis_real_kind)
 
  public :: prg_metispartition
  public :: prg_costpartition
  public :: update_prg_costpartition
  public :: prg_simannealing
  public :: prg_check_arrays
  public :: prg_kernlin
  public :: prg_kernlin2
  public :: prg_kernlin_queue
  public :: prg_update_gp
  public :: prg_simannealing_old
 
#ifdef DO_GRAPHLIB
  interface
 
    integer function metis_setdefaultoptions(options) &
         bind(C, name="METIS_SetDefaultOptions")
 
      import metis_index_kind
 
      integer(kind=metis_index_kind), intent(in) :: options(*)
 
    end function metis_setdefaultoptions
 
    integer function metis_partgraphkway(nvtxs, ncon, xadj, adjncy, vwgt, &
         vsize, adjwgt, nparts, tpwgts, ubvec, options, objval, part) &
         bind(C, name="METIS_PartGraphKway")
 
      import metis_index_kind
      import metis_real_kind
 
      integer(kind=metis_index_kind), intent(in) :: nvtxs(*)
      integer(kind=metis_index_kind), intent(in) :: ncon(*)
      integer(kind=metis_index_kind), intent(in) :: xadj(*)
      integer(kind=metis_index_kind), intent(in) :: adjncy(*)
      integer(kind=metis_index_kind), intent(in) :: vwgt(*)
      integer(kind=metis_index_kind), intent(in) :: vsize(*)
      integer(kind=metis_index_kind), intent(in) :: adjwgt(*)
      integer(kind=metis_index_kind), intent(in) :: nparts(*)
      real(kind=metis_real_kind), intent(in) :: tpwgts(*)
      real(kind=metis_real_kind), intent(in) :: ubvec(*)
      integer(kind=metis_index_kind), intent(in) :: options(*)
      integer(kind=metis_index_kind), intent(inout) :: objval(*)
      integer(kind=metis_index_kind), intent(inout) :: part(*)
 
    end function metis_partgraphkway
 
  end interface
#endif
 
contains
 
#ifdef DO_GRAPHLIB
  subroutine metis_setdefaultoptions_wrapper(options)
 
    integer(kind=metis_index_kind), intent(in) :: options(:)
    integer :: result
 
    result = metis_setdefaultoptions(options)
    if (result /= 1) then
      write(*, *) "error calling METIS_SetDefaultOptions"
      stop
    end if
 
  end subroutine metis_setdefaultoptions_wrapper
 
  subroutine metis_partgraphkway_wrapper(nvtxs, ncon, xadj, adjncy, vwgt, &
       vsize, adjwgt, nparts, tpwgts, ubvec, options, objval, part)
 
    integer, intent(in) :: nvtxs
    integer, intent(in) :: ncon
    integer, intent(in) :: xadj(:)
    integer, intent(in) :: adjncy(:)
    integer, pointer, intent(in) :: vwgt(:)
    integer, pointer, intent(in) :: vsize(:)
    integer, pointer, intent(in) :: adjwgt(:)
    integer, intent(in) :: nparts
    double precision, pointer, intent(in) :: tpwgts(:)
    double precision, pointer, intent(in) :: ubvec(:)
    integer(kind=metis_index_kind), intent(in) :: options(:)
    integer, intent(inout) :: objval
    integer, intent(inout) :: part(:)
 
    integer(kind=metis_index_kind) :: nvtxs_metis(1)
    integer(kind=metis_index_kind) :: ncon_metis(1)
    integer(kind=metis_index_kind), allocatable :: xadj_metis(:)
    integer(kind=metis_index_kind), allocatable :: adjncy_metis(:)
    integer(kind=metis_index_kind), pointer :: vwgt_metis(:) => null()
    integer(kind=metis_index_kind), pointer :: vsize_metis(:) => null()
    integer(kind=metis_index_kind), pointer :: adjwgt_metis(:) => null()
    integer(kind=metis_index_kind) :: nparts_metis(1)
    real(kind=metis_real_kind), pointer :: tpwgts_metis(:) => null()
    real(kind=metis_real_kind), pointer :: ubvec_metis(:) => null()
    integer(kind=metis_index_kind) :: objval_metis(1)
    integer(kind=metis_index_kind), allocatable :: part_metis(:)
 
    integer :: result
 
    nvtxs_metis(1) = nvtxs
    ncon_metis(1) = ncon
 
    allocate(xadj_metis(size(xadj)))
    xadj_metis = xadj
 
    allocate(adjncy_metis(size(adjncy)))
    adjncy_metis = adjncy
 
    if (associated(vwgt)) then
      allocate(vwgt_metis(size(vwgt)))
      vwgt_metis = vwgt
    end if
 
    if (associated(vsize)) then
      allocate(vsize_metis(size(vsize)))
      vsize_metis = vsize
    end if
 
    if (associated(adjwgt)) then
      allocate(adjwgt_metis(size(adjwgt)))
      adjwgt_metis = adjwgt
    end if
 
    nparts_metis(1) = nparts
 
    if (associated(tpwgts)) then
      allocate(tpwgts_metis(size(tpwgts)))
      tpwgts_metis = tpwgts
    end if
 
    if (associated(ubvec)) then
      allocate(ubvec_metis(size(ubvec)))
      ubvec_metis = ubvec
    end if
 
    objval_metis(1) = objval
    part_metis = part
 
    result = metis_partgraphkway(nvtxs_metis, ncon_metis, xadj_metis, adjncy_metis, vwgt_metis, vsize_metis, adjwgt_metis, &
         nparts_metis, tpwgts_metis, ubvec_metis, options, objval_metis, part_metis)
    if (result /= 1) then
      write(*, *) "error calling METIS_PartGraphKway"
      stop
    end if
 
    if (associated(vwgt_metis)) then
      deallocate(vwgt_metis)
    end if
 
    if (associated(vsize_metis)) then
      deallocate(vsize_metis)
    end if
 
    if (associated(adjwgt_metis)) then
      deallocate(adjwgt_metis)
    end if
 
    if (associated(tpwgts_metis)) then
      deallocate(tpwgts_metis)
    end if
 
    if (associated(ubvec_metis)) then
      deallocate(ubvec_metis)
    end if
 
    objval = objval_metis(1)
    part = part_metis
 
  end subroutine metis_partgraphkway_wrapper
#endif
 
  subroutine prg_metispartition(gp, ngroups, nnodes, xadj, adjncy, nparts, part, core_count, CH_count, Halo_count, sumCubes, &
       maxCH, smooth_maxCH, pnorm)
 
    implicit none
 
    type (graph_partitioning_t), intent(inout) :: gp
 
    integer(kind=metis_index_kind), allocatable :: options(:)  ! options for metis (64-bit wide)
    integer, allocatable, intent(inout)         :: xadj(:), adjncy(:), part(:)
    integer, intent(inout)                      :: nparts
    integer                                     :: ncon, objval
    integer                                     :: i, j
    integer, target                             :: dummy_vwgt, dummy_vsize, dummy_adjwgt
    real(8), target                             :: dummy_tpwgts, dummy_ubvec
    real(dp), intent (inout)                    :: sumcubes, maxch, smooth_maxch, pnorm
    integer, intent (in)                        :: ngroups, nnodes
    integer, allocatable, intent(inout)         :: ch_count(:), core_count(:)
    integer, allocatable, intent(inout)         :: halo_count(:,:)
    integer, allocatable                        :: copy_core_count(:)
    integer, pointer                            :: vwgt(:) => null(), vsize(:) => null(), adjwgt(:) => null()
    ! type(c_ptr)                                 :: vwgt, vsize, adjwgt
    ! type(c_ptr)                                 :: tpwgts, ubvec
    real(8), pointer                            :: tpwgts(:) => null(), ubvec(:) => null()
    character(len=100)                          :: pname
 
    allocate(options(40))
    allocate(copy_core_count(nparts))
 
    write(pname, '("metisParts")')
 
    options=0
#ifdef DO_GRAPHLIB
    call metis_setdefaultoptions_wrapper(options)
#endif
 
    ncon        = 1
    objval      = 1
 
    options(1)  = 1  !METIS_PTYPE_KWAY
    options(2)  = 0  !METIS_OBJTYPE_CUT
    options(9)  = 10 !METIS_OPTION_SEED
    options(18) = 1  !Fortran-style numbering is assumed that starts from 1
 
    halo_count  = 0
    ch_count    = 0
    core_count  = 0
 
    if (printrank() .eq. 1) then
      write(*,*) "prg_metisPartition_test start ..."
    endif
    ! prg_initialize gp
    call prg_initgraphpartitioning(gp, pname, nparts, ngroups, nnodes)
 
    !write(*,*) "The number of nodes in the graph is:", gp%totalNodes, &
    !     gp%totalNodes2, ncon, nparts, objval
 
#ifdef DO_GRAPHLIB
    call metis_partgraphkway_wrapper(gp%totalNodes, ncon, xadj, adjncy, vwgt, &
         vsize, adjwgt, nparts, tpwgts, ubvec, options, objval, part)
#endif
 
    call prg_costpartition(gp, xadj, adjncy, part, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
 
    !prg_initialize and fill up subgraph structure
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    do i = 1, nparts
      gp%nnodesInPartAll(i) = core_count(i)
      copy_core_count(i) = core_count(i)
      call prg_initsubgraph(gp%sgraph(i), i, nnodes)
      allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
      gp%nnodesInPart(i) = core_count(i)
    enddo
 
    !Assign node ids to sgraph
    do i = 1, gp%totalNodes
      copy_core_count(part(i)) = copy_core_count(part(i)) - 1
      ! --> core_count((part(i))) - copy_core_count(part(i)) <--- = 1,2,3, ..., core_count( partnumber ) with respect to i
 
      ! NOTE: nodes in gp%sgraph()%nodeInPart() are currently 0 based!
      gp%sgraph(part(i))%nodeInPart(core_count((part(i))) - copy_core_count(part(i)) ) = i -1
    end do
    do i = 1, nparts
      do j = 1, core_count(i)
        if( part( gp%sgraph(i)%nodeInPart(j)+1 ) /= i) then
          write(*,*) "ERROR: subgraph struc incorrect!!", "node=",gp%sgraph(i)%nodeInPart(j)+1 , &
               "part=",i, "actual_part=", part(gp%sgraph(i)%nodeInPart(j)+1 )
          stop
        end if
 
      end do
    end do
 
  end subroutine prg_metispartition
!
  subroutine prg_costpartition(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm)
 
    type (graph_partitioning_t), intent(inout)  :: gp
    integer, allocatable, intent(inout)         :: xadj(:), adjncy(:)
    integer, allocatable, intent(in)            :: partnumber(:)
    integer, allocatable, intent(inout)         :: core_count(:)
    integer                                     :: totalparts, totalnodes,  i, j, neighbor
    real(dp), intent (inout)                    :: sumcubes, maxch, smooth_maxch, pnorm
    integer, allocatable, intent(inout)         ::  ch_count(:)
    integer, allocatable, intent(inout)         :: halo_count(:,:)
    real(dp)                                    :: temp
 
    maxch        = 0
    smooth_maxch = 0
    sumcubes     = 0
    totalparts   = gp%totalParts
    totalnodes   = gp%totalNodes
 
    halo_count = 0
    ch_count   = 0
    core_count = 0
 
    do i = 1, totalnodes
      ch_count(partnumber(i)) = ch_count(partnumber(i)) + 1 !core count
      core_count(partnumber(i)) = core_count(partnumber(i)) + 1 !core count
      do j = xadj(i), xadj(i + 1) - 1
        neighbor = adjncy(j)
        if (partnumber(i) /= partnumber(neighbor)) then
          if (halo_count(partnumber(i) ,neighbor) == 0) then
            ch_count(partnumber(i)) = ch_count(partnumber(i)) + 1 !halo count
            halo_count(partnumber(i), neighbor) = 1
          else
            halo_count(partnumber(i), neighbor) = halo_count(partnumber(i), neighbor) + 1
          end if
        end if
      end do
    end do
 
    do i = 1, totalparts
      if (core_count(i) <= 1) then
        print *, "core count <= 1 for partition "//to_string(i)//"!"
        stop
      end if
      temp = real(ch_count(i), dp)
      sumcubes = sumcubes+  temp*temp*temp
      smooth_maxch = smooth_maxch + temp**int(pnorm)
      if (ch_count(i) > maxch) then
        maxch = ch_count(i)
      end if
    end do
    smooth_maxch = smooth_maxch**(1/pnorm)
 
  end subroutine prg_costpartition
 
  subroutine update_prg_costpartition(gp, xadj, adjncy, partNumber,core_count, CH_count, Halo_count, sumCubes, maxCH,smooth_maxCH, pnorm, node, new_part)
 
    type (graph_partitioning_t), intent(inout)  :: gp
    integer, allocatable, intent(inout)         :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)         :: partnumber(:), core_count(:)
    integer                                     :: totalparts, totalnodes,  i, j, neighbor
    real(dp), intent (inout)                    :: sumcubes, maxch, smooth_maxch, pnorm
    integer,  intent (in)                       :: node, new_part
    integer, allocatable, intent(inout)         :: ch_count(:)
    integer, allocatable, intent(inout)         :: halo_count(:,:)
    integer                                     :: old_part
    real(dp)                                    :: temp
 
    totalparts = gp%totalParts
    totalnodes = gp%totalNodes
    old_part   = partnumber(node)
 
    if (old_part /= new_part) then
      core_count(new_part) = core_count(new_part)  + 1
      core_count(old_part) = core_count(old_part)  - 1
      ch_count(old_part) = ch_count(old_part) - 1 !core--
      ch_count(new_part) = ch_count(new_part) + 1 !core--
      do i=xadj(node), xadj(node+1) -1
        neighbor = adjncy(i)
        if (node /= neighbor) then
          if(partnumber(neighbor) == old_part) then !case 1
            halo_count(old_part, node) = halo_count(old_part, node) + 1
            if (halo_count(old_part, node) == 1) then
              ch_count(old_part) = ch_count(old_part) + 1 !halo++
            end if
            halo_count(new_part, neighbor) = halo_count(new_part, neighbor) + 1
            if (halo_count(new_part, neighbor) == 1) then
              ch_count(new_part) = ch_count(new_part) + 1 !halo++
            end if
          else if (partnumber(neighbor) == new_part) then !case 2
            halo_count(old_part, neighbor) = halo_count(old_part, neighbor) - 1
            if (halo_count(old_part, neighbor) == 0) then
              ch_count(old_part) = ch_count(old_part) - 1 !halo--
            else if (halo_count(old_part, neighbor) < 0) then
              write(*,*) "ERROR: Halo_count value cannot be negative, case 2i"
              write(*,*) "input matrix should be perfectly symmetric"
              stop
            end if
            halo_count(new_part, node) = halo_count(new_part, node) - 1
            if (halo_count(new_part, node) == 0) then
              ch_count(new_part) = ch_count(new_part) - 1 !halo--, halo has become a core
            else if (halo_count(new_part, node) < 0) then
              write(*,*) "ERROR: Halo_count value cannot be negative, case 2ii"
              write(*,*) "input matrix should be perfectly symmetric"
              stop
            end if
          else !case 3
            halo_count(old_part, neighbor) = halo_count(old_part, neighbor) - 1
            if (halo_count(old_part, neighbor) == 0) then
              ch_count(old_part) = ch_count(old_part) - 1 !halo--
            else if (halo_count(old_part, neighbor) < 0) then
              write(*,*) "ERROR: Halo_count value cannot be negative, case 3"
              write(*,*) "input matrix should be perfectly symmetric"
              stop
            end if
            halo_count(new_part, neighbor) = halo_count(new_part, neighbor) + 1
            if (halo_count(new_part, neighbor) == 1) then
              ch_count(new_part) = ch_count(new_part) + 1 !halo++
            end if
          end if
        end if
      end do
      partnumber(node) = new_part
      sumcubes     = 0
      maxch        = 0
      smooth_maxch = 0
      do i=1, totalparts
        temp = real(ch_count(i), dp)
        sumcubes = sumcubes+  temp*temp*temp
        smooth_maxch = smooth_maxch + temp**int(pnorm)
        if (ch_count(i) > maxch) then
          maxch = real(ch_count(i),dp)
        end if
      end do
      smooth_maxch = smooth_maxch**(1/pnorm)
    end if
 
  end subroutine update_prg_costpartition
 
  subroutine prg_accept_prob(it, prg_delta, r)
    integer, intent(in)     :: it
    real(dp), intent(in)    :: prg_delta
    real, intent(inout)     :: r
    real                    :: temp
    temp = 1.0/it
    r = 1
    if (prg_delta > 0) then
      r = exp(-(prg_delta/10.0)/temp)
    end if
 
  end subroutine prg_accept_prob
 
  subroutine prg_costindex(cost, sumCubes, maxCH, smooth_maxCH, obj_fun)
    real(dp), intent (inout) :: cost, maxCH,smooth_maxCH, sumCubes
    integer, intent (inout) ::  obj_fun
 
    cost = -1
 
    if (obj_fun .eq. 0) then
      cost = sumcubes
    else if (obj_fun .eq. 1) then
      cost = maxch
    else if (obj_fun .eq. 2) then
      cost = smooth_maxch
    end if
 
  end subroutine prg_costindex
 
  subroutine prg_rand_node(gp,node, seed)
    type (graph_partitioning_t), intent(inout) :: gp
    integer, intent(inout) :: node, seed
    integer ::  totalNodes,  i, ssize
    integer, allocatable :: seedin(:)
    real :: u
    !call srand(seed)
    call random_seed()
    call random_seed(size=ssize)
    allocate(seedin(ssize))
    seedin = seed
    call random_seed(put=seedin)
 
    call random_number(u)
    node = floor(gp%totalNodes*u) + 1
 
  end subroutine prg_rand_node
 
  subroutine prg_simannealing(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH,smooth_maxCH,pnorm, niter, seed)
 
    type (graph_partitioning_t), intent(inout)  :: gp
    integer, allocatable, intent(inout)         :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)         :: partnumber(:), core_count(:)
    integer                                     :: totalparts, totalnodes,  it, i,j,k, neighbor,  node, part_backup
    integer                                     :: totalnodes2
    real(dp)                                    :: cost, prev_cost, prg_delta, prev_maxch
    real(dp), intent (inout)                    :: sumcubes, maxch, smooth_maxch, pnorm
    integer,  intent (in)                       :: niter
    integer, intent(inout)                      :: seed
    integer, allocatable, intent(inout)         :: ch_count(:)
    integer, allocatable, intent(inout)         :: halo_count(:,:)
    integer, allocatable                        :: copy_core_count(:), empty_parts(:)
    integer                                     :: obj_fun = 2, min_ch_part, no_empty_parts,  new_part
    real                                        :: r, u
    character(len=100)                          :: pname
    integer :: ssize
    integer, allocatable :: seedin(:)
 
    totalparts = gp%totalParts
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
    allocate(copy_core_count(totalparts))
    allocate(empty_parts(totalparts))
    !call srand(seed)
    call random_seed()
    call random_seed(size=ssize)
    allocate(seedin(ssize))
    seedin = seed
    call random_seed(put=seedin)
 
    write(*,*) "SA called..."
    if (totalnodes .lt. totalparts) then
      write(*,*) "ERROR: Number of parts cannot be greater than number of nodes."
      stop
    end if
 
    call prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
 
    call prg_costindex(cost, sumcubes, maxch,smooth_maxch, obj_fun)
    prev_cost = cost
 
    do it=1, niter
      call prg_rand_node(gp, node, seed)
      min_ch_part = 1
      do k =1, totalparts
        if (ch_count(min_ch_part) .gt. ch_count(k)) then
          min_ch_part = k
        end if
      end do
 
      if (ch_count( partnumber(node) ) .eq. maxch) then
        do j = xadj(node), xadj(node+1)-1
          neighbor = adjncy(j)
          part_backup = partnumber(neighbor)
          call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, &
               smooth_maxch, pnorm, neighbor, min_ch_part)
          call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
          prg_delta = cost - prev_cost
          call prg_accept_prob(it, prg_delta, r)
          call random_number(u)
          if (u > r) then ! reject
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup)
          else
            prev_cost = cost
          end if
        end do
      else
        if (ch_count( min_ch_part) .eq. 0) then
          do j = xadj(node), xadj(node+1)-1
            neighbor = adjncy(j)
            part_backup = partnumber(neighbor)
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, min_ch_part)
            call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
            prg_delta = cost - prev_cost
            call prg_accept_prob(it, prg_delta, r)
            call random_number(u)
            if (u > r) then ! reject
              call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm,  neighbor, part_backup)
            else
              prev_cost = cost
            end if
          end do
        else
          do j = xadj(node), xadj(node+1)-1
            neighbor = adjncy(j)
            part_backup = partnumber(neighbor)
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, partnumber(node))
            call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
            prg_delta = cost - prev_cost
            call prg_accept_prob(it, prg_delta, r)
            call random_number(u)
            if (u > r) then ! reject
              call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup)
            else
              prev_cost = cost
            end if
          end do
        end if
      end if
    end do
 
    no_empty_parts = 0
    do i=1,totalparts
      if (ch_count(i) .eq. 0) then
        no_empty_parts = no_empty_parts + 1
        empty_parts(no_empty_parts) = i
      end if
    end do
    prev_maxch = maxch
    do node=1,totalnodes
      if (no_empty_parts .le. 0) then
        exit
      end if
 
      if (ch_count(partnumber(node)) .eq. maxch) then
        new_part = empty_parts(no_empty_parts)
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, node,new_part )
        prev_maxch = maxch
        no_empty_parts = no_empty_parts - 1
 
        do j = xadj(node), xadj(node+1)-1
          neighbor = adjncy(j)
          if (ch_count(partnumber(neighbor)) .eq. maxch) then
            part_backup = partnumber(neighbor)
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, neighbor, new_part)
            call prg_costindex(cost, sumcubes, maxch,smooth_maxch, obj_fun)
            if (maxch .ge. prev_maxch) then ! reject
              call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup )
              call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
            else
              prev_cost = cost
              prev_maxch = maxch
            end if
          end if
        end do
      end if
    end do
 
    write(*,*) "Cost of meTIS+SA", sumcubes, maxch, smooth_maxch
 
 
 
    !prg_initialize and fill up subgraph structure
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    call prg_destroygraphpartitioning(gp)
    write(pname, '("SAParts")')
    call prg_initgraphpartitioning(gp, pname, totalparts, totalnodes, totalnodes2)
    do i = 1, totalparts
      gp%nnodesInPartAll(i) = core_count(i)
      copy_core_count(i) = core_count(i)
      call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
      allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
      gp%nnodesInPart(i) = core_count(i)
    enddo
 
    !Assign node ids to sgraph
    do i=1, gp%totalNodes
      copy_core_count(partnumber(i)) =copy_core_count(partnumber(i)) - 1
      ! core_count((part(i))) - copy_core_count(part(i)) = node postion in array
      gp%sgraph(partnumber(i))%nodeInPart(core_count((partnumber(i))) - copy_core_count(partnumber(i)) ) = i -1 !NOTE: nodes in gp%sgraph()%nodeInPart() are currently 0 based!
    end do
 
 
    do i = 1, totalparts
      do j = 1, core_count(i)
        if( partnumber( gp%sgraph(i)%nodeInPart(j)+1 ) /= i) then
          write(*,*) "ERROR: subgraph struc incorrect!!", "node=",gp%sgraph(i)%nodeInPart(j)+1 , "part=",i, "actual_part=", partnumber(gp%sgraph(i)%nodeInPart(j)+1 )
          stop
        end if
      end do
    end do
    do i=1, totalparts
      write(*,*) "part=",i, "C=", core_count(i), "CH=", ch_count(i)
      if (ch_count(i) .eq. 0) then
        write(*,*) "ERROR: SA produced an empty part"
        stop
      end if
    end do
  end subroutine prg_simannealing
 
 
 
 
  subroutine prg_kernlin(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm, nconverg, seed)
 
    type (graph_partitioning_t), intent(inout)    :: gp
    integer, allocatable, intent(inout)           :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)           :: partnumber(:), core_count(:)
    integer                                       :: totalparts, totalnodes, i, iit, j,k, neighbor,  node, part_backup, h, node2
    real(dp), intent(inout)                       :: sumcubes, maxch, smooth_maxch, pnorm
    real(dp)                                      :: cost, prev_cost, prev_iteration_cost, prev_maxch, minch
    integer, intent(inout)                        :: seed
    integer, intent(in)                           :: nconverg
    integer, allocatable, intent(inout)           :: ch_count(:)
    integer, allocatable, intent(inout)           :: halo_count(:,:)
    integer, allocatable                          :: copy_core_count(:)
    integer                                       :: obj_fun = 2, counter, min_part, max_climb = 1, climb_counter, temp, convg_counter, converge, no_locked_nodes, empty_counter, backup, best_node, best_part,no_empty_parts, new_part
    integer                                       :: totalnodes2
    real                                          :: r, u
    integer, allocatable                          :: vertex_locked(:), hedge_span(:), node_backup(:), node_part_backup(:), nodes(:), empty_parts(:)
    character(len=100)                            :: pname
 
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
    totalparts = gp%totalParts
 
    allocate(vertex_locked(totalnodes))
    allocate(hedge_span(totalparts))
    allocate(node_backup(totalnodes))
    allocate(node_part_backup(totalnodes))
    allocate(nodes(totalnodes))
    allocate(copy_core_count(totalparts))
    allocate(empty_parts(totalparts))
 
    vertex_locked   = 0
    hedge_span      = 0
    counter         = 0
    climb_counter   = 1
    converge        = 0
    convg_counter   = 0
    iit             = 0
    core_count      = 0
    ch_count        = 0
    halo_count      = 0
    no_locked_nodes = 0
 
    do  i = 1, totalnodes
      nodes(i) = i
    end do
 
    call prg_rand_shuffle(nodes, seed)
 
 
    call prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
 
    call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
    prev_cost           = cost
    prev_iteration_cost = cost
    prev_maxch          = maxch
 
 
 
    do while ( converge .eq. 0 .and. iit .lt. nconverg)
      iit = iit + 1
      vertex_locked=0 ! Free vertices
      no_locked_nodes = 0
      call prg_rand_shuffle(nodes, seed)
 
      do i=1, gp%totalNodes
        h = nodes(i) !h represents a hyperedge
 
        minch = totalnodes + 1
        do j = 1, totalparts
          if (ch_count(j) .lt. minch) then
            min_part = j
            minch = ch_count(j)
          end if
        end do
        if (min_part .eq. -1) then
          min_part = partnumber(h) !hyperedge h contains node h
          do j = xadj(h), xadj(h+1)-1
            node = adjncy(j)
            if (hedge_span( partnumber(node))==0) then
              counter = counter + 1
              hedge_span( partnumber(node)) = ch_count(partnumber(node))
              if (ch_count(partnumber(node)) .le. ch_count(min_part)) then
                min_part = partnumber(node)
              end if
            end if
            if (counter == totalparts) then
              counter = 0
              exit
            end if
          end do
        end if
        hedge_span    = 0
 
        if(h .eq. 0) then
          write(*,*) "error h =0"
          stop
        end if
        climb_counter = 1
        do j = xadj(h), xadj(h+1)-1
          node = adjncy(j)
 
          if (vertex_locked(node) .eq. 0  ) then
            part_backup = partnumber(node)
            node_part_backup(climb_counter) = part_backup
            node_backup(climb_counter) = node
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, node, min_part)
            call prg_costindex(cost, sumcubes, maxch,smooth_maxch, obj_fun)
            if (cost .le. prev_cost) then !accept
              prev_cost = cost
              no_locked_nodes = no_locked_nodes + climb_counter
              !write(*,*) maxCH
 
 
              temp = climb_counter
              do k=1, temp
                vertex_locked(node_backup(climb_counter)) = 1
                climb_counter = climb_counter - 1
              end do
 
              climb_counter = 1
              node_backup = -1 !for debug purposes
              node_part_backup = -1 !for debug purposes
 
            else
              if (climb_counter .lt. max_climb) then !climb
                climb_counter = climb_counter + 1
              else !undo climb
                do k=1, max_climb
                  call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, node_backup(climb_counter), node_part_backup(climb_counter))
                  call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
 
                  climb_counter = climb_counter - 1
                end do
                if (prev_cost .ne. cost) then
                  write(*,*) "ERROR: There was an error in undo step 2", node, cost, prev_cost, j
                  stop
                end if
                climb_counter = 1 !reset
              end if
            end if
          end if
 
          !            end if
        end do
        !end of hyperedge, undo any climbs
        if(prev_cost .ne. cost) then
          temp = climb_counter
 
          !if climb_counter = 1, we cannot have prev_cost /= cost
          do k=1, temp-1
            climb_counter = climb_counter - 1
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, node_backup(climb_counter), node_part_backup(climb_counter))
            call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
 
          end do
          if (prev_cost .ne. cost) then
            write(*,*) "ERROR: Undo after hyperedge"
            stop
          end if
        end if
 
        if (no_locked_nodes .eq. gp%totalNodes) then
          exit
        end if
 
 
        empty_counter = 0
        do j=1,totalparts
          if (core_count(j) .eq. 0) then
            empty_counter = empty_counter +1
            empty_parts(empty_counter) = j
          end if
        end do
        if (empty_counter .gt. 0) then
          do j= 1,totalnodes
            if (ch_count(partnumber(j) ) .eq. maxch) then
              do k = xadj(j), xadj(j+1)-1
                node2 = adjncy(k)
                backup = partnumber(j)
                if (partnumber(node2) .eq. partnumber(j) ) then
                  call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, node2,empty_parts(empty_counter) )
                  call prg_costindex(cost, sumcubes, maxch,smooth_maxch,  obj_fun)
                  prev_cost = cost
                  !prev_maxCH = maxCH
                  if (prev_maxch .lt. maxch) then !undo
                    call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, node2,backup )
                    call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
                    prev_cost = cost
                    prev_maxch = maxch
                  end if
                end if
              end do
              empty_counter = empty_counter - 1
              if (empty_counter .eq. 0) then
                exit
              end if
            end if
          end do
        end if
      end do !End KL iterations
 
      if (prev_iteration_cost .eq. cost) then
        convg_counter = convg_counter + 1
        if (convg_counter .eq. nconverg) then
          converge =1
        end if
      else
        prev_iteration_cost = cost
        convg_counter = 0
      end if
 
    end do !while loop
 
    no_empty_parts = 0
 
    do i=1,totalparts
      if (ch_count(i) .eq. 0) then
        no_empty_parts = no_empty_parts + 1
        empty_parts(no_empty_parts) = i
      end if
    end do
    prev_maxch = maxch
    do node=1,totalnodes
      if (no_empty_parts .le. 0) then
        exit
      end if
 
      if (ch_count(partnumber(node)) .eq. maxch .and. core_count(partnumber(node)) .ne. 1 ) then
        new_part = empty_parts(no_empty_parts)
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, node,new_part )
        prev_maxch = maxch
        no_empty_parts = no_empty_parts - 1
 
        do j = xadj(node), xadj(node+1)-1
          neighbor = adjncy(j)
          if (ch_count(partnumber(neighbor)) .eq. maxch) then
            part_backup = partnumber(neighbor)
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, neighbor, new_part)
            call prg_costindex(cost, sumcubes, maxch,smooth_maxch, obj_fun)
            if (maxch .ge. prev_maxch) then ! reject
              call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup )
              call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
 
            else
              prev_cost = cost
              prev_maxch = maxch
            end if
          end if
        end do
      end if
    end do
 
 
 
 
 
 
    write(*,*) "Cost of KL:", cost, maxch, "No iterations:", iit
 
    !prg_initialize and fill up subgraph structure
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    call prg_destroygraphpartitioning(gp)
    write(pname, '("KLParts")')
    call prg_initgraphpartitioning(gp, pname, totalparts, totalnodes, totalnodes2)
 
    do i = 1, totalparts
      gp%nnodesInPartAll(i) = core_count(i)
      copy_core_count(i) = core_count(i)
      call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
      allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
      gp%nnodesInPart(i) = core_count(i)
    enddo
 
    do i=1, gp%totalNodes
      copy_core_count(partnumber(i)) =copy_core_count(partnumber(i)) - 1
      !!>core_count((part(i))) - copy_core_count(part(i)) = node postion in array
      gp%sgraph(partnumber(i))%nodeInPart(core_count((partnumber(i))) - copy_core_count(partnumber(i)) ) = i -1 !NOTE: nodes in gp%sgraph()%nodeInPart() are currently 0 based!
    end do
 
 
    do i=1, totalparts
      write(*,*) "part=",i, "C=", core_count(i), "CH=", ch_count(i)
      if (ch_count(i) .eq. 0) then
        write(*,*) "ERROR: KL produced an empty part"
        stop
      end if
    end do
    deallocate(vertex_locked)
    deallocate(hedge_span)
    deallocate(node_backup)
    deallocate(node_part_backup)
    deallocate(nodes)
    deallocate(copy_core_count)
 
  end subroutine prg_kernlin
 
 
 
  subroutine prg_update_gp(gp, partNumber, core_count)
    type (graph_partitioning_t), intent(inout)    :: gp
    integer, allocatable, intent(inout)           :: partnumber(:), core_count(:)
    integer                                       :: totalparts, totalnodes, i
    integer                                       :: totalnodes2
    integer, allocatable                          :: copy_core_count(:)
    character(len=100)                            :: pname
 
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
    totalparts = gp%totalParts
    allocate(copy_core_count(totalparts))
    !prg_initialize and fill up subgraph structure
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    call prg_destroygraphpartitioning(gp)
    write(pname, '("Parts")')
    call prg_initgraphpartitioning(gp, pname, totalparts, totalnodes, totalnodes2)
 
    do i = 1, totalparts
      gp%nnodesInPartAll(i) = core_count(i)
      copy_core_count(i) = core_count(i)
      call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
      allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
      gp%nnodesInPart(i) = core_count(i)
    enddo
 
    do i=1, gp%totalNodes
      copy_core_count(partnumber(i)) =copy_core_count(partnumber(i)) - 1
      !!>core_count((part(i))) - copy_core_count(part(i)) = node postion in array
      gp%sgraph(partnumber(i))%nodeInPart(core_count((partnumber(i))) - copy_core_count(partnumber(i)) ) = i -1 !NOTE: nodes in gp%sgraph()%nodeInPart() are currently 0 based!
    end do
 
 
    deallocate(copy_core_count)
  end subroutine prg_update_gp
 
 
  subroutine prg_rand_shuffle(array, seed)
 
    integer, intent(inout)  :: array(:), seed
    integer                 :: i, randpos, temp, ssize
    integer, allocatable :: seedin(:)
    real                    :: r
 
    !call srand(seed)
    call random_seed()
    call random_seed(size=ssize)
    allocate(seedin(ssize))
    seedin = seed
    call random_seed(put=seedin)
 
 
    do i = size(array), 2, -1
      call random_number(r)
      randpos = int(r * i) + 1
      temp = array(randpos)
      array(randpos) = array(i)
      array(i) = temp
    end do
 
  end subroutine prg_rand_shuffle
 
 
  subroutine prg_check_arrays(gp, core_count, CH_count, Halo_count)
    type (graph_partitioning_t), intent(inout) :: gp
    integer, allocatable, intent(inout)  :: core_count(:)
    integer, allocatable, intent(inout) ::  ch_count(:)
    integer, allocatable, intent(inout)  :: halo_count(:,:)
    integer                              :: i,j, check
 
    do i=1,gp%totalParts
      check = core_count(i)
      do j=1, gp%totalNodes
        if (halo_count(i,j) >0) then
          check = check + 1
        end if
      end do
      if (check /= ch_count(i)) then
 
        write(*,*) "ERROR: Halo_count is incorrect!"
        write(*,*) "check=", check, "CH_count(i)", ch_count(i)
        stop
      end if
    end do
 
    write(*,*) "prg_check_arrays PASSED!"
  end subroutine prg_check_arrays
 
  subroutine prg_kernlin_queue(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm)
 
    type (graph_partitioning_t), intent(inout)    :: gp
    integer, allocatable, intent(inout)           :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)           :: partnumber(:), core_count(:)
    integer                                       :: totalparts, totalnodes, i, iit, j,k, neighbor,  node, part_backup, h, node2
    real(dp), intent(inout)                       :: sumcubes, maxch, smooth_maxch, pnorm
    real(dp)                                      :: cost, prev_cost, prev_iteration_cost, prev_maxch, minch, best_obj_val, current_cost
    integer, allocatable, intent(inout)           :: ch_count(:)
    integer, allocatable, intent(inout)           :: halo_count(:,:)
    integer, allocatable                          :: copy_core_count(:)
    integer                                       :: backup, best_node, best_part
    real                                          :: r, u
    integer, allocatable                          :: vertex_locked(:), hedge_span(:), node_backup(:), node_part_backup(:), nodes(:), emptyparts(:)
    character(len=100)                            :: pname
 
 
 
    do i =1, 5
      call prg_find_best_move(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, best_node, best_part )
      if(best_node .eq. 0) then
        write(*,*) "error: node is 0"
        stop
      end if
      call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, best_node, best_part)
 
 
    end do
    write(*,*) "kl finished", gp%totalParts
    do i=1,gp%totalParts
      write(*,*) "part=", i, core_count(i), "ch=", ch_count(i)
    end do
 
  end subroutine prg_kernlin_queue
 
  subroutine prg_find_best_move(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm, best_node, best_part )
 
    type (graph_partitioning_t), intent(inout)    :: gp
    integer, allocatable, intent(inout)           :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)           :: partNumber(:), core_count(:)
    integer                                       :: totalParts, totalNodes, i, iit, j,k, neighbor,  node, part_backup, h, node2
    integer                                       :: totalNodes2
    real(dp), intent(inout)                       :: sumCubes, maxCH, smooth_maxCH, pnorm
    real(dp)                                      :: cost, prev_cost, prev_iteration_cost, prev_maxCH, minCH, best_obj_val, current_cost
    integer, intent(inout)                        :: best_node, best_part
    integer, allocatable, intent(inout)           :: CH_count(:)
    integer, allocatable, intent(inout)           :: Halo_count(:,:)
    integer, allocatable                          :: copy_core_count(:)
    integer                                       :: obj_fun = 2,  backup
    real                                          :: r, u
    integer, allocatable                          :: vertex_locked(:), hedge_span(:), node_backup(:), node_part_backup(:), nodes(:), emptyParts(:)
    character(len=100)                            :: pname
    !!  iterate through hyperedges
    !!  find min part and empty part
    !!  move vertices into part incident to h with smalles CH_count or empty part if it exists
 
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
    totalparts = gp%totalParts
 
    best_node = 0
    best_part = 0
    call prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
    call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
    best_obj_val = cost
 
    do i=1, totalnodes
      do j=1,totalparts
        backup = partnumber(i)
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, i, j)
        call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
 
        if (cost .le. best_obj_val) then
          best_node = i
          best_part = j
          best_obj_val = cost
        end if
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, i, backup)
      end do
    end do
 
  end subroutine prg_find_best_move
 
  subroutine prg_kernlin2(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm)
 
    type (graph_partitioning_t), intent(inout)    :: gp
    integer, allocatable, intent(inout)           :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)           :: partnumber(:), core_count(:)
    integer                                       :: totalparts, totalnodes, i, it, iit, j,k, neighbor,  node, part_backup
    integer                                       :: totalnodes2
    real(dp), intent(inout)                       :: sumcubes, maxch, smooth_maxch, pnorm
    real(dp)                                      :: cost, prev_cost, prev_maxch, minch
    integer, allocatable, intent(inout)           :: ch_count(:)
    real                                          :: r
    integer, allocatable, intent(inout)           :: halo_count(:,:)
    integer, allocatable                          :: copy_core_count(:),  empty_parts(:)
    integer                                       :: obj_fun = 2, largest_hedge = -1, search_part, min_ch_part, new_part, no_empty_parts, seed =1
    character(len=100)                            :: pname
    !!  iterate through hyperedges
    !!  find min part and empty part
    !!  move vertices into part incident to h with smalles CH_count or empty part if it exists
 
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
    totalparts = gp%totalParts
 
    allocate(copy_core_count(totalparts))
    allocate(empty_parts(totalparts))
 
    call prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
    call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
    prev_cost = cost
    do iit = 1, 40
      do i = 1, totalparts
        call random_number(r)
        if (r .ge. 0.7) then
          it = i
          call prg_get_largest_hedge_in_part(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, it, largest_hedge)
        else
          call prg_rand_node(gp, largest_hedge, seed)
 
        end if
        min_ch_part = 1
        do k =1, totalparts
          if (ch_count(min_ch_part) .gt. ch_count(k)) then
            min_ch_part = k
          end if
        end do
 
        if (ch_count(partnumber(largest_hedge)) .eq. maxch) then
          new_part = min_ch_part
        else
          new_part = partnumber(largest_hedge)
        end if
 
        do j = xadj(largest_hedge), xadj(largest_hedge + 1)-1
          neighbor = adjncy(j)
          part_backup = partnumber(neighbor)
          call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, new_part)
          call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
          if (cost .gt. prev_cost) then ! reject
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup)
            call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
          else
            prev_cost = cost
          end if
        end do
      end do
    end do
 
    !call prg_Kernlin_queue(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm)
 
    no_empty_parts = 0
 
    do i=1,totalparts
      if (ch_count(i) .eq. 0) then
        no_empty_parts = no_empty_parts + 1
        empty_parts(no_empty_parts) = i
      end if
    end do
 
    prev_maxch = maxch
    do node=1,totalnodes
      if (no_empty_parts .le. 0) then
        exit
      end if
 
      if (ch_count(partnumber(node)) .eq. maxch .and. core_count(partnumber(node)) .ne. 1 ) then
        new_part = empty_parts(no_empty_parts)
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, node,new_part )
        prev_maxch = maxch
        no_empty_parts = no_empty_parts - 1
 
        do j = xadj(node), xadj(node+1)-1
          neighbor = adjncy(j)
          if (ch_count(partnumber(neighbor)) .eq. maxch) then
            part_backup = partnumber(neighbor)
            call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch,smooth_maxch, pnorm, neighbor, new_part)
            call prg_costindex(cost, sumcubes, maxch,smooth_maxch, obj_fun)
            if (maxch .ge. prev_maxch) then ! reject
              call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup )
              call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
 
            else
              prev_cost = cost
              prev_maxch = maxch
            end if
          end if
        end do
      else if (core_count(partnumber(node)) .ne. 1 ) then
        new_part = empty_parts(no_empty_parts)
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, node,new_part )
        prev_maxch = maxch
        no_empty_parts = no_empty_parts - 1
      end if
    end do
 
 
    !prg_initialize and fill up subgraph structure
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    call prg_destroygraphpartitioning(gp)
    write(pname, '("KLParts")')
    call prg_initgraphpartitioning(gp, pname, totalparts, totalnodes, totalnodes2)
 
    do i = 1, totalparts
      gp%nnodesInPartAll(i) = core_count(i)
      copy_core_count(i) = core_count(i)
      call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
      allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
      gp%nnodesInPart(i) = core_count(i)
    enddo
 
    do i=1, gp%totalNodes
      copy_core_count(partnumber(i)) =copy_core_count(partnumber(i)) - 1
      !!>core_count((part(i))) - copy_core_count(part(i)) = node postion in array
      gp%sgraph(partnumber(i))%nodeInPart(core_count((partnumber(i))) - copy_core_count(partnumber(i)) ) = i -1 !NOTE: nodes in gp%sgraph()%nodeInPart() are currently 0 based!
    end do
 
 
    do i=1, totalparts
      write(*,*) "part=",i, "C=", core_count(i), "CH=", ch_count(i)
      if (ch_count(i) .eq. 0) then
        write(*,*) "ERROR: KL produced an empty part"
        stop
      end if
    end do
 
  end subroutine prg_kernlin2
 
 
  subroutine prg_get_largest_hedge_in_part(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH, smooth_maxCH, pnorm, search_part, largest_Hedge)
    type (graph_partitioning_t), intent(inout)    :: gp
    integer, allocatable, intent(inout)           :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)           :: partNumber(:), core_count(:)
    integer                                       :: totalParts, totalNodes, i, iit, j,k, neighbor,  node
    integer                                       :: totalNodes2
    real(dp), intent(inout)                       :: sumCubes, maxCH, smooth_maxCH, pnorm
    integer, allocatable, intent(inout)           :: CH_count(:)
    integer, allocatable, intent(inout)           :: Halo_count(:,:)
    integer, allocatable                          :: copy_core_count(:)
    integer                                       :: obj_fun = 2, largest_hsize
    integer, intent(inout)                        :: search_part, largest_Hedge
    real                                          :: r, u
 
 
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
    totalparts = gp%totalParts
    largest_hsize = 0
    do i=1, totalnodes
      if (partnumber(i) .eq. search_part) then
        if ( xadj(i + 1) - xadj(i) .gt. largest_hsize) then
          largest_hsize = xadj(i + 1) - xadj(i)
          largest_hedge = i
        end if
      end if
    end do
 
  end subroutine prg_get_largest_hedge_in_part
 
 
  subroutine prg_simannealing_old(gp, xadj, adjncy, partNumber, core_count, CH_count, Halo_count, sumCubes, maxCH,smooth_maxCH,pnorm, niter, seed)
 
    type (graph_partitioning_t), intent(inout)  :: gp
    integer, allocatable, intent(inout)         :: xadj(:), adjncy(:)
    integer, allocatable, intent(inout)         :: partnumber(:), core_count(:)
    integer                                     :: totalparts, totalnodes,  it, i,j,k, neighbor,  node, part_backup, obj_fun=0
    integer                                     :: totalnodes2
    real(dp)                                    :: cost, prev_cost, prg_delta, prev_maxch
    real(dp), intent (inout)                    :: sumcubes, maxch, smooth_maxch, pnorm
    integer,  intent (in)                       :: niter
    integer, intent(inout)                      :: seed
    integer, allocatable, intent(inout)         :: ch_count(:)
    integer, allocatable, intent(inout)         :: halo_count(:,:)
    integer, allocatable                        :: copy_core_count(:), empty_parts(:)
    real                                        :: r, u
    character(len=100)                          :: pname
 
    totalparts = gp%totalParts
    totalnodes = gp%totalNodes
    totalnodes2 = gp%totalNodes2
 
    allocate(copy_core_count(totalparts))
    allocate(empty_parts(totalparts))
    !call srand(seed)
    write(*,*) "SA called..."
    if (totalnodes .lt. totalparts) then
      write(*,*) "ERROR: Number of parts cannot be greater than number of nodes."
      stop
    end if
 
    call prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm)
 
    call prg_costindex(cost, sumcubes, maxch,smooth_maxch, obj_fun)
    prev_cost = cost
 
    do it=1, niter
      call prg_rand_node(gp, node, seed)
      do j = xadj(node), xadj(node+1)-1
        neighbor = adjncy(j)
        part_backup = partnumber(neighbor)
        call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, partnumber(node))
        call prg_costindex(cost, sumcubes, maxch, smooth_maxch, obj_fun)
        prg_delta = cost - prev_cost
        call prg_accept_prob(it, prg_delta, r)
        call random_number(u)
        if (u > r) then ! reject
          call update_prg_costpartition(gp, xadj, adjncy, partnumber, core_count, ch_count, halo_count, sumcubes, maxch, smooth_maxch, pnorm, neighbor, part_backup)
        else
          prev_cost = cost
          !write(*,*) maxCH
        end if
      end do
 
    end do
 
 
 
    write(*,*) "Cost of meTIS+SA", sumcubes, maxch, smooth_maxch
 
 
 
    !prg_initialize and fill up subgraph structure
    !! Assign node ids (mapped to orbitals as rows) to each node in each
    call prg_destroygraphpartitioning(gp)
    write(pname, '("SAParts")')
    call prg_initgraphpartitioning(gp, pname, totalparts, totalnodes, totalnodes2)
    do i = 1, totalparts
      gp%nnodesInPartAll(i) = core_count(i)
      copy_core_count(i) = core_count(i)
      call prg_initsubgraph(gp%sgraph(i), i, gp%totalNodes2)
      allocate(gp%sgraph(i)%nodeInPart(core_count(i)))
      gp%nnodesInPart(i) = core_count(i)
    enddo
 
    !Assign node ids to sgraph
    do i=1, gp%totalNodes
      copy_core_count(partnumber(i)) =copy_core_count(partnumber(i)) - 1
      ! core_count((part(i))) - copy_core_count(part(i)) = node postion in array
      gp%sgraph(partnumber(i))%nodeInPart(core_count((partnumber(i))) - copy_core_count(partnumber(i)) ) = i -1 !NOTE: nodes in gp%sgraph()%nodeInPart() are currently 0 based!
    end do
 
 
    do i = 1, totalparts
      do j = 1, core_count(i)
        if( partnumber( gp%sgraph(i)%nodeInPart(j)+1 ) /= i) then
          write(*,*) "ERROR: subgraph struc incorrect!!", "node=",gp%sgraph(i)%nodeInPart(j)+1 , "part=",i, "actual_part=", partnumber(gp%sgraph(i)%nodeInPart(j)+1 )
          stop
        end if
      end do
    end do
    do i=1, totalparts
      write(*,*) "part=",i, "C=", core_count(i), "CH=", ch_count(i)
      if (ch_count(i) .eq. 0) then
        write(*,*) "ERROR: SA produced an empty part"
 
      end if
    end do
  end subroutine prg_simannealing_old
 
 
end module prg_partition_mod