hecmw_precond_SSOR_66.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_66
!C***
!C
module hecmw_precond_ssor_66
  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_66_setup
  public:: hecmw_precond_ssor_66_apply
  public:: hecmw_precond_ssor_66_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.
 
contains
 
  subroutine hecmw_precond_ssor_66_setup(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(in) :: hecmat
    integer(kind=kint ) :: npl, npu
    integer(kind=kint ) :: ncolor_in
    real   (kind=kreal) :: sigma_diag
    real   (kind=kreal) :: alutmp(6,6), pw(6)
    integer(kind=kint ) :: ii, i, j, k
    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
 
#ifndef _OPENACC
    !$ nthreads = omp_get_max_threads()
#endif
 
    n = hecmat%N
    ! N = hecMAT%NP
    ncolor_in = hecmw_mat_get_ncolor_in(hecmat)
    sigma_diag = hecmw_mat_get_sigma_diag(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)
    !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)
 
 
    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(36*n), al(36*npl), au(36*npu))
    call hecmw_matrix_reorder_values(n, 6, 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)
#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
 
      d => hecmat%D
      al => hecmat%AL
      au => hecmat%AU
      indexl => hecmat%indexL
      indexu => hecmat%indexU
      iteml => hecmat%itemL
      itemu => hecmat%itemU
    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)
 
 
      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(36*n), al(36*npl), au(36*npu))
      call hecmw_matrix_reorder_values(n, 6, 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)
    end if
#endif
 
    allocate(alu(36*n))
    alu  = 0.d0
 
    do ii= 1, 36*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,SIGMA_DIAG)
    !$omp do
#endif
    do ii= 1, n
      alutmp(1,1)= alu(36*ii-35) * sigma_diag
      alutmp(1,2)= alu(36*ii-34)
      alutmp(1,3)= alu(36*ii-33)
      alutmp(1,4)= alu(36*ii-32)
      alutmp(1,5)= alu(36*ii-31)
      alutmp(1,6)= alu(36*ii-30)
 
      alutmp(2,1)= alu(36*ii-29)
      alutmp(2,2)= alu(36*ii-28) * sigma_diag
      alutmp(2,3)= alu(36*ii-27)
      alutmp(2,4)= alu(36*ii-26)
      alutmp(2,5)= alu(36*ii-25)
      alutmp(2,6)= alu(36*ii-24)
 
      alutmp(3,1)= alu(36*ii-23)
      alutmp(3,2)= alu(36*ii-22)
      alutmp(3,3)= alu(36*ii-21) * sigma_diag
      alutmp(3,4)= alu(36*ii-20)
      alutmp(3,5)= alu(36*ii-19)
      alutmp(3,6)= alu(36*ii-18)
 
      alutmp(4,1)= alu(36*ii-17)
      alutmp(4,2)= alu(36*ii-16)
      alutmp(4,3)= alu(36*ii-15)
      alutmp(4,4)= alu(36*ii-14) * sigma_diag
      alutmp(4,5)= alu(36*ii-13)
      alutmp(4,6)= alu(36*ii-12)
 
      alutmp(5,1)= alu(36*ii-11)
      alutmp(5,2)= alu(36*ii-10)
      alutmp(5,3)= alu(36*ii-9 )
      alutmp(5,4)= alu(36*ii-8 )
      alutmp(5,5)= alu(36*ii-7 ) * sigma_diag
      alutmp(5,6)= alu(36*ii-6 )
 
      alutmp(6,1)= alu(36*ii-5 )
      alutmp(6,2)= alu(36*ii-4 )
      alutmp(6,3)= alu(36*ii-3 )
      alutmp(6,4)= alu(36*ii-2 )
      alutmp(6,5)= alu(36*ii-1 )
      alutmp(6,6)= alu(36*ii   ) * sigma_diag
 
      do k= 1, 6
        alutmp(k,k)= 1.d0/alutmp(k,k)
        do i= k+1, 6
          alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
          do j= k+1, 6
            pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
          enddo
          do j= k+1, 6
            alutmp(i,j)= pw(j)
          enddo
        enddo
      enddo
 
      alu(36*ii-35)= alutmp(1,1)
      alu(36*ii-34)= alutmp(1,2)
      alu(36*ii-33)= alutmp(1,3)
      alu(36*ii-32)= alutmp(1,4)
      alu(36*ii-31)= alutmp(1,5)
      alu(36*ii-30)= alutmp(1,6)
      alu(36*ii-29)= alutmp(2,1)
      alu(36*ii-28)= alutmp(2,2)
      alu(36*ii-27)= alutmp(2,3)
      alu(36*ii-26)= alutmp(2,4)
      alu(36*ii-25)= alutmp(2,5)
      alu(36*ii-24)= alutmp(2,6)
      alu(36*ii-23)= alutmp(3,1)
      alu(36*ii-22)= alutmp(3,2)
      alu(36*ii-21)= alutmp(3,3)
      alu(36*ii-20)= alutmp(3,4)
      alu(36*ii-19)= alutmp(3,5)
      alu(36*ii-18)= alutmp(3,6)
      alu(36*ii-17)= alutmp(4,1)
      alu(36*ii-16)= alutmp(4,2)
      alu(36*ii-15)= alutmp(4,3)
      alu(36*ii-14)= alutmp(4,4)
      alu(36*ii-13)= alutmp(4,5)
      alu(36*ii-12)= alutmp(4,6)
      alu(36*ii-11)= alutmp(5,1)
      alu(36*ii-10)= alutmp(5,2)
      alu(36*ii-9 )= alutmp(5,3)
      alu(36*ii-8 )= alutmp(5,4)
      alu(36*ii-7 )= alutmp(5,5)
      alu(36*ii-6 )= alutmp(5,6)
      alu(36*ii-5 )= alutmp(6,1)
      alu(36*ii-4 )= alutmp(6,2)
      alu(36*ii-3 )= alutmp(6,3)
      alu(36*ii-2 )= alutmp(6,4)
      alu(36*ii-1 )= alutmp(6,5)
      alu(36*ii   )= alutmp(6,6)
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
 
    isfirst = .true.
 
    !write(*,*) 'DEBUG: SSOR setup done', hecmw_Wtime()-t0
 
  end subroutine hecmw_precond_ssor_66_setup
 
  subroutine hecmw_precond_ssor_66_apply(ZP)
    use hecmw_tuning_fx
    implicit none
    real(kind=kreal), intent(inout) :: zp(:)
    integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k
    real(kind=kreal) :: x1, x2, x3, x4, x5, x6
    real(kind=kreal) :: sw1, sw2, sw3, sw4, sw5, sw6
 
    ! added 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) :: blockindex, elementcount, numofelement, ii
    real(kind=kreal) :: numofelementperblock
    integer(kind=kint) :: my_rank
 
#ifndef _OPENACC
    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, 3, &
        sectorcachesize0, sectorcachesize1 )
 
      isfirst = .false.
    endif
#endif
    ! <<< added for tuning
 
#ifndef _OPENACC
    !call start_collection("loopInPrecond66")
 
    !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) &
      !$omp&private(SW1,SW2,SW3,SW4,SW5,SW6,X1,X2,X3,X4,X5,X6,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= zp(6*iold-5)
          sw2= zp(6*iold-4)
          sw3= zp(6*iold-3)
          sw4= zp(6*iold-2)
          sw5= zp(6*iold-1)
          sw6= zp(6*iold  )
          isl= indexl(i-1)+1
          iel= indexl(i)
          do j= isl, iel
            !k= perm(itemL(j))
            k= iteml(j)
            x1= zp(6*k-5)
            x2= zp(6*k-4)
            x3= zp(6*k-3)
            x4= zp(6*k-2)
            x5= zp(6*k-1)
            x6= zp(6*k  )
            sw1= sw1 -al(36*j-35)*x1 -al(36*j-34)*x2 -al(36*j-33)*x3 -al(36*j-32)*x4 -al(36*j-31)*x5 -al(36*j-30)*x6
            sw2= sw2 -al(36*j-29)*x1 -al(36*j-28)*x2 -al(36*j-27)*x3 -al(36*j-26)*x4 -al(36*j-25)*x5 -al(36*j-24)*x6
            sw3= sw3 -al(36*j-23)*x1 -al(36*j-22)*x2 -al(36*j-21)*x3 -al(36*j-20)*x4 -al(36*j-19)*x5 -al(36*j-18)*x6
            sw4= sw4 -al(36*j-17)*x1 -al(36*j-16)*x2 -al(36*j-15)*x3 -al(36*j-14)*x4 -al(36*j-13)*x5 -al(36*j-12)*x6
            sw5= sw5 -al(36*j-11)*x1 -al(36*j-10)*x2 -al(36*j-9 )*x3 -al(36*j-8 )*x4 -al(36*j-7 )*x5 -al(36*j-6 )*x6
            sw6= sw6 -al(36*j-5 )*x1 -al(36*j-4 )*x2 -al(36*j-3 )*x3 -al(36*j-2 )*x4 -al(36*j-1 )*x5 -al(36*j   )*x6
          enddo ! j
 
          x1= sw1
          x2= sw2
          x3= sw3
          x4= sw4
          x5= sw5
          x6= sw6
          x2= x2 -alu(36*i-29)*x1
          x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
          x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-5)*x3
          x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9)*x3 -alu(36*i-8)*x4
          x6= x6 -alu(36*i-5 )*x1 -alu(36*i-4 )*x2 -alu(36*i-3)*x3 -alu(36*i-2)*x4 -alu(36*i-1)*x5
          x6= alu(36*i   )*  x6
          x5= alu(36*i-7 )*( x5 -alu(36*i-6)*x6 )
          x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
          x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
          x2= alu(36*i-28)*( x2 -alu(36*i-24)*x6 -alu(36*i-25)*x5 -alu(36*i-26)*x4 -alu(36*i-27)*x3)
          x1= alu(36*i-35)*( x1 -alu(36*i-30)*x6 -alu(36*i-31)*x5 -alu(36*i-32)*x4 -alu(36*i-33)*x3 -alu(36*i-34)*x2)
          zp(6*iold-5)= x1
          zp(6*iold-4)= x2
          zp(6*iold-3)= x3
          zp(6*iold-2)= x4
          zp(6*iold-1)= x5
          zp(6*iold  )= x6
#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
          sw4= 0.d0
          sw5= 0.d0
          sw6= 0.d0
          isu= indexu(i-1) + 1
          ieu= indexu(i)
          do j= ieu, isu, -1
            !k= perm(itemU(j))
            k= itemu(j)
            x1= zp(6*k-5)
            x2= zp(6*k-4)
            x3= zp(6*k-3)
            x4= zp(6*k-2)
            x5= zp(6*k-1)
            x6= zp(6*k  )
            sw1= sw1 +au(36*j-35)*x1 +au(36*j-34)*x2 +au(36*j-33)*x3 +au(36*j-32)*x4 +au(36*j-31)*x5 +au(36*j-30)*x6
            sw2= sw2 +au(36*j-29)*x1 +au(36*j-28)*x2 +au(36*j-27)*x3 +au(36*j-26)*x4 +au(36*j-25)*x5 +au(36*j-24)*x6
            sw3= sw3 +au(36*j-23)*x1 +au(36*j-22)*x2 +au(36*j-21)*x3 +au(36*j-20)*x4 +au(36*j-19)*x5 +au(36*j-18)*x6
            sw4= sw4 +au(36*j-17)*x1 +au(36*j-16)*x2 +au(36*j-15)*x3 +au(36*j-14)*x4 +au(36*j-13)*x5 +au(36*j-12)*x6
            sw5= sw5 +au(36*j-11)*x1 +au(36*j-10)*x2 +au(36*j-9 )*x3 +au(36*j-8 )*x4 +au(36*j-7 )*x5 +au(36*j-6 )*x6
            sw6= sw6 +au(36*j-5 )*x1 +au(36*j-4 )*x2 +au(36*j-3 )*x3 +au(36*j-2 )*x4 +au(36*j-1 )*x5 +au(36*j   )*x6
          enddo ! j
 
          x1= sw1
          x2= sw2
          x3= sw3
          x4= sw4
          x5= sw5
          x6= sw6
          x2= x2 -alu(36*i-29)*x1
          x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
          x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-5)*x3
          x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9)*x3 -alu(36*i-8)*x4
          x6= x6 -alu(36*i-5 )*x1 -alu(36*i-4 )*x2 -alu(36*i-3)*x3 -alu(36*i-2)*x4 -alu(36*i-1)*x5
          x6= alu(36*i   )*  x6
          x5= alu(36*i-7 )*( x5 -alu(36*i-6)*x6 )
          x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
          x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
          x2= alu(36*i-28)*( x2 -alu(36*i-24)*x6 -alu(36*i-25)*x5 -alu(36*i-26)*x4 -alu(36*i-27)*x3)
          x1= alu(36*i-35)*( x1 -alu(36*i-30)*x6 -alu(36*i-31)*x5 -alu(36*i-32)*x4 -alu(36*i-33)*x3 -alu(36*i-34)*x2)
          iold = perm(i)
          zp(6*iold-5)= zp(6*iold-5) -x1
          zp(6*iold-4)= zp(6*iold-4) -x2
          zp(6*iold-3)= zp(6*iold-3) -x3
          zp(6*iold-2)= zp(6*iold-2) -x4
          zp(6*iold-1)= zp(6*iold-1) -x5
          zp(6*iold  )= zp(6*iold  ) -x6
#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("loopInPrecond66")
#endif
 
  end subroutine hecmw_precond_ssor_66_apply
 
  subroutine hecmw_precond_ssor_66_clear(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: nthreads = 1
#ifndef _OPENACC
    !$ nthreads = omp_get_max_threads()
#endif
    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)
  end subroutine hecmw_precond_ssor_66_clear
 
 
  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_66