hecmw_matrix_ordering_CM.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module m_hecmw_matrix_ordering_cm
  use hecmw_util
  implicit none
 
  private
  public :: hecmw_matrix_ordering_cm
  public :: hecmw_matrix_ordering_rcm
 
  integer(kind=kint), parameter :: DEBUG = 0
 
contains
 
  subroutine hecmw_matrix_ordering_cm(N, indexL, itemL, indexU, itemU, &
      perm, iperm)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    integer(kind=kint), intent(out) :: perm(:), iperm(:)
    integer(kind=kint), parameter :: nminmax = 5
    integer(kind=kint) :: i, nmin
    integer(kind=kint) :: mins(nminmax)
    integer(kind=kint), allocatable :: degs(:)
    integer(kind=kint), allocatable :: nlevel(:)
    integer(kind=kint), allocatable :: lv_index(:,:), lv_item(:,:)
    integer(kind=kint) :: nlevel_max, max_id
    allocate(degs(n))
    call count_degrees(n, indexl, indexu, itemu, degs)
    call find_minimum_degrees(n, degs, nminmax, nmin, mins)
    allocate(nlevel(nmin), lv_index(0:n,nmin), lv_item(n,nmin))
    ! perform CM ordering starting from each minimum degree node
    !$omp parallel default(none),private(i), &
      !$omp&  shared(nmin,N,indexL,itemL,indexU,itemU,degs,mins,nlevel,lv_index,lv_item)
    !$omp do
    do i=1,nmin
      call ordering_cm_inner(n, indexl, iteml, indexu, itemu, degs, mins(i), &
        nlevel(i), lv_index(:,i), lv_item(:,i))
      if (debug > 0) write(*,*) 'DEBUG:: hecmw_matrix_ordering_CM: i, nstart, nlevel = ', i, mins(i), nlevel(i)
    end do
    !$omp end do
    !$omp end parallel
    deallocate(degs)
    ! choose CM ordering with maximum nlevel
    nlevel_max = nlevel(1)
    max_id = 1
    do i=2,nmin
      if (nlevel(i) > nlevel_max) then
        nlevel_max = nlevel(i)
        max_id = i
      end if
    end do
    if (debug > 0) write(*,*) 'DEBUG:: hecmw_matrix_ordering_CM: chose ordering',max_id
    do i=1,n
      perm(i) = lv_item(i,max_id)
      iperm(perm(i)) = i
    end do
    deallocate(nlevel, lv_index, lv_item)
  end subroutine hecmw_matrix_ordering_cm
 
  subroutine hecmw_matrix_ordering_rcm(N, indexL, itemL, indexU, itemU, &
      perm, iperm)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    integer(kind=kint), intent(out) :: perm(:), iperm(:)
    call hecmw_matrix_ordering_cm(n, indexl, iteml, indexu, itemu, perm, iperm)
    call reverse_ordering(n, perm, iperm)
    if (debug > 0) then
      call write_nonzero_profile(n, indexl, iteml, indexu, itemu, perm, iperm)
      call write_perm(n, perm, iperm)
    endif
  end subroutine hecmw_matrix_ordering_rcm
 
  subroutine ordering_cm_inner(N, indexL, itemL, indexU, itemU, degs, nstart, &
      nlevel, lv_index, lv_item)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    integer(kind=kint), intent(in) :: degs(:)
    integer(kind=kint), intent(in) :: nstart
    integer(kind=kint), intent(out) :: nlevel
    integer(kind=kint), intent(out) :: lv_index(0:)
    integer(kind=kint), intent(out) :: lv_item(:)
    integer(kind=kint), allocatable :: iwk(:)
    integer(kind=kint) :: level, cntall, cnt, j, jnode, k, knode
    allocate(iwk(n))
    iwk(:) = 0
    lv_index(0) = 0
    ! first level
    iwk(nstart) = 1
    cntall = 1
    lv_item(1) = nstart
    lv_index(1) = 1
    ! other levels
    do level=2,n
      cnt = 0
      ! all nodes in previous level
      prlv: do j = lv_index(level-2)+1, lv_index(level-1)
        jnode = lv_item(j)
        ! all connected nodes
        do k = indexl(jnode-1)+1, indexl(jnode)
          knode = iteml(k)
          if (iwk(knode) == 0) then
            iwk(knode) = level
            cnt = cnt + 1
            cntall = cntall + 1
            lv_item(cntall) = knode
            !if (cntall == N) exit PRLV
          end if
        end do
        do k = indexu(jnode-1)+1, indexu(jnode)
          knode = itemu(k)
          if (knode > n) cycle
          if (iwk(knode) == 0) then
            iwk(knode) = level
            cnt = cnt + 1
            cntall = cntall + 1
            lv_item(cntall) = knode
            !if (cntall == N) exit PRLV
          end if
        end do
      end do prlv
      if (cnt == 0) then
        if (debug > 0) write(*,*) 'DEBUG: choose any uncolored node..'
        do knode = 1, n
          if (iwk(knode) == 0) then
            iwk(knode) = level
            cnt = cnt + 1
            cntall = cntall + 1
            lv_item(cntall) = knode
            exit
          end if
        end do
      endif
      lv_index(level) = cntall
      call sort_nodes_by_degree(lv_item, lv_index(level-1)+1, lv_index(level), degs)
      if (cntall == n) then
        nlevel = level
        exit
      end if
    end do
    if (debug > 0) then
      level = 0
      do j = 1, n
        jnode = lv_item(j)
        if (jnode <= 0 .or. jnode > n) stop 'ERROR: ordering_CM_inner: out of range'
        if (iwk(jnode) < 0) stop 'ERROR: ordering_CM_inner: duplicated node found'
        if (iwk(jnode) < level) stop 'ERROR: ordering_CM_inner: not in level order'
        level = iwk(jnode)
        iwk(jnode) = -1
      enddo
      do j = 1, n
        if (iwk(j) /= -1) stop 'ERROR: ordering_CM_inner: non-numbered node found'
      enddo
    endif
    deallocate(iwk)
  end subroutine ordering_cm_inner
 
  subroutine count_degrees(N, indexL, indexU, itemU, degs)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: itemu(:)
    integer(kind=kint), intent(out) :: degs(:)
    integer(kind=kint) :: i, j
    do i=1,n
      degs(i) = indexl(i) - indexl(i-1)
      do j = indexu(i-1)+1, indexu(i)
        if (itemu(j) <= n) degs(i) = degs(i) + 1
      end do
    end do
  end subroutine count_degrees
 
  subroutine find_minimum_degrees(N, degs, nminmax, nmin, mins)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: degs(:)
    integer(kind=kint), intent(in) :: nminmax
    integer(kind=kint), intent(out) :: nmin
    integer(kind=kint), intent(out) :: mins(nminmax)
    integer(kind=kint) :: degmin, i, j
    degmin = n
    nmin = 0
    do i=1,n
      ! skip unconnected nodes
      if (degs(i) == 0) cycle
      if (degs(i) < degmin) then
        degmin = degs(i)
        nmin = 1
        mins(1) = i
      else if (degs(i) == degmin) then
        nmin = nmin + 1
        if (nmin <= nminmax) mins(nmin) = i
      end if
    end do
    if (debug > 0) write(*,*) 'DEBUG:: find_minimum_degrees: nmin, deg = ', nmin, degmin
    if (nmin > nminmax) nmin = nminmax
  end subroutine find_minimum_degrees
 
  recursive subroutine sort_nodes_by_degree(lv_item, istart, iend, degs)
    implicit none
    integer(kind=kint), intent(inout) :: lv_item(:)
    integer(kind=kint), intent(in) :: istart, iend
    integer(kind=kint), intent(in) :: degs(:)
    integer(kind=kint) :: pivot, left, right, itmp
    if (istart >= iend) return
    pivot = degs(lv_item((istart + iend) / 2))
    left = istart
    right = iend
    do
      do while (degs(lv_item(left)) < pivot)
        left = left + 1
      end do
      do while (pivot < degs(lv_item(right)))
        right = right - 1
      end do
      if (left >= right) exit
      itmp = lv_item(left)
      lv_item(left) = lv_item(right)
      lv_item(right) = itmp
      left = left + 1
      right = right - 1
    end do
    call sort_nodes_by_degree(lv_item, istart, left-1, degs)
    call sort_nodes_by_degree(lv_item, right+1, iend, degs)
  end subroutine sort_nodes_by_degree
 
  subroutine reverse_ordering(N, perm, iperm)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(inout) :: perm(:), iperm(:)
    integer(kind=kint) :: i, n1
    n1 = n + 1
    do i=1,n
      iperm(i) = n1 - iperm(i)
      perm(iperm(i)) = i
    end do
  end subroutine reverse_ordering
 
  subroutine write_nonzero_profile(N, indexL, itemL, indexU, itemU, perm, iperm)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    integer(kind=kint), intent(in) :: perm(:), iperm(:)
    integer(kind=kint), parameter :: f_org = 901
    integer(kind=kint), parameter :: f_new = 902
    integer(kind=kint) :: i, j, irow, jcol
    open(f_org, file='nzprof_org.txt', status='replace')
    do irow = 1, n
      i = irow
      do j = indexl(i-1)+1, indexl(i)
        jcol = iteml(j)
        write(f_org,*) irow, jcol
      end do
      do j = indexu(i-1)+1, indexu(i)
        jcol = itemu(j)
        if (jcol > n) cycle
        write(f_org,*) irow, jcol
      end do
    end do
    close(f_org)
    open(f_new, file='nzprof_new.txt', status='replace')
    do irow = 1, n
      i = perm(irow)
      do j = indexl(i-1)+1, indexl(i)
        jcol = iteml(j)
        write(f_new,*) irow, iperm(jcol)
      end do
      do j = indexu(i-1)+1, indexu(i)
        jcol = itemu(j)
        if (jcol > n) cycle
        write(f_new,*) irow, iperm(jcol)
      end do
    end do
    close(f_new)
  end subroutine write_nonzero_profile
 
  subroutine write_perm(N, perm, iperm)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: perm(:), iperm(:)
    integer(kind=kint), parameter :: f_perm = 903
    integer(kind=kint) :: i
    open(f_perm, file='perm_iperm.txt', status='replace')
    do i = 1, n
      write(f_perm,*) i, perm(i), iperm(i)
    end do
    close(f_perm)
  end subroutine write_perm
 
end module m_hecmw_matrix_ordering_cm