hecmw_matrix_reorder.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 hecmw_matrix_reorder
  use hecmw_util
  use hecmw_array_util
  implicit none
 
  private
  public :: hecmw_matrix_reorder_profile
  public :: hecmw_matrix_reorder_values
  public :: hecmw_matrix_reorder_vector
  public :: hecmw_matrix_reorder_back_vector
  public :: hecmw_matrix_reorder_renum_item
 
contains
 
  subroutine hecmw_matrix_reorder_profile(N, perm, iperm, &
      indexL, indexU, itemL, itemU, indexLp, indexUp, itemLp, itemUp)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: perm(:), iperm(:)
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    integer(kind=kint), intent(out) :: indexlp(0:), indexup(0:)
    integer(kind=kint), intent(out) :: itemlp(:), itemup(:)
    integer(kind=kint) :: cntl, cntu, inew, iold, j, jold, jnew
    cntl = 0
    cntu = 0
    indexlp(0) = 0
    indexup(0) = 0
    do inew=1,n
      iold = perm(inew)
      ! original L
      do j = indexl(iold-1)+1, indexl(iold)
        jold = iteml(j)
        jnew = iperm(jold)
        if (jnew < inew) then
          cntl = cntl + 1
          itemlp(cntl) = jnew
        else
          cntu = cntu + 1
          itemup(cntu) = jnew
        end if
      end do
      ! original U
      do j = indexu(iold-1)+1, indexu(iold)
        jold = itemu(j)
        if (jold > n) cycle
        jnew = iperm(jold)
        if (jnew < inew) then
          cntl = cntl + 1
          itemlp(cntl) = jnew
        else
          cntu = cntu + 1
          itemup(cntu) = jnew
        end if
      end do
      indexlp(inew) = cntl
      indexup(inew) = cntu
      call hecmw_qsort_int_array(itemlp, indexlp(inew-1)+1, indexlp(inew))
      call hecmw_qsort_int_array(itemup, indexup(inew-1)+1, indexup(inew))
    end do
  end subroutine hecmw_matrix_reorder_profile
 
  subroutine hecmw_matrix_reorder_values(N, NDOF, perm, iperm, &
      indexL, indexU, itemL, itemU, AL, AU, D, &
      indexLp, indexUp, itemLp, itemUp, ALp, AUp, Dp)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: perm(:), iperm(:)
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    real(kind=kreal), intent(in) :: al(:), au(:), d(:)
    integer(kind=kint), intent(in) :: indexlp(0:), indexup(0:)
    integer(kind=kint), intent(in) :: itemlp(:), itemup(:)
    real(kind=kreal), intent(out) :: alp(:), aup(:), dp(:)
    dp = 0.d0
    alp = 0.d0
    aup = 0.d0
    ! call reorder_diag(N, NDOF, perm, D, Dp)
    ! call reorder_off_diag(N, NDOF, perm, &
      !    indexL, indexU, itemL, itemU, AL, AU, &
      !    indexLp, itemLp, ALp)
    ! call reorder_off_diag(N, NDOF, perm, &
      !    indexL, indexU, itemL, itemU, AL, AU, &
      !    indexUp, itemUp, AUp)
    call reorder_diag2(n, ndof, iperm, d, dp)
    call reorder_off_diag2(n, ndof, iperm, &
      indexl, iteml, al, &
      indexlp, indexup, itemlp, itemup, alp, aup)
    call reorder_off_diag2(n, ndof, iperm, &
      indexu, itemu, au, &
      indexlp, indexup, itemlp, itemup, alp, aup)
  end subroutine hecmw_matrix_reorder_values
 
  subroutine hecmw_matrix_reorder_vector(N, NDOF, perm, X, Xp)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: perm(:)
    real(kind=kreal), intent(in) :: x(:)
    real(kind=kreal), intent(out) :: xp(:)
    integer(kind=kint) :: inew, iold, j0new, j0old, j
    !$omp parallel default(none),private(inew,iold,j0new,j0old,j), &
      !$omp   shared(N,perm,NDOF,Xp,X)
    !$omp do
    do inew=1,n
      iold = perm(inew)
      j0new = (inew-1)*ndof
      j0old = (iold-1)*ndof
      do j=1,ndof
        xp(j0new + j) = x(j0old + j)
      end do
    end do
    !$omp end do
    !$omp end parallel
  end subroutine hecmw_matrix_reorder_vector
 
  subroutine hecmw_matrix_reorder_back_vector(N, NDOF, perm, Xp, X)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: perm(:)
    real(kind=kreal), intent(in) :: xp(:)
    real(kind=kreal), intent(out) :: x(:)
    integer(kind=kint) :: inew, iold, j0new, j0old, j
    !$omp parallel default(none),private(inew,iold,j0new,j0old,j), &
      !$omp&  shared(N,perm,NDOF,X,Xp)
    !$omp do
    do inew=1,n
      iold = perm(inew)
      j0new = (inew-1)*ndof
      j0old = (iold-1)*ndof
      do j=1,ndof
        x(j0old + j) = xp(j0new + j)
      end do
    end do
    !$omp end do
    !$omp end parallel
  end subroutine hecmw_matrix_reorder_back_vector
 
  subroutine hecmw_matrix_reorder_renum_item(N, perm, indexXp, itemXp)
    implicit none
    integer(kind=kint), intent(in) :: n
    integer(kind=kint), intent(in) :: perm(:)
    integer(kind=kint), intent(in) :: indexxp(0:)
    integer(kind=kint), intent(inout) :: itemxp(:)
    integer(kind=kint) :: npx, i
    npx = indexxp(n)
    !$omp parallel default(none),private(i),shared(NPX,itemXp,perm)
    !$omp do
    do i=1,npx
      itemxp(i) = perm( itemxp(i) )
    end do
    !$omp end do
    !$omp end parallel
  end subroutine hecmw_matrix_reorder_renum_item
 
  subroutine reorder_diag(N, NDOF, perm, D, Dp)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: perm(:)
    real(kind=kreal), intent(in) :: d(:)
    real(kind=kreal), intent(out) :: dp(:)
    integer(kind=kint) :: ndof2, inew, iold, j0new, j0old, j
    ndof2 = ndof*ndof
    ! diagonal
    do inew=1,n
      iold = perm(inew)
      j0new = (inew-1)*ndof2
      j0old = (iold-1)*ndof2
      do j=1,ndof2
        dp(j0new + j) = d(j0old + j)
      end do
    end do
  end subroutine reorder_diag
 
  subroutine reorder_diag2(N, NDOF, iperm, D, Dp)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: iperm(:)
    real(kind=kreal), intent(in) :: d(:)
    real(kind=kreal), intent(out) :: dp(:)
    integer(kind=kint) :: ndof2, inew, iold, j0new, j0old, j
    ndof2 = ndof*ndof
    ! diagonal
    !$omp parallel default(none),private(iold,inew,j0old,j0new,j), &
      !$omp&         shared(N,iperm,NDOF2,Dp,D)
    !$omp do
    do iold=1,n
      inew = iperm(iold)
      j0old = (iold-1)*ndof2
      j0new = (inew-1)*ndof2
      do j=1,ndof2
        dp(j0new + j) = d(j0old + j)
      end do
    end do
    !$omp end do
    !$omp end parallel
  end subroutine reorder_diag2
 
  subroutine reorder_off_diag(N, NDOF, perm, &
      indexL, indexU, itemL, itemU, AL, AU, &
      indexXp, itemXp, AXp)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: perm(:)
    integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
    integer(kind=kint), intent(in) :: iteml(:), itemu(:)
    real(kind=kreal), intent(in) :: al(:), au(:)
    integer(kind=kint), intent(in) :: indexxp(0:)
    integer(kind=kint), intent(in) :: itemxp(:)
    real(kind=kreal), intent(out) :: axp(:)
    integer(kind=kint) :: ndof2, inew, iold
    integer(kind=kint) :: jsoldl, jeoldl, jsoldu, jeoldu
    integer(kind=kint) :: jnew, knew, kold, jold, l0new, l0old, l
    ndof2 = ndof*ndof
    ! new L
    do inew=1,n
      iold = perm(inew)
      jsoldl = indexl(iold-1)+1
      jeoldl = indexl(iold)
      jsoldu = indexu(iold-1)+1
      jeoldu = indexu(iold)
      do jnew = indexxp(inew-1)+1, indexxp(inew)
        knew = itemxp(jnew)
        kold = perm(knew)
        if (kold < iold) then
          call hecmw_bsearch_int_array(iteml, jsoldl, jeoldl, kold, jold)
          if (jold < 0) then
            write(0,*) 'DEBUG:: jold < 0 in reorder_off_diag'
            cycle
          end if
          l0new = (jnew-1)*ndof2
          l0old = (jold-1)*ndof2
          do l=1,ndof2
            axp(l0new + l) = al(l0old + l)
          end do
        else
          call hecmw_bsearch_int_array(itemu, jsoldu, jeoldu, kold, jold)
          if (jold < 0) then
            write(0,*) 'DEBUG:: jold < 0 in reorder_off_diag'
            cycle
          end if
          l0new = (jnew-1)*ndof2
          l0old = (jold-1)*ndof2
          do l=1,ndof2
            axp(l0new + l) = au(l0old + l)
          end do
        end if
      end do
    end do
  end subroutine reorder_off_diag
 
  subroutine reorder_off_diag2(N, NDOF, iperm, &
      indexX, itemX, AX, &
      indexLp, indexUp, itemLp, itemUp, ALp, AUp)
    implicit none
    integer(kind=kint), intent(in) :: n, ndof
    integer(kind=kint), intent(in) :: iperm(:)
    integer(kind=kint), intent(in) :: indexx(0:)
    integer(kind=kint), intent(in) :: itemx(:)
    real(kind=kreal), intent(in) :: ax(:)
    integer(kind=kint), intent(in) :: indexlp(0:), indexup(0:)
    integer(kind=kint), intent(in) :: itemlp(:), itemup(:)
    real(kind=kreal), intent(out) :: alp(:), aup(:)
    integer(kind=kint) :: ndof2, iold, inew
    integer(kind=kint) :: jsnewl, jenewl, jsnewu, jenewu
    integer(kind=kint) :: jold, kold, knew, jnew, l0old, l0new, l
    ndof2 = ndof*ndof
    ! new L
    !$omp parallel default(none), &
      !$omp   private(iold,inew,jsnewL,jenewL,jsnewU,jenewU,jold,kold,knew,jnew,l0old,l0new,l), &
      !$omp   shared(N,iperm,indexLp,indexUp,indexX,itemX,itemLp,NDOF2,ALp,AX,itemUp,AUp)
    !$omp do
    do iold=1,n
      inew = iperm(iold)
      jsnewl = indexlp(inew-1)+1
      jenewl = indexlp(inew)
      jsnewu = indexup(inew-1)+1
      jenewu = indexup(inew)
      do jold = indexx(iold-1)+1, indexx(iold)
        kold = itemx(jold)
        if (kold > n) cycle
        knew = iperm(kold)
        if (knew < inew) then
          call hecmw_bsearch_int_array(itemlp, jsnewl, jenewl, knew, jnew)
          if (jnew < 0) then
            write(0,*) 'ERROR:: jnew < 0 in reorder_off_diag2'
            call hecmw_abort( hecmw_comm_get_comm() )
          end if
          l0old = (jold-1)*ndof2
          l0new = (jnew-1)*ndof2
          do l=1,ndof2
            alp(l0new + l) = ax(l0old + l)
          end do
        else
          call hecmw_bsearch_int_array(itemup, jsnewu, jenewu, knew, jnew)
          if (jnew < 0) then
            write(0,*) 'ERROR:: jnew < 0 in reorder_off_diag2'
            call hecmw_abort( hecmw_comm_get_comm() )
          end if
          l0old = (jold-1)*ndof2
          l0new = (jnew-1)*ndof2
          do l=1,ndof2
            aup(l0new + l) = ax(l0old + l)
          end do
        end if
      end do
    end do
    !$omp end do
    !$omp end parallel
  end subroutine reorder_off_diag2
 
end module hecmw_matrix_reorder