hecmw_precond_SSOR_nn.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
!C
!C***
!C*** module hecmw_precond_SSOR_nn
!C***
!C
module hecmw_precond_ssor_nn
  use hecmw_util
  use hecmw_matrix_misc
  use m_hecmw_matrix_ordering_cm
  use m_hecmw_matrix_ordering_mc
  use hecmw_matrix_reorder
  !$ use omp_lib
 
  private
 
  public:: hecmw_precond_ssor_nn_setup
  public:: hecmw_precond_ssor_nn_apply
  public:: hecmw_precond_ssor_nn_clear
 
  integer(kind=kint) :: N
  real(kind=kreal), pointer :: d(:) => null()
  real(kind=kreal), pointer :: al(:) => null()
  real(kind=kreal), pointer :: au(:) => null()
  integer(kind=kint), pointer :: indexL(:) => null()
  integer(kind=kint), pointer :: indexU(:) => null()
  integer(kind=kint), pointer :: itemL(:) => null()
  integer(kind=kint), pointer :: itemU(:) => null()
  real(kind=kreal), pointer :: alu(:) => null()
 
  integer(kind=kint) :: NColor
  integer(kind=kint), pointer :: COLORindex(:) => null()
  integer(kind=kint), pointer :: perm(:) => null()
  integer(kind=kint), pointer :: iperm(:) => null()
 
  logical, save :: isFirst = .true.
 
  logical, save :: INITIALIZED = .false.
 
contains
 
  subroutine hecmw_precond_ssor_nn_setup(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: npl, npu, npcl, npcu
    real   (kind=kreal), allocatable :: cd(:)
    integer(kind=kint ) :: ncolor_in
    real   (kind=kreal) :: sigma_diag
    real   (kind=kreal) :: alutmp(hecmat%NDOF,hecmat%NDOF), pw(hecmat%NDOF)
    integer(kind=kint ) :: ii, i, j, k, ndof, ndof2
    integer(kind=kint ) :: nthreads = 1
    integer(kind=kint ), allocatable :: perm_tmp(:)
    !real   (kind=kreal) :: t0
 
    !t0 = hecmw_Wtime()
    !write(*,*) 'DEBUG: SSOR setup start', hecmw_Wtime()-t0
 
    if (initialized) then
      if (hecmat%Iarray(98) == 1) then ! need symbolic and numerical setup
        call hecmw_precond_ssor_nn_clear(hecmat)
      else if (hecmat%Iarray(97) == 1) then ! need numerical setup only
        call hecmw_precond_ssor_nn_clear(hecmat) ! TEMPORARY
      else
        return
      endif
    endif
 
    !$ nthreads = omp_get_max_threads()
 
    n = hecmat%N
    ndof=hecmat%NDOF
    ndof2=ndof*ndof
 
    ncolor_in = hecmw_mat_get_ncolor_in(hecmat)
    sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
 
    if (nthreads == 1) then
      ncolor = 1
      allocate(colorindex(0:1), perm(n), iperm(n))
      colorindex(0) = 0
      colorindex(1) = n
      do i=1,n
        perm(i) = i
        iperm(i) = i
      end do
    else
      allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
      call hecmw_matrix_ordering_rcm(n, hecmat%indexL, hecmat%itemL, &
        hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
      !write(*,*) 'DEBUG: RCM ordering done', hecmw_Wtime()-t0
      call hecmw_matrix_ordering_mc(n, hecmat%indexL, hecmat%itemL, &
        hecmat%indexU, hecmat%itemU, perm_tmp, &
        ncolor_in, ncolor, colorindex, perm, iperm)
      !write(*,*) 'DEBUG: MC ordering done', hecmw_Wtime()-t0
      deallocate(perm_tmp)
 
      !call write_debug_info
    endif
 
    npl = hecmat%indexL(n)
    npu = hecmat%indexU(n)
    allocate(indexl(0:n), indexu(0:n), iteml(npl), itemu(npu))
    call hecmw_matrix_reorder_profile(n, perm, iperm, &
      hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
      indexl, indexu, iteml, itemu)
    !write(*,*) 'DEBUG: reordering profile done', hecmw_Wtime()-t0
 
    !call check_ordering
 
    allocate(d(ndof2*n), al(ndof2*npl), au(ndof2*npu))
    call hecmw_matrix_reorder_values(n, ndof, perm, iperm, &
      hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
      hecmat%AL, hecmat%AU, hecmat%D, &
      indexl, indexu, iteml, itemu, al, au, d)
    !write(*,*) 'DEBUG: reordering values done', hecmw_Wtime()-t0
 
    call hecmw_matrix_reorder_renum_item(n, perm, indexl, iteml)
    call hecmw_matrix_reorder_renum_item(n, perm, indexu, itemu)
 
    allocate(alu(ndof2*n))
    alu  = 0.d0
 
    do ii= 1, ndof2*n
      alu(ii) = d(ii)
    enddo
 
    !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,NDOF,NDOF2,ALU,SIGMA_DIAG)
    !$omp do
    do ii= 1, n
      do i = 1, ndof
        do j =  1, ndof
          alutmp(i,j) = alu(ndof2*(ii-1)+(i-1)*ndof+j)
          if (i==j) alutmp(i,j)=alutmp(i,j)*sigma_diag
        end do
      end do
      do k= 1, ndof
        alutmp(k,k)= 1.d0/alutmp(k,k)
        do i= k+1, ndof
          alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
          do j= k+1, ndof
            pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
          enddo
          do j= k+1, ndof
            alutmp(i,j)= pw(j)
          enddo
        enddo
      enddo
      do i = 1, ndof
        do j =  1, ndof
          alu(ndof2*(ii-1)+(i-1)*ndof+j)= alutmp(i,j)
        end do
      end do
    enddo
    !$omp end do
    !$omp end parallel
 
    isfirst = .true.
 
    initialized = .true.
    hecmat%Iarray(98) = 0 ! symbolic setup done
    hecmat%Iarray(97) = 0 ! numerical setup done
 
    !write(*,*) 'DEBUG: SSOR setup done', hecmw_Wtime()-t0
 
  end subroutine hecmw_precond_ssor_nn_setup
 
  subroutine hecmw_precond_ssor_nn_apply(ZP, NDOF)
    use hecmw_tuning_fx
    implicit none
    real(kind=kreal), intent(inout) :: zp(:)
    integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k, ndof, ndof2, idof,jdof
    real(kind=kreal) :: sw(ndof), x(ndof)
 
    ! added for turning >>>
    integer(kind=kint), parameter :: numofblockperthread = 100
    integer(kind=kint), save :: numofthread = 1, numofblock
    integer(kind=kint), save, allocatable :: ictoblockindex(:)
    integer(kind=kint), save, allocatable :: blockindextocolorindex(:)
    integer(kind=kint), save :: sectorcachesize0, sectorcachesize1
    integer(kind=kint) :: blockindex, elementcount, numofelement, ii
    real(kind=kreal) :: numofelementperblock
    integer(kind=kint) :: my_rank
 
    ndof2=ndof*ndof
    if (isfirst) then
      !$ numOfThread = omp_get_max_threads()
      numofblock = numofthread * numofblockperthread
      if (allocated(ictoblockindex)) deallocate(ictoblockindex)
      if (allocated(blockindextocolorindex)) deallocate(blockindextocolorindex)
      allocate (ictoblockindex(0:ncolor), &
        blockindextocolorindex(0:numofblock + ncolor))
      numofelement = n + indexl(n) + indexu(n)
      numofelementperblock = dble(numofelement) / numofblock
      blockindex = 0
      ictoblockindex = -1
      ictoblockindex(0) = 0
      blockindextocolorindex = -1
      blockindextocolorindex(0) = 0
      my_rank = hecmw_comm_get_rank()
      ! write(9000+my_rank,*) &
        !      '# numOfElementPerBlock =', numOfElementPerBlock
      ! write(9000+my_rank,*) &
        !      '# ic, blockIndex, colorIndex, elementCount'
      do ic = 1, ncolor
        elementcount = 0
        ii = 1
        do i = colorindex(ic-1)+1, colorindex(ic)
          elementcount = elementcount + 1
          elementcount = elementcount + (indexl(i) - indexl(i-1))
          elementcount = elementcount + (indexu(i) - indexu(i-1))
          if (elementcount > ii * numofelementperblock &
              .or. i == colorindex(ic)) then
            ii = ii + 1
            blockindex = blockindex + 1
            blockindextocolorindex(blockindex) = i
            ! write(9000+my_rank,*) ic, blockIndex, &
              !      blockIndexToColorIndex(blockIndex), elementCount
          endif
        enddo
        ictoblockindex(ic) = blockindex
      enddo
      numofblock = blockindex
 
      call hecmw_tuning_fx_calc_sector_cache( n, ndof, &
        sectorcachesize0, sectorcachesize1 )
 
      isfirst = .false.
    endif
    ! <<< added for turning
 
    !call start_collection("loopInPrecond33")
 
    !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
    !OCL CACHE_SUBSECTOR_ASSIGN(ZP)
 
    !$omp parallel default(none) &
      !$omp&shared(NColor,indexL,itemL,indexU,itemU,AL,AU,D,ALU,perm,&
      !$omp&       ZP,icToBlockIndex,blockIndexToColorIndex,NDOF,NDOF2) &
      !$omp&private(SW,X,ic,i,iold,isL,ieL,isU,ieU,j,k,blockIndex,idof,jdof)
 
    !C-- FORWARD
    do ic=1,ncolor
      !$omp do schedule (static, 1)
      do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
        do i = blockindextocolorindex(blockindex-1)+1, &
            blockindextocolorindex(blockindex)
          iold = perm(i)
          do idof = 1, ndof
            sw(idof) = zp(ndof*(iold-1)+idof)
          end do
          isl= indexl(i-1)+1
          iel= indexl(i)
          do j= isl, iel
            k= iteml(j)
            do idof = 1, ndof
              x(idof) = zp(ndof*(k-1)+idof)
            end do
            do idof = 1, ndof
              do jdof = 1, ndof
                sw(idof) = sw(idof) - al(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
              end do
            end do
          enddo ! j
 
          x = sw
          do idof = 2,ndof
            do jdof = 1, idof-1
              x(idof) = x(idof) - alu(ndof2*(i-1)+ndof*(idof-1)+jdof )*x(jdof)
            end do
          end do
          do idof = ndof, 1, -1
            do jdof = ndof, idof+1, -1
              x(idof) = x(idof) - alu(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
            end do
            x(idof) = alu(ndof2*(i-1)+(ndof+1)*(idof-1)+1)*x(idof)
          end do
          zp(ndof*(iold-1)+1:ndof*(iold-1)+ndof) = x(1:ndof)
 
        enddo ! i
      enddo ! blockIndex
      !$omp end do
    enddo ! ic
 
    !C-- BACKWARD
    do ic=ncolor, 1, -1
      !$omp do schedule (static, 1)
      do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
        do i = blockindextocolorindex(blockindex), &
            blockindextocolorindex(blockindex-1)+1, -1
          ! do blockIndex = icToBlockIndex(ic-1)+1, icToBlockIndex(ic)
          !   do i = blockIndexToColorIndex(blockIndex-1)+1, &
            !        blockIndexToColorIndex(blockIndex)
          !   do i = endPos(threadNum, ic), startPos(threadNum, ic), -1
          sw= 0.d0
          isu= indexu(i-1) + 1
          ieu= indexu(i)
          do j= ieu, isu, -1
            k= itemu(j)
            do idof = 1, ndof
              x(idof) = zp(ndof*(k-1)+idof)
            end do
            do idof = 1, ndof
              do jdof = 1, ndof
                sw(idof) = sw(idof) + au(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
              end do
            end do
          enddo ! j
 
          x = sw
          do idof = 2, ndof
            do k = 1,idof-1
              x(idof) = x(idof) - alu(ndof2*(i-1)+ndof*(idof-1)+k)*x(k)
            end do
          end do
          do idof = ndof, 1, -1
            do k = ndof, idof+1, -1
              x(idof) = x(idof) - alu(ndof2*(i-1)+ndof*(idof-1)+k)*x(k)
            end do
            x(idof) = alu(ndof2*(i-1)+(ndof+1)*(idof-1)+1)*x(idof)
          end do
          iold = perm(i)
          do idof = 1, ndof
            zp(ndof*(iold-1)+idof) = zp(ndof*(iold-1)+idof) - x(idof)
          end do
        enddo ! i
      enddo ! blockIndex
      !$omp end do
    enddo ! ic
    !$omp end parallel
 
    !OCL END_CACHE_SUBSECTOR
    !OCL END_CACHE_SECTOR_SIZE
 
    !call stop_collection("loopInPrecond33")
 
  end subroutine hecmw_precond_ssor_nn_apply
 
  subroutine hecmw_precond_ssor_nn_clear(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: nthreads = 1
    !$ nthreads = omp_get_max_threads()
    if (associated(colorindex)) deallocate(colorindex)
    if (associated(perm)) deallocate(perm)
    if (associated(iperm)) deallocate(iperm)
    if (associated(alu)) deallocate(alu)
    if (nthreads >= 1) then
      if (associated(d)) deallocate(d)
      if (associated(al)) deallocate(al)
      if (associated(au)) deallocate(au)
      if (associated(indexl)) deallocate(indexl)
      if (associated(indexu)) deallocate(indexu)
      if (associated(iteml)) deallocate(iteml)
      if (associated(itemu)) deallocate(itemu)
    end if
    nullify(colorindex)
    nullify(perm)
    nullify(iperm)
    nullify(alu)
    nullify(d)
    nullify(al)
    nullify(au)
    nullify(indexl)
    nullify(indexu)
    nullify(iteml)
    nullify(itemu)
    initialized = .false.
  end subroutine hecmw_precond_ssor_nn_clear
 
  subroutine write_debug_info
    implicit none
    integer(kind=kint) :: my_rank, ic, in
    my_rank = hecmw_comm_get_rank()
    !--------------------> debug: shizawa
    if (my_rank.eq.0) then
      write(*,*) 'DEBUG: Output fort.19000+myrank and fort.29000+myrank for coloring information'
    endif
    write(19000+my_rank,'(a)') '#NCOLORTot'
    write(19000+my_rank,*) ncolor
    write(19000+my_rank,'(a)') '#ic  COLORindex(ic-1)+1  COLORindex(ic)'
    do ic=1,ncolor
      write(19000+my_rank,*) ic, colorindex(ic-1)+1,colorindex(ic)
    enddo ! ic
    write(29000+my_rank,'(a)') '#n_node'
    write(29000+my_rank,*) n
    write(29000+my_rank,'(a)') '#in  OLDtoNEW(in)  NEWtoOLD(in)'
    do in=1,n
      write(29000+my_rank,*) in, iperm(in), perm(in)
      if (perm(iperm(in)) .ne. in) then
        write(29000+my_rank,*) '** WARNING **: NEWtoOLD and OLDtoNEW: ',in
      endif
    enddo
  end subroutine write_debug_info
 
  subroutine check_ordering
    implicit none
    integer(kind=kint) :: ic, i, j, k
    integer(kind=kint), allocatable :: iicolor(:)
    ! check color dependence of neighbouring nodes
    if (ncolor.gt.1) then
      allocate(iicolor(n))
      do ic=1,ncolor
        do i= colorindex(ic-1)+1, colorindex(ic)
          iicolor(i) = ic
        enddo ! i
      enddo ! ic
      ! FORWARD: L-part
      do ic=1,ncolor
        do i= colorindex(ic-1)+1, colorindex(ic)
          do j= indexl(i-1)+1, indexl(i)
            k= iteml(j)
            if (iicolor(i).eq.iicolor(k)) then
              write(*,*) .eq.'** ERROR **: L-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
            endif
          enddo ! j
        enddo ! i
      enddo ! ic
      ! BACKWARD: U-part
      do ic=ncolor, 1, -1
        do i= colorindex(ic), colorindex(ic-1)+1, -1
          do j= indexu(i-1)+1, indexu(i)
            k= itemu(j)
            if (iicolor(i).eq.iicolor(k)) then
              write(*,*) .eq.'** ERROR **: U-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
            endif
          enddo ! j
        enddo ! i
      enddo ! ic
      deallocate(iicolor)
    endif ! if (NColor.gt.1)
    !--------------------< debug: shizawa
  end subroutine check_ordering
 
end module     hecmw_precond_ssor_nn