hecmw_solver_scaling_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
!-------------------------------------------------------------------------------
 
module hecmw_solver_scaling_nn
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
  implicit none
 
  private
  real(kind=kreal), private, allocatable :: scale(:)
 
  public :: hecmw_solver_scaling_fw_nn
  public :: hecmw_solver_scaling_bk_nn
 
contains
 
  subroutine hecmw_solver_scaling_fw_nn(hecMESH, hecMAT, COMMtime)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout) :: hecmat
    real(kind=kreal), intent(inout) :: commtime
    integer(kind=kint) :: n, np, ndof, ndof2
    real(kind=kreal), pointer :: d(:), al(:), au(:), b(:)
    integer(kind=kint), pointer :: inl(:), ial(:), inu(:), iau(:)
    integer(kind=kint) :: i,j,k,ii,ij, ip(hecmat%ndof),iq(hecmat%ndof)
    integer(kind=kint) :: isl, iel, isu, ieu, inod
    real(kind=kreal) :: start_time, end_time
 
    if (hecmw_mat_get_scaling(hecmat).eq.0) return
 
    n = hecmat%N
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
    d => hecmat%D
    al => hecmat%AL
    au => hecmat%AU
    inl => hecmat%indexL
    ial => hecmat%itemL
    inu => hecmat%indexU
    iau => hecmat%itemU
    b => hecmat%B
 
    allocate(scale(ndof*np))
 
    do i= 1, n
      do k=1, ndof
        scale(ndof*(i-1)+k)= 1.d0/dsqrt(dabs(d(ndof*ndof*(i-1)+(k-1)*(ndof+1)+1)))
      end do
    enddo
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, scale, hecmesh%n_node, ndof)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
    do i= 1, np
      do j = 1, ndof
        ip(j)=ndof*(i-1)+j
      end do
      do j = 1, ndof
        do k = 1, ndof
          d(ndof2*(i-1)+ndof*(j-1)+k)=d(ndof2*(i-1)+ndof*(j-1)+k)*scale(ip(j))*scale(ip(k))
        end do
      end do
 
      isl= inl(i-1) + 1
      iel= inl(i  )
      !*voption indep (IAL,AL,SCALE)
      do k= isl, iel
        inod= ial(k)
        do ii = 1, ndof
          iq(ii) = ndof*(inod-1)+ii
        end do
        do ii = 1, ndof
          do ij = 1, ndof
            al(ndof2*(k-1)+ndof*(ii-1)+ij)=al(ndof2*(k-1)+ndof*(ii-1)+ij)*scale(ip(ii))*scale(iq(ij))
          end do
        end do
      enddo
 
      isu= inu(i-1) + 1
      ieu= inu(i  )
      !*voption indep (IAU,AU,SCALE)
      do k= isu, ieu
        inod= iau(k)
        do ii = 1, ndof
          iq(ii) = ndof*(inod-1)+ii
        end do
        do ii = 1, ndof
          do ij = 1, ndof
            au(ndof2*(k-1)+ndof*(ii-1)+ij)=au(ndof2*(k-1)+ndof*(ii-1)+ij)*scale(ip(ii))*scale(iq(ij))
          end do
        end do
      enddo
    enddo
    !*voption indep (B,SCALE)
    do i= 1, n
      do k = 1, ndof
        b(ndof*(i-1)+k)=b(ndof*(i-1)+k)*scale(ndof*(i-1)+k)
      end do
    enddo
  end subroutine hecmw_solver_scaling_fw_nn
 
  subroutine hecmw_solver_scaling_bk_nn(hecMAT)
    use hecmw_util
    implicit none
    type (hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint) :: n, np, ndof, ndof2
    real(kind=kreal), pointer :: d(:), al(:), au(:), b(:), x(:)
    integer(kind=kint), pointer :: inl(:), ial(:), inu(:), iau(:)
    integer(kind=kint) :: i,j,k,ii,ij, ip(hecmat%ndof),iq(hecmat%ndof)
    integer(kind=kint) :: isl, iel, isu, ieu, inod
 
    if (hecmw_mat_get_scaling(hecmat).eq.0) return
 
    n = hecmat%N
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
    d => hecmat%D
    al => hecmat%AL
    au => hecmat%AU
    inl => hecmat%indexL
    ial => hecmat%itemL
    inu => hecmat%indexU
    iau => hecmat%itemU
    b => hecmat%B
    x => hecmat%X
 
    !*voption indep (X,B,SCALE)
    do i= 1, n
      do k=1,ndof
        x(ndof*(i-1)+k)=x(ndof*(i-1)+k)*scale(ndof*(i-1)+k)
        b(ndof*(i-1)+k)=b(ndof*(i-1)+k)/scale(ndof*(i-1)+k)
      end do
    enddo
 
    do i= 1, np
      do j = 1, ndof
        ip(j)=ndof*(i-1)+j
      end do
      do j = 1, ndof
        do k = 1, ndof
          d(ndof2*(i-1)+ndof*(j-1)+k)=d(ndof2*(i-1)+ndof*(j-1)+k)/(scale(ip(j))*scale(ip(k)))
        end do
      end do
 
      isl= inl(i-1) + 1
      iel= inl(i  )
      !*voption indep (IAL,AL,SCALE)
      do k= isl, iel
        inod= ial(k)
        do ii = 1, ndof
          iq(ii) = ndof*(inod-1)+ii
        end do
        do ii = 1, ndof
          do ij = 1, ndof
            al(ndof2*(k-1)+ndof*(ii-1)+ij)=al(ndof2*(k-1)+ndof*(ii-1)+ij)/(scale(ip(ii))*scale(iq(ij)))
          end do
        end do
      enddo
 
      isu= inu(i-1) + 1
      ieu= inu(i  )
      !*voption indep (IAU,AU,SCALE)
      do k= isu, ieu
        inod= iau(k)
        do ii = 1, ndof
          iq(ii) = ndof*(inod-1)+ii
        end do
        do ii = 1, ndof
          do ij = 1, ndof
            au(ndof2*(k-1)+ndof*(ii-1)+ij)=au(ndof2*(k-1)+ndof*(ii-1)+ij)/(scale(ip(ii))*scale(iq(ij)))
          end do
        end do
      enddo
    enddo
 
    deallocate(scale)
  end subroutine hecmw_solver_scaling_bk_nn
 
end module hecmw_solver_scaling_nn