hecmw_precond_SSOR_33.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_33
!C***
!C
module hecmw_precond_ssor_33
  use hecmw_util
  use hecmw_matrix_misc
  use m_hecmw_matrix_ordering_cm
  use m_hecmw_matrix_ordering_mc
  use hecmw_matrix_reorder
#ifndef _OPENACC
  !$ use omp_lib
#endif
 
  private
 
  public:: hecmw_precond_ssor_33_setup
  public:: hecmw_precond_ssor_33_apply
  public:: hecmw_precond_ssor_33_clear
 
  integer(kind=kint) :: N
#ifdef HECMW_ARCH_FX64
  integer(kind=kint) :: NP
#endif
  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()
 
  ! taken from hecMAT at setup: _apply receives only ZP and cannot read it back
  integer(kind=kint) :: precond_impl = hecmw_precond_impl_bsr
 
  ! a format may be selected on a build that carries no SSOR for it; saying so once
  ! per run keeps the fallback from looking like the requested format took effect
  logical, save :: precond_impl_missing_reported = .false.
 
  logical, save :: isFirst = .true.
 
  logical, save :: INITIALIZED = .false.
 
  ! relaxation parameter of
  !   M = 1/(2-OMEGA) (D/OMEGA + L) (D/OMEGA)^-1 (D/OMEGA + U)
  ! OMEGA = 1 is the classic (D+L) D^-1 (D+U).  _apply is handed only ZP, so setup
  ! leaves the value here for it.
  real(kind=kreal), save :: omega = 1.d0
 
  ! for tuning
  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), parameter :: SSOR_DEBUG = 0
 
contains
 
  subroutine hecmw_precond_ssor_33_setup(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: npl, npu
    integer(kind=kint ) :: ncolor_in
    real   (kind=kreal) :: omega_inv
    real   (kind=kreal) :: alutmp(3,3), pw(3)
    integer(kind=kint ) :: ii, i, j, k
    integer(kind=kint ) :: nthreads = 1
    integer(kind=kint ), allocatable :: perm_tmp(:)
    real   (kind=kreal) :: t0
 
    if (ssor_debug >= 1) then
      t0 = hecmw_wtime()
      write(*,*) 'DEBUG: SSOR setup start', hecmw_wtime()-t0
    endif
 
    if (initialized) then
      if (hecmat%Iarray(98) == 1) then ! need symbolic and numerical setup
        call hecmw_precond_ssor_33_clear(hecmat)
      else if (hecmat%Iarray(97) == 1) then ! need numerical setup only
        call hecmw_precond_ssor_33_clear(hecmat) ! TEMPORARY
      else
        return
      endif
    endif
 
#ifndef _OPENACC
    !$ nthreads = omp_get_max_threads()
#endif
 
    n = hecmat%N
    ! N = hecMAT%NP
    ncolor_in = hecmw_mat_get_ncolor_in(hecmat)
    omega = hecmw_mat_get_omega(hecmat)
    omega_inv = 1.d0 / omega
    precond_impl = hecmw_mat_get_precond_impl(hecmat)
 
#ifdef _OPENACC
    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)
    if (ssor_debug >= 1) 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)
    if (ssor_debug >= 1) write(*,*) 'DEBUG: MC ordering done', hecmw_wtime()-t0
    deallocate(perm_tmp)
 
#else
    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)
      if (ssor_debug >= 1) 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)
      if (ssor_debug >= 1) write(*,*) 'DEBUG: MC ordering done', hecmw_wtime()-t0
      deallocate(perm_tmp)
 
    endif
#endif
 
    npl = 0
    do i=1,n
      do j=hecmat%indexU(i-1)+1,hecmat%indexU(i)
        if( hecmat%itemU(j) > n ) exit
        npl = npl + 1
      enddo
    enddo
    npl = max(hecmat%indexL(n),npl)
    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)
    if (ssor_debug >= 1) write(*,*) 'DEBUG: reordering profile done', hecmw_wtime()-t0
 
 
    allocate(d(9*n), al(9*npl), au(9*npu))
    call hecmw_matrix_reorder_values(n, 3, perm, iperm, &
      hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
      hecmat%AL, hecmat%AU, hecmat%D, &
      indexl, indexu, iteml, itemu, al, au, d)
    if (ssor_debug >= 1) 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(9*n))
    alu  = 0.d0
 
    do ii= 1, 9*n
      alu(ii) = d(ii)
    enddo
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent private(ALUtmp,PW)
#else
    !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,OMEGA_INV)
    !$omp do
#endif
    do ii= 1, n
      alutmp(1,1)= alu(9*ii-8) * omega_inv
      alutmp(1,2)= alu(9*ii-7) * omega_inv
      alutmp(1,3)= alu(9*ii-6) * omega_inv
      alutmp(2,1)= alu(9*ii-5) * omega_inv
      alutmp(2,2)= alu(9*ii-4) * omega_inv
      alutmp(2,3)= alu(9*ii-3) * omega_inv
      alutmp(3,1)= alu(9*ii-2) * omega_inv
      alutmp(3,2)= alu(9*ii-1) * omega_inv
      alutmp(3,3)= alu(9*ii  ) * omega_inv
      do k= 1, 3
        alutmp(k,k)= 1.d0/alutmp(k,k)
        do i= k+1, 3
          alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
          do j= k+1, 3
            pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
          enddo
          do j= k+1, 3
            alutmp(i,j)= pw(j)
          enddo
        enddo
      enddo
      alu(9*ii-8)= alutmp(1,1)
      alu(9*ii-7)= alutmp(1,2)
      alu(9*ii-6)= alutmp(1,3)
      alu(9*ii-5)= alutmp(2,1)
      alu(9*ii-4)= alutmp(2,2)
      alu(9*ii-3)= alutmp(2,3)
      alu(9*ii-2)= alutmp(3,1)
      alu(9*ii-1)= alutmp(3,2)
      alu(9*ii  )= alutmp(3,3)
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
 
    isfirst = .true.
 
    initialized = .true.
    hecmat%Iarray(98) = 0 ! symbolic setup done
    hecmat%Iarray(97) = 0 ! numerical setup done
 
    if (ssor_debug >= 1) write(*,*) 'DEBUG: SSOR setup done', hecmw_wtime()-t0
 
    if (precond_impl == hecmw_precond_impl_csr) call hecmw_precond_ssor_33_setup_csr(hecmat)
 
  end subroutine hecmw_precond_ssor_33_setup
 
  subroutine hecmw_precond_ssor_33_setup_csr(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(in) :: hecmat
#ifdef HECMW_ARCH_FX64
    integer(kind=kint ) :: npl, npu
    real   (kind=kreal) :: atmp(9)
    integer(kind=kint ) :: ii
    integer(kind=kint ), pointer :: item_tmp(:)
 
    np = hecmat%NP
 
    !$omp parallel default(none),private(ii,Atmp),shared(N,ALU)
    !$omp do
    do ii= 1, n
      atmp(1)= alu(9*ii-8)
      atmp(2)= alu(9*ii-7)
      atmp(3)= alu(9*ii-6)
      atmp(4)= alu(9*ii-5)
      atmp(5)= alu(9*ii-4)
      atmp(6)= alu(9*ii-3)
      atmp(7)= alu(9*ii-2)
      atmp(8)= alu(9*ii-1)
      atmp(9)= alu(9*ii  )
 
      alu(9*ii  )= atmp(9)
      alu(9*ii-1)= -atmp(8)*atmp(9)
      alu(9*ii-2)= (atmp(4)*atmp(8)-atmp(7))*atmp(9)
      alu(9*ii-3)= -atmp(5)*atmp(6)*alu(9*ii)
      alu(9*ii-4)= atmp(5)*(1.d0-atmp(6)*alu(9*ii-1))
      alu(9*ii-5)= -atmp(5)*(atmp(4)+atmp(6)*alu(9*ii-2))
      alu(9*ii-6)= -atmp(1)*(atmp(3)*alu(9*ii)+atmp(2)*alu(9*ii-3))
      alu(9*ii-7)= -atmp(1)*(atmp(3)*alu(9*ii-1)+atmp(2)*alu(9*ii-4))
      alu(9*ii-8)= atmp(1)*(1.d0-atmp(3)*alu(9*ii-2)-atmp(2)*alu(9*ii-5))
    enddo
    !$omp end do
    !$omp end parallel
 
    ! setup sized itemL and itemU for BSR; the expansion below takes three entries per block
    npl = indexl(n)
    npu = indexu(n)
    allocate(item_tmp(3*npl))
    item_tmp(1:npl) = iteml(1:npl)
    deallocate(iteml)
    iteml => item_tmp
    allocate(item_tmp(3*npu))
    item_tmp(1:npu) = itemu(1:npu)
    deallocate(itemu)
    itemu => item_tmp
 
    do ii=1,n
      indexl(ii) = 3*indexl(ii)
      indexu(ii) = 3*indexu(ii)
    end do
    do ii=npl,1,-1
      iteml(3*ii  ) = 3*iteml(ii)
      iteml(3*ii-1) = 3*iteml(ii)-1
      iteml(3*ii-2) = 3*iteml(ii)-2
    end do
    do ii=npu,1,-1
      itemu(3*ii  ) = 3*itemu(ii)
      itemu(3*ii-1) = 3*itemu(ii)-1
      itemu(3*ii-2) = 3*itemu(ii)-2
    end do
 
    !$omp parallel default(none),private(ii),shared(NPL,AL)
    !$omp do
    do ii=1,npl
      call hecmw_mat_block_transpose_33(al(9*ii-8:9*ii))
    end do
    !$omp end do
    !$omp end parallel
    !$omp parallel default(none),private(ii),shared(NPU,AU)
    !$omp do
    do ii=1,npu
      call hecmw_mat_block_transpose_33(au(9*ii-8:9*ii))
    end do
    !$omp end do
    !$omp end parallel
 
#endif
  end subroutine hecmw_precond_ssor_33_setup_csr
 
  subroutine setup_tuning_parameters
    use hecmw_tuning_fx
    implicit none
    integer(kind=kint) :: blockindex, elementcount, numofelement, ii
    real(kind=kreal) :: numofelementperblock
    integer(kind=kint) :: my_rank
    integer(kind=kint) :: ic, i
    if (ssor_debug >= 1) write(*,*) 'DEBUG: setting up tuning parameters for SSOR'
#ifndef _OPENACC
    !$ numOfThread = omp_get_max_threads()
#endif
 
    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, 3, &
         sectorcachesize0, sectorcachesize1 )
  end subroutine setup_tuning_parameters
 
  subroutine hecmw_precond_ssor_33_apply(ZP)
    use m_hecmw_comm_f
    implicit none
    real(kind=kreal), intent(inout) :: zp(:)
 
    select case (precond_impl)
      case (hecmw_precond_impl_bsr)
        call hecmw_precond_ssor_33_apply_generic(zp)
#ifdef HECMW_ARCH_FX64
      case (hecmw_precond_impl_csr)
        call hecmw_precond_ssor_33_apply_csr(zp)
#endif
      case default
        if (.not. precond_impl_missing_reported) then
          precond_impl_missing_reported = .true.
          if (hecmw_comm_get_rank() == 0) write(*,'(a)') &
            '#### MATRIXFORMAT: this build has no SSOR preconditioner for the selected '// &
            'format -- running the BSR implementation'
        endif
        call hecmw_precond_ssor_33_apply_generic(zp)
    end select
  end subroutine hecmw_precond_ssor_33_apply
 
  subroutine hecmw_precond_ssor_33_apply_generic(ZP)
    implicit none
    real(kind=kreal), intent(inout) :: zp(:)
    integer(kind=kint) :: ic, i, iold, j, isL, ieL, isU, ieU, k
    real(kind=kreal) :: sw1, sw2, sw3, x1, x2, x3
 
    real(kind=kreal) :: omega_fac
 
    ! added for tuning >>>
    integer(kind=kint) :: blockIndex
 
#ifndef _OPENACC
    if (isfirst) then
      call setup_tuning_parameters
      isfirst = .false.
    endif
#endif
    ! <<< added for tuning
 
    omega_fac = 2.d0 - omega
 
#ifndef _OPENACC
    !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,OMEGA_FAC) &
      !$omp&private(SW1,SW2,SW3,X1,X2,X3,ic,i,iold,isL,ieL,isU,ieU,j,k,blockIndex)
#endif
 
    !C-- FORWARD
    do ic=1,ncolor
#ifdef _OPENACC
      !$acc kernels
      !$acc loop independent
      do i = colorindex(ic-1)+1, colorindex(ic)
#else
      !$omp do schedule (static, 1)
      do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
        do i = blockindextocolorindex(blockindex-1)+1, &
            blockindextocolorindex(blockindex)
#endif
          ! do i = startPos(threadNum, ic), endPos(threadNum, ic)
          iold = perm(i)
          sw1= omega_fac * zp(3*iold-2)
          sw2= omega_fac * zp(3*iold-1)
          sw3= omega_fac * zp(3*iold  )
          isl= indexl(i-1)+1
          iel= indexl(i)
          do j= isl, iel
            !k= perm(itemL(j))
            k= iteml(j)
            x1= zp(3*k-2)
            x2= zp(3*k-1)
            x3= zp(3*k  )
            sw1= sw1 - al(9*j-8)*x1 - al(9*j-7)*x2 - al(9*j-6)*x3
            sw2= sw2 - al(9*j-5)*x1 - al(9*j-4)*x2 - al(9*j-3)*x3
            sw3= sw3 - al(9*j-2)*x1 - al(9*j-1)*x2 - al(9*j  )*x3
          enddo ! j
 
          x1= sw1
          x2= sw2
          x3= sw3
          x2= x2 - alu(9*i-5)*x1
          x3= x3 - alu(9*i-2)*x1 - alu(9*i-1)*x2
          x3= alu(9*i  )*  x3
          x2= alu(9*i-4)*( x2 - alu(9*i-3)*x3 )
          x1= alu(9*i-8)*( x1 - alu(9*i-6)*x3 - alu(9*i-7)*x2)
          zp(3*iold-2)= x1
          zp(3*iold-1)= x2
          zp(3*iold  )= x3
#ifdef _OPENACC
      enddo
      !$acc end kernels
#else
        enddo ! i
      enddo ! blockIndex
      !$omp end do
#endif
    enddo ! ic
 
    !C-- BACKWARD
    do ic=ncolor, 1, -1
#ifdef _OPENACC
      !$acc kernels
      !$acc loop independent
      do i = colorindex(ic-1)+1, colorindex(ic)
#else
      !$omp do schedule (static, 1)
      do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
        do i = blockindextocolorindex(blockindex), &
            blockindextocolorindex(blockindex-1)+1, -1
#endif
          ! 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
          sw1= 0.d0
          sw2= 0.d0
          sw3= 0.d0
          isu= indexu(i-1) + 1
          ieu= indexu(i)
          do j= ieu, isu, -1
            !k= perm(itemU(j))
            k= itemu(j)
            x1= zp(3*k-2)
            x2= zp(3*k-1)
            x3= zp(3*k  )
            sw1= sw1 + au(9*j-8)*x1 + au(9*j-7)*x2 + au(9*j-6)*x3
            sw2= sw2 + au(9*j-5)*x1 + au(9*j-4)*x2 + au(9*j-3)*x3
            sw3= sw3 + au(9*j-2)*x1 + au(9*j-1)*x2 + au(9*j  )*x3
          enddo ! j
 
          x1= sw1
          x2= sw2
          x3= sw3
          x2= x2 - alu(9*i-5)*x1
          x3= x3 - alu(9*i-2)*x1 - alu(9*i-1)*x2
          x3= alu(9*i  )*  x3
          x2= alu(9*i-4)*( x2 - alu(9*i-3)*x3 )
          x1= alu(9*i-8)*( x1 - alu(9*i-6)*x3 - alu(9*i-7)*x2)
          iold = perm(i)
          zp(3*iold-2)=  zp(3*iold-2) - x1
          zp(3*iold-1)=  zp(3*iold-1) - x2
          zp(3*iold  )=  zp(3*iold  ) - x3
#ifdef _OPENACC
      enddo
      !$acc end kernels
#else
        enddo ! i
      enddo ! blockIndex
      !$omp end do
#endif
    enddo ! ic
#ifndef _OPENACC
    !$omp end parallel
 
    !OCL END_CACHE_SUBSECTOR
    !OCL END_CACHE_SECTOR_SIZE
 
    !call stop_collection("loopInPrecond33")
#endif
 
  end subroutine hecmw_precond_ssor_33_apply_generic
 
  subroutine hecmw_precond_ssor_33_apply_csr(ZP)
    use hecmw_tuning_fx
    implicit none
    real(kind=kreal), intent(inout) :: zp(:)
#ifdef HECMW_ARCH_FX64
    integer(kind=kint) :: ic, i
 
    real(kind=kreal) :: omega_fac
 
    ! added for turning >>>
    integer(kind=kint) :: blockIndex
 
    if (isfirst) then
      call setup_tuning_parameters
      isfirst = .false.
    endif
    ! <<< added for turning
 
    omega_fac = 2.d0 - omega
 
    call hecmw_precond_ssor_33_apply_inner( &
    n,np,zp,al,au,d,alu,iteml,itemu,indexl,indexu,perm,ictoblockindex,blockindextocolorindex, &
    indexl(n),indexu(n),ncolor,numofblock,omega_fac)
 
#endif
  end subroutine hecmw_precond_ssor_33_apply_csr
 
  subroutine hecmw_precond_ssor_33_apply_inner( &
    N,NP,ZP,AL,AU,D,ALU,itemL,itemU,indexL,indexU,perm,icToBlockIndex,blockIndexToColorIndex, &
    NPL,NPU,NColor,numOfBlock,OMEGA_FAC)
    integer(kind=kint), intent(in)  :: N
    integer(kind=kint), intent(in)  :: NP
    real(kind=kreal), intent(inout) :: zp(3*np)
    real(kind=kreal), intent(in)    :: al(3*npl)
    real(kind=kreal), intent(in)    :: au(3*npu)
    real(kind=kreal), intent(in)    :: d(9*n)
    real(kind=kreal), intent(in)    :: alu(9*n)
    integer(kind=kint), intent(in)  :: itemL(NPL)
    integer(kind=kint), intent(in)  :: itemU(NPU)
    integer(kind=kint), intent(in)  :: indexL(0:N)
    integer(kind=kint), intent(in)  :: indexU(0:N)
    integer(kind=kint), intent(in)  :: perm(N)
    integer(kind=kint), intent(in)  :: icToBlockIndex(0:NColor)
    integer(kind=kint), intent(in)  :: blockIndexToColorIndex(0:numOfBlock+NColor)
    integer(kind=kint), intent(in)  :: NPL
    integer(kind=kint), intent(in)  :: NPU
    integer(kind=kint), intent(in)  :: NColor
    integer(kind=kint), intent(in)  :: numOfBlock
    real(kind=kreal), intent(in)    :: omega_fac
#ifdef HECMW_ARCH_FX64
 
    integer(kind=kint) :: ic, i, iold, j, isL, ieL, isU, ieU, k
    integer(kind=kint) :: blockIndex
    real(kind=kreal) :: sw1, sw2, sw3, x1, x2, x3
 
    !call start_collection("loopInPrecond33")
 
    !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
    !OCL CACHE_SUBSECTOR_ASSIGN(ZP)
 
    if( ncolor > 1 ) then
 
    !$omp parallel default(none) &
      !$omp&shared(NColor,indexL,itemL,indexU,itemU,AL,AU,D,ALU,perm,&
      !$omp&       ZP,icToBlockIndex,blockIndexToColorIndex,OMEGA_FAC) &
      !$omp&private(SW1,SW2,SW3,X1,X2,X3,ic,i,iold,isL,ieL,isU,ieU,j,k,blockIndex)
 
    !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)
          ! do i = startPos(threadNum, ic), endPos(threadNum, ic)
          iold = perm(i)
          sw1= omega_fac * zp(3*iold-2)
          sw2= omega_fac * zp(3*iold-1)
          sw3= omega_fac * zp(3*iold  )
          isl= indexl(i-1)+1
          iel= indexl(i)
          do j= isl, iel
            !k= perm(itemL(j))
            x1= zp(iteml(j))
            sw1= sw1 - al(3*j-2)*x1
            sw2= sw2 - al(3*j-1)*x1
            sw3= sw3 - al(3*j  )*x1
          enddo ! j
 
          zp(3*iold-2)= alu(9*i-8)*sw1+alu(9*i-7)*sw2+alu(9*i-6)*sw3
          zp(3*iold-1)= alu(9*i-5)*sw1+alu(9*i-4)*sw2+alu(9*i-3)*sw3
          zp(3*iold  )= alu(9*i-2)*sw1+alu(9*i-1)*sw2+alu(9*i  )*sw3
        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
          sw1= 0.d0
          sw2= 0.d0
          sw3= 0.d0
          isu= indexu(i-1) + 1
          ieu= indexu(i)
          do j= ieu, isu, -1
            !k= perm(itemU(j))
            x1= zp(itemu(j))
            sw1= sw1 + au(3*j-2)*x1
            sw2= sw2 + au(3*j-1)*x1
            sw3= sw3 + au(3*j  )*x1
          enddo ! j
 
          iold = perm(i)
          zp(3*iold-2)= zp(3*iold-2) - alu(9*i-8)*sw1 - alu(9*i-7)*sw2 - alu(9*i-6)*sw3
          zp(3*iold-1)= zp(3*iold-1) - alu(9*i-5)*sw1 - alu(9*i-4)*sw2 - alu(9*i-3)*sw3
          zp(3*iold  )= zp(3*iold  ) - alu(9*i-2)*sw1 - alu(9*i-1)*sw2 - alu(9*i  )*sw3
        enddo ! i
      enddo ! blockIndex
      !$omp end do
    enddo ! ic
    !$omp end parallel
 
    else
 
      !C-- FORWARD
      do i = 1, n
        iold = perm(i)
        sw1= omega_fac * zp(3*iold-2)
        sw2= omega_fac * zp(3*iold-1)
        sw3= omega_fac * zp(3*iold  )
        isl= indexl(i-1)+1
        iel= indexl(i)
        do j= isl, iel
          !k= perm(itemL(j))
          x1= zp(iteml(j))
          sw1= sw1 - al(3*j-2)*x1
          sw2= sw2 - al(3*j-1)*x1
          sw3= sw3 - al(3*j  )*x1
        enddo ! j
 
        zp(3*iold-2)= alu(9*i-8)*sw1+alu(9*i-7)*sw2+alu(9*i-6)*sw3
        zp(3*iold-1)= alu(9*i-5)*sw1+alu(9*i-4)*sw2+alu(9*i-3)*sw3
        zp(3*iold  )= alu(9*i-2)*sw1+alu(9*i-1)*sw2+alu(9*i  )*sw3
      enddo ! i
 
      !C-- BACKWARD
      do i = n, 1, -1
        sw1= 0.d0
        sw2= 0.d0
        sw3= 0.d0
        isu= indexu(i-1) + 1
        ieu= indexu(i)
        do j= ieu, isu, -1
          x1= zp(itemu(j))
          sw1= sw1 + au(3*j-2)*x1
          sw2= sw2 + au(3*j-1)*x1
          sw3= sw3 + au(3*j  )*x1
        enddo ! j
 
        iold = perm(i)
        zp(3*iold-2)= zp(3*iold-2) - alu(9*i-8)*sw1 - alu(9*i-7)*sw2 - alu(9*i-6)*sw3
        zp(3*iold-1)= zp(3*iold-1) - alu(9*i-5)*sw1 - alu(9*i-4)*sw2 - alu(9*i-3)*sw3
        zp(3*iold  )= zp(3*iold  ) - alu(9*i-2)*sw1 - alu(9*i-1)*sw2 - alu(9*i  )*sw3
      enddo ! i
 
    end if
 
    !OCL END_CACHE_SUBSECTOR
    !OCL END_CACHE_SECTOR_SIZE
 
    !call stop_collection("loopInPrecond33")
 
#endif
  end subroutine
 
  subroutine hecmw_precond_ssor_33_clear(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    if (associated(colorindex)) deallocate(colorindex)
    if (associated(perm)) deallocate(perm)
    if (associated(iperm)) deallocate(iperm)
    if (associated(alu)) deallocate(alu)
    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)
    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_33_clear
 
end module     hecmw_precond_ssor_33