hecmw_precond_DIAG_44.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_DIAG_44
!C***
!C
module hecmw_precond_diag_44
  use  hecmw_util
 
  private
 
  public:: hecmw_precond_diag_44_setup
  public:: hecmw_precond_diag_44_apply
  public:: hecmw_precond_diag_44_clear
 
  integer(kind=kint) :: N
  real(kind=kreal), pointer :: alu(:) => null()
 
  logical, save :: INITIALIZED = .false.
 
contains
 
  subroutine hecmw_precond_diag_44_setup(hecMAT)
    use hecmw_matrix_misc
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: np
    real   (kind=kreal) :: sigma_diag
    real(kind=kreal), pointer:: d(:)
 
    real   (kind=kreal):: alutmp(4,4), pw(4)
    integer(kind=kint ):: ii, i, j, k
 
    if (initialized) then
      if (hecmat%Iarray(98) == 0 .and. hecmat%Iarray(97) == 0) return
      call hecmw_precond_diag_44_clear()
    endif
 
    n = hecmat%N
    np = hecmat%NP
    d => hecmat%D
    sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
 
    allocate(alu(16*np))
    alu = 0.d0
 
    do ii= 1, n
      alu(16*ii-15) = d(16*ii-15)
      alu(16*ii-14) = d(16*ii-14)
      alu(16*ii-13) = d(16*ii-13)
      alu(16*ii-12) = d(16*ii-12)
      alu(16*ii-11) = d(16*ii-11)
      alu(16*ii-10) = d(16*ii-10)
      alu(16*ii-9 ) = d(16*ii-9 )
      alu(16*ii-8 ) = d(16*ii-8 )
      alu(16*ii-7 ) = d(16*ii-7 )
      alu(16*ii-6 ) = d(16*ii-6 )
      alu(16*ii-5 ) = d(16*ii-5 )
      alu(16*ii-4 ) = d(16*ii-4 )
      alu(16*ii-3 ) = d(16*ii-3 )
      alu(16*ii-2 ) = d(16*ii-2 )
      alu(16*ii-1 ) = d(16*ii-1 )
      alu(16*ii   ) = d(16*ii   )
    enddo
 
    !$omp parallel default(none),private(ii,ALUtmp,k,i,j,PW),shared(N,ALU,SIGMA_DIAG)
    !$omp do
    do ii= 1, n
      alutmp(1,1)= alu(16*ii-15) * sigma_diag
      alutmp(1,2)= alu(16*ii-14)
      alutmp(1,3)= alu(16*ii-13)
      alutmp(1,4)= alu(16*ii-12)
 
      alutmp(2,1)= alu(16*ii-11)
      alutmp(2,2)= alu(16*ii-10) * sigma_diag
      alutmp(2,3)= alu(16*ii- 9)
      alutmp(2,4)= alu(16*ii- 8)
 
      alutmp(3,1)= alu(16*ii- 7)
      alutmp(3,2)= alu(16*ii- 6)
      alutmp(3,3)= alu(16*ii- 5) * sigma_diag
      alutmp(3,4)= alu(16*ii- 4)
 
      alutmp(4,1)= alu(16*ii- 3)
      alutmp(4,2)= alu(16*ii- 2)
      alutmp(4,3)= alu(16*ii- 1)
      alutmp(4,4)= alu(16*ii   ) * sigma_diag
 
      do k= 1, 4
        alutmp(k,k)= 1.d0/alutmp(k,k)
        do i= k+1, 4
          alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
          do j= k+1, 4
            pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
          enddo
          do j= k+1, 4
            alutmp(i,j)= pw(j)
          enddo
        enddo
      enddo
 
      alu(16*ii-15)= alutmp(1,1)
      alu(16*ii-14)= alutmp(1,2)
      alu(16*ii-13)= alutmp(1,3)
      alu(16*ii-12)= alutmp(1,4)
      alu(16*ii-11)= alutmp(2,1)
      alu(16*ii-10)= alutmp(2,2)
      alu(16*ii- 9)= alutmp(2,3)
      alu(16*ii- 8)= alutmp(2,4)
      alu(16*ii- 7)= alutmp(3,1)
      alu(16*ii- 6)= alutmp(3,2)
      alu(16*ii- 5)= alutmp(3,3)
      alu(16*ii- 4)= alutmp(3,4)
      alu(16*ii- 3)= alutmp(4,1)
      alu(16*ii- 2)= alutmp(4,2)
      alu(16*ii- 1)= alutmp(4,3)
      alu(16*ii   )= alutmp(4,4)
    enddo
    !$omp end do
    !$omp end parallel
 
    initialized = .true.
    hecmat%Iarray(98) = 0 ! symbolic setup done
    hecmat%Iarray(97) = 0 ! numerical setup done
  end subroutine hecmw_precond_diag_44_setup
 
  subroutine hecmw_precond_diag_44_apply(WW)
    implicit none
    real(kind=kreal), intent(inout) :: ww(:)
    integer(kind=kint) :: i
    real(kind=kreal) :: x1, x2, x3, x4
 
    !C
    !C== Block SCALING
    !$omp parallel default(none),private(i,X1,X2,X3,X4),shared(N,WW,ALU)
    !$omp do
    do i= 1, n
      x1= ww(4*i-3)
      x2= ww(4*i-2)
      x3= ww(4*i-1)
      x4= ww(4*i  )
      x2= x2 -alu(16*i-11)*x1
      x3= x3 -alu(16*i- 7)*x1 -alu(16*i- 6)*x2
      x4= x4 -alu(16*i- 3)*x1 -alu(16*i- 2)*x2 -alu(16*i- 1)*x3
      x4= alu(16*i   )* x4
      x3= alu(16*i- 5)*(x3 -alu(16*i- 4)*x4)
      x2= alu(16*i-10)*(x2 -alu(16*i- 8)*x4 -alu(16*i- 9)*x3)
      x1= alu(16*i-15)*(x1 -alu(16*i-12)*x4 -alu(16*i-13)*x3 -alu(16*i-14)*x2)
      ww(4*i-3)= x1
      ww(4*i-2)= x2
      ww(4*i-1)= x3
      ww(4*i  )= x4
    enddo
    !$omp end do
    !$omp end parallel
 
  end subroutine hecmw_precond_diag_44_apply
 
  subroutine hecmw_precond_diag_44_clear()
    implicit none
    if (associated(alu)) deallocate(alu)
    nullify(alu)
    initialized = .false.
  end subroutine hecmw_precond_diag_44_clear
 
end module     hecmw_precond_diag_44