hecmw_precond_DIAG_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_DIAG_nn
!C***
!C
module hecmw_precond_diag_nn
  use  hecmw_util
 
  private
 
  public:: hecmw_precond_diag_nn_setup
  public:: hecmw_precond_diag_nn_apply
  public:: hecmw_precond_diag_nn_clear
 
  integer(kind=kint) :: N
  real(kind=kreal), pointer :: alu(:) => null()
 
  logical, save :: INITIALIZED = .false.
 
contains
 
  subroutine hecmw_precond_diag_nn_setup(hecMAT)
    use hecmw_matrix_misc
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: np,ndof,ndof2
    real   (kind=kreal) :: sigma_diag
    real(kind=kreal), pointer:: d(:)
 
    real   (kind=kreal):: alutmp(hecmat%NDOF,hecmat%NDOF), pw(hecmat%NDOF)
    integer(kind=kint ):: i, j, k, ii
 
    if (initialized) then
      if (hecmat%Iarray(98) == 0 .and. hecmat%Iarray(97) == 0) return
      call hecmw_precond_diag_nn_clear()
    endif
 
    n = hecmat%N
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
    np = hecmat%NP
    d => hecmat%D
    sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
 
    allocate(alu(ndof2*np))
    alu = 0.d0
 
    do ii= 1, n
      do i = 1, ndof2
        alu(ndof2*(ii-1)+i)=d(ndof2*(ii-1)+i)
      end do
    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
    initialized = .true.
    hecmat%Iarray(98) = 0 ! symbolic setup done
    hecmat%Iarray(97) = 0 ! numerical setup done
 
  end subroutine hecmw_precond_diag_nn_setup
 
  subroutine hecmw_precond_diag_nn_apply(WW,NDOF)
    implicit none
    real(kind=kreal), intent(inout) :: ww(:)
    integer(kind=kint) :: i,j,k,ndof
    real(kind=kreal) :: x(ndof)
 
    !C
    !C== Block SCALING
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#else
    !$omp parallel default(none),private(i,j,k,X),shared(N,WW,ALU,NDOF)
    !$omp do
#endif
    do i= 1, n
      do j=1,ndof
        x(j)=ww(ndof*(i-1)+j)
      end do
      do j=2,ndof
        do k = 1,j-1
          x(j)=x(j)-alu(ndof*ndof*(i-1)+ndof*(j-1)+k )*x(k)
        end do
      end do
      do j=ndof,1,-1
        do k = ndof,j+1,-1
          x(j)=x(j)-alu(ndof*ndof*(i-1)+ndof*(j-1)+k )*x(k)
        end do
        x(j)=alu(ndof*ndof*(i-1)+(ndof+1)*(j-1)+1 )*x(j)
      end do
      do j=1,ndof
        ww(ndof*(i-1)+j)=x(j)
      end do
 
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
 
  end subroutine hecmw_precond_diag_nn_apply
 
  subroutine hecmw_precond_diag_nn_clear()
    implicit none
    if (associated(alu)) deallocate(alu)
    nullify(alu)
    initialized = .false.
  end subroutine hecmw_precond_diag_nn_clear
 
end module     hecmw_precond_diag_nn