hecmw_precond_BILU_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_BILU_44
!C***
!C
module hecmw_precond_bilu_44
  use hecmw_util
  use hecmw_matrix_misc
 
  private
 
  public:: hecmw_precond_bilu_44_setup
  public:: hecmw_precond_bilu_44_apply
  public:: hecmw_precond_bilu_44_clear
 
  integer(kind=kint) :: N
  real(kind=kreal), pointer :: dlu0(:) => null()
  real(kind=kreal), pointer :: allu0(:) => null()
  real(kind=kreal), pointer :: aulu0(:) => null()
  integer(kind=kint), pointer :: inumFI1L(:) => null()
  integer(kind=kint), pointer :: inumFI1U(:) => null()
  integer(kind=kint), pointer :: FI1L(:) => null()
  integer(kind=kint), pointer :: FI1U(:) => null()
 
  logical, save :: INITIALIZED = .false.
 
contains
 
  subroutine hecmw_precond_bilu_44_setup(hecMAT)
    implicit none
    type(hecmwst_matrix), intent(inout) :: hecmat
    integer(kind=kint ) :: np, npu, npl
    integer(kind=kint ) :: precond
    real   (kind=kreal) :: sigma, sigma_diag
 
    real(kind=kreal), pointer :: d(:)
    real(kind=kreal), pointer :: al(:)
    real(kind=kreal), pointer :: au(:)
 
    integer(kind=kint ), pointer :: inl(:), inu(:)
    integer(kind=kint ), pointer :: ial(:)
    integer(kind=kint ), pointer :: iau(:)
 
    if (initialized) then
      if (hecmat%Iarray(98) == 1) then ! need symbolic and numerical setup
        call hecmw_precond_bilu_44_clear()
      else if (hecmat%Iarray(97) == 1) then ! need numerical setup only
        call hecmw_precond_bilu_44_clear() ! TEMPORARY
      else
        return
      endif
    endif
 
    n = hecmat%N
    np = hecmat%NP
    npl = hecmat%NPL
    npu = hecmat%NPU
    d => hecmat%D
    al => hecmat%AL
    au => hecmat%AU
    inl => hecmat%indexL
    inu => hecmat%indexU
    ial => hecmat%itemL
    iau => hecmat%itemU
    precond = hecmw_mat_get_precond(hecmat)
    sigma = hecmw_mat_get_sigma(hecmat)
    sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
 
    !if (PRECOND.eq.10) call FORM_ILU0_44 &
      call form_ilu0_44 &
      &   (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
      &    sigma, sigma_diag)
 
    initialized = .true.
    hecmat%Iarray(98) = 0 ! symbolic setup done
    hecmat%Iarray(97) = 0 ! numerical setup done
 
  end subroutine hecmw_precond_bilu_44_setup
 
  subroutine hecmw_precond_bilu_44_apply(WW)
    implicit none
    real(kind=kreal), intent(inout) :: ww(:)
    integer(kind=kint) :: i, j, isl, iel, isu, ieu, k
    real(kind=kreal) :: sw1, sw2, sw3, sw4, x1, x2, x3, x4
    !C
    !C-- FORWARD
 
    do i= 1, n
      sw1= ww(4*i-3)
      sw2= ww(4*i-2)
      sw3= ww(4*i-1)
      sw4= ww(4*i  )
      isl= inumfi1l(i-1)+1
      iel= inumfi1l(i)
      do j= isl, iel
        k= fi1l(j)
        x1= ww(4*k-3)
        x2= ww(4*k-2)
        x3= ww(4*k-1)
        x4= ww(4*k  )
        sw1= sw1 - allu0(16*j-15)*x1-allu0(16*j-14)*x2-allu0(16*j-13)*x3-allu0(16*j-12)*x4
        sw2= sw2 - allu0(16*j-11)*x1-allu0(16*j-10)*x2-allu0(16*j- 9)*x3-allu0(16*j- 8)*x4
        sw3= sw3 - allu0(16*j- 7)*x1-allu0(16*j- 6)*x2-allu0(16*j- 5)*x3-allu0(16*j- 4)*x4
        sw4= sw4 - allu0(16*j- 3)*x1-allu0(16*j- 2)*x2-allu0(16*j- 1)*x3-allu0(16*j   )*x4
      enddo
 
      x1= sw1
      x2= sw2
      x3= sw3
      x4= sw4
      x2= x2 - dlu0(16*i-11)*x1
      x3= x3 - dlu0(16*i- 7)*x1 - dlu0(16*i-6)*x2
      x4= x4 - dlu0(16*i- 3)*x1 - dlu0(16*i-2)*x2 - dlu0(16*i-1)*x3
      x4= dlu0(16*i   )* x4
      x3= dlu0(16*i- 5)*(x3 - dlu0(16*i- 4)*x4)
      x2= dlu0(16*i-10)*(x2 - dlu0(16*i- 8)*x4 - dlu0(16*i- 9)*x3 )
      x1= dlu0(16*i-15)*(x1 - dlu0(16*i-12)*x4 - dlu0(16*i-13)*x3 - dlu0(16*i-14)*x2)
 
      ww(4*i-3)= x1
      ww(4*i-2)= x2
      ww(4*i-1)= x3
      ww(4*i  )= x4
    enddo
 
    !C
    !C-- BACKWARD
 
    do i= n, 1, -1
      isu= inumfi1u(i-1) + 1
      ieu= inumfi1u(i)
      sw1= 0.d0
      sw2= 0.d0
      sw3= 0.d0
      sw4= 0.d0
      do j= ieu, isu, -1
        k= fi1u(j)
        x1= ww(4*k-3)
        x2= ww(4*k-2)
        x3= ww(4*k-1)
        x4= ww(4*k  )
        sw1= sw1 + aulu0(16*j-15)*x1+aulu0(16*j-14)*x2+aulu0(16*j-13)*x3+aulu0(16*j-12)*x4
        sw2= sw2 + aulu0(16*j-11)*x1+aulu0(16*j-10)*x2+aulu0(16*j- 9)*x3+aulu0(16*j- 8)*x4
        sw3= sw3 + aulu0(16*j- 7)*x1+aulu0(16*j- 6)*x2+aulu0(16*j- 5)*x3+aulu0(16*j- 4)*x4
        sw4= sw4 + aulu0(16*j- 3)*x1+aulu0(16*j- 2)*x2+aulu0(16*j- 1)*x3+aulu0(16*j   )*x4
      enddo
      x1= sw1
      x2= sw2
      x3= sw3
      x4= sw4
      x2= x2 - dlu0(16*i-11)*x1
      x3= x3 - dlu0(16*i- 7)*x1 - dlu0(16*i-6)*x2
      x4= x4 - dlu0(16*i- 3)*x1 - dlu0(16*i-2)*x2 - dlu0(16*i-1)*x3
      x4= dlu0(16*i   )*  x4
      x3= dlu0(16*i- 5)*( x3 - dlu0(16*i- 4)*x4 )
      x2= dlu0(16*i-10)*( x2 - dlu0(16*i- 8)*x4 - dlu0(16*i- 9)*x3 )
      x1= dlu0(16*i-15)*( x1 - dlu0(16*i-12)*x4 - dlu0(16*i-13)*x3 - dlu0(16*i-14)*x2)
      ww(4*i-3)=  ww(4*i-3) - x1
      ww(4*i-2)=  ww(4*i-2) - x2
      ww(4*i-1)=  ww(4*i-1) - x3
      ww(4*i  )=  ww(4*i  ) - x4
    enddo
  end subroutine hecmw_precond_bilu_44_apply
 
  subroutine hecmw_precond_bilu_44_clear()
    implicit none
    if (associated(dlu0)) deallocate(dlu0)
    if (associated(allu0)) deallocate(allu0)
    if (associated(aulu0)) deallocate(aulu0)
    if (associated(inumfi1l)) deallocate(inumfi1l)
    if (associated(inumfi1u)) deallocate(inumfi1u)
    if (associated(fi1l)) deallocate(fi1l)
    if (associated(fi1u)) deallocate(fi1u)
    nullify(dlu0)
    nullify(allu0)
    nullify(aulu0)
    nullify(inumfi1l)
    nullify(inumfi1u)
    nullify(fi1l)
    nullify(fi1u)
    initialized = .false.
  end subroutine hecmw_precond_bilu_44_clear
 
  !C
  !C***
  !C*** FORM_ILU0_44
  !C***
  !C
  !C    form ILU(0) matrix
  !C
  subroutine form_ilu0_44                                   &
      &   (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
      &    sigma, sigma_diag)
    implicit none
    integer(kind=kint ), intent(in):: n, np, npu, npl
    real   (kind=kreal), intent(in):: sigma, sigma_diag
 
    real(kind=kreal), dimension(16*NPL), intent(in):: al
    real(kind=kreal), dimension(16*NPU), intent(in):: au
    real(kind=kreal), dimension(16*NP ), intent(in):: d
 
    integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
    integer(kind=kint ), dimension(  NPL),intent(in) :: ial
    integer(kind=kint ), dimension(  NPU),intent(in) :: iau
 
    integer(kind=kint), dimension(:), allocatable :: iw1, iw2
    real (kind=kreal),  dimension(4,4) :: rhs_aij, dkinv, aik, akj
    integer(kind=kint) :: i,jj,jj1,ij0,kk,kk1
    integer(kind=kint) :: j,k
    allocate (iw1(np) , iw2(np))
    allocate(dlu0(9*np), allu0(9*npl), aulu0(9*npu))
    allocate(inumfi1l(0:np), inumfi1u(0:np), fi1l(npl), fi1u(npu))
 
    do i=1,16*np
      dlu0(i) = d(i)
    end do
    do i=1,16*npl
      allu0(i) = al(i)
    end do
    do i=1,16*npu
      aulu0(i) = au(i)
    end do
    do i=0,np
      inumfi1l(i) = inl(i)
      inumfi1u(i) = inu(i)
    end do
    do i=1,npl
      fi1l(i) = ial(i)
    end do
    do i=1,npu
      fi1u(i) = iau(i)
    end do
 
    !C
    !C +----------------------+
    !C | ILU(0) factorization |
    !C +----------------------+
    !C===
    do i=1,np
      dlu0(16*i-15)=dlu0(16*i-15)*sigma_diag
      dlu0(16*i-10)=dlu0(16*i-10)*sigma_diag
      dlu0(16*i- 5)=dlu0(16*i- 5)*sigma_diag
      dlu0(16*i   )=dlu0(16*i   )*sigma_diag
    enddo
 
    i = 1
    call ilu1a44 (dkinv, &
      dlu0(16*i-15), dlu0(16*i-14), dlu0(16*i-13), dlu0(16*i-12), &
      dlu0(16*i-11), dlu0(16*i-10), dlu0(16*i- 9), dlu0(16*i- 8), &
      dlu0(16*i- 7), dlu0(16*i- 6), dlu0(16*i- 5), dlu0(16*i- 4), &
      dlu0(16*i- 3), dlu0(16*i- 2), dlu0(16*i- 1), dlu0(16*i   ) )
    dlu0(16*i-15)= dkinv(1,1)
    dlu0(16*i-14)= dkinv(1,2)
    dlu0(16*i-13)= dkinv(1,3)
    dlu0(16*i-12)= dkinv(1,4)
    dlu0(16*i-11)= dkinv(2,1)
    dlu0(16*i-10)= dkinv(2,2)
    dlu0(16*i- 9)= dkinv(2,3)
    dlu0(16*i- 8)= dkinv(2,4)
    dlu0(16*i- 7)= dkinv(3,1)
    dlu0(16*i- 6)= dkinv(3,2)
    dlu0(16*i- 5)= dkinv(3,3)
    dlu0(16*i- 4)= dkinv(3,4)
    dlu0(16*i- 3)= dkinv(4,1)
    dlu0(16*i- 2)= dkinv(4,2)
    dlu0(16*i- 1)= dkinv(4,3)
    dlu0(16*i   )= dkinv(4,4)
 
    do i= 2, np
      iw1= 0
      iw2= 0
 
      do k= inumfi1l(i-1)+1, inumfi1l(i)
        iw1(fi1l(k))= k
      enddo
 
      do k= inumfi1u(i-1)+1, inumfi1u(i)
        iw2(fi1u(k))= k
      enddo
 
      do kk= inl(i-1)+1, inl(i)
        k= ial(kk)
 
        dkinv(1,1)= dlu0(16*k-15)
        dkinv(1,2)= dlu0(16*k-14)
        dkinv(1,3)= dlu0(16*k-13)
        dkinv(1,4)= dlu0(16*k-12)
        dkinv(2,1)= dlu0(16*k-11)
        dkinv(2,2)= dlu0(16*k-10)
        dkinv(2,3)= dlu0(16*k- 9)
        dkinv(2,4)= dlu0(16*k- 8)
        dkinv(3,1)= dlu0(16*k- 7)
        dkinv(3,2)= dlu0(16*k- 6)
        dkinv(3,3)= dlu0(16*k- 5)
        dkinv(3,4)= dlu0(16*k- 4)
        dkinv(4,1)= dlu0(16*k- 3)
        dkinv(4,2)= dlu0(16*k- 2)
        dkinv(4,3)= dlu0(16*k- 1)
        dkinv(4,4)= dlu0(16*k   )
 
        aik(1,1)= allu0(16*kk-15)
        aik(1,2)= allu0(16*kk-14)
        aik(1,3)= allu0(16*kk-13)
        aik(1,4)= allu0(16*kk-12)
        aik(2,1)= allu0(16*kk-11)
        aik(2,2)= allu0(16*kk-10)
        aik(2,3)= allu0(16*kk- 9)
        aik(2,4)= allu0(16*kk- 8)
        aik(3,1)= allu0(16*kk- 7)
        aik(3,2)= allu0(16*kk- 6)
        aik(3,3)= allu0(16*kk- 5)
        aik(3,4)= allu0(16*kk- 4)
        aik(4,1)= allu0(16*kk- 3)
        aik(4,2)= allu0(16*kk- 2)
        aik(4,3)= allu0(16*kk- 1)
        aik(4,4)= allu0(16*kk   )
 
        do jj= inu(k-1)+1, inu(k)
          j= iau(jj)
          if (iw1(j).eq.0.and.iw2(j).eq.0) cycle
 
          akj(1,1)= aulu0(16*jj-15)
          akj(1,2)= aulu0(16*jj-14)
          akj(1,3)= aulu0(16*jj-13)
          akj(1,4)= aulu0(16*jj-12)
          akj(2,1)= aulu0(16*jj-11)
          akj(2,2)= aulu0(16*jj-10)
          akj(2,3)= aulu0(16*jj- 9)
          akj(2,4)= aulu0(16*jj- 8)
          akj(3,1)= aulu0(16*jj- 7)
          akj(3,2)= aulu0(16*jj- 6)
          akj(3,3)= aulu0(16*jj- 5)
          akj(3,4)= aulu0(16*jj- 4)
          akj(4,1)= aulu0(16*jj- 3)
          akj(4,2)= aulu0(16*jj- 2)
          akj(4,3)= aulu0(16*jj- 1)
          akj(4,4)= aulu0(16*jj   )
 
          call ilu1b44 (rhs_aij, dkinv, aik, akj)
 
          if (j.eq.i) then
            dlu0(16*i-15)= dlu0(16*i-15) - rhs_aij(1,1)
            dlu0(16*i-14)= dlu0(16*i-14) - rhs_aij(1,2)
            dlu0(16*i-13)= dlu0(16*i-13) - rhs_aij(1,3)
            dlu0(16*i-12)= dlu0(16*i-12) - rhs_aij(1,4)
            dlu0(16*i-11)= dlu0(16*i-11) - rhs_aij(2,1)
            dlu0(16*i-10)= dlu0(16*i-10) - rhs_aij(2,2)
            dlu0(16*i- 9)= dlu0(16*i- 9) - rhs_aij(2,3)
            dlu0(16*i- 8)= dlu0(16*i- 8) - rhs_aij(2,4)
            dlu0(16*i- 7)= dlu0(16*i- 7) - rhs_aij(3,1)
            dlu0(16*i- 6)= dlu0(16*i- 6) - rhs_aij(3,2)
            dlu0(16*i- 5)= dlu0(16*i- 5) - rhs_aij(3,3)
            dlu0(16*i- 4)= dlu0(16*i- 4) - rhs_aij(3,4)
            dlu0(16*i- 3)= dlu0(16*i- 3) - rhs_aij(4,1)
            dlu0(16*i- 2)= dlu0(16*i- 2) - rhs_aij(4,2)
            dlu0(16*i- 1)= dlu0(16*i- 1) - rhs_aij(4,3)
            dlu0(16*i   )= dlu0(16*i   ) - rhs_aij(4,4)
          endif
 
          if (j.lt.i) then
            ij0= iw1(j)
            allu0(16*ij0-15)= allu0(16*ij0-15) - rhs_aij(1,1)
            allu0(16*ij0-14)= allu0(16*ij0-14) - rhs_aij(1,2)
            allu0(16*ij0-13)= allu0(16*ij0-13) - rhs_aij(1,3)
            allu0(16*ij0-12)= allu0(16*ij0-12) - rhs_aij(1,4)
            allu0(16*ij0-11)= allu0(16*ij0-11) - rhs_aij(2,1)
            allu0(16*ij0-10)= allu0(16*ij0-10) - rhs_aij(2,2)
            allu0(16*ij0- 9)= allu0(16*ij0- 9) - rhs_aij(2,3)
            allu0(16*ij0- 8)= allu0(16*ij0- 8) - rhs_aij(2,4)
            allu0(16*ij0- 7)= allu0(16*ij0- 7) - rhs_aij(3,1)
            allu0(16*ij0- 6)= allu0(16*ij0- 6) - rhs_aij(3,2)
            allu0(16*ij0- 5)= allu0(16*ij0- 5) - rhs_aij(3,3)
            allu0(16*ij0- 4)= allu0(16*ij0- 4) - rhs_aij(3,4)
            allu0(16*ij0- 3)= allu0(16*ij0- 3) - rhs_aij(4,1)
            allu0(16*ij0- 2)= allu0(16*ij0- 2) - rhs_aij(4,2)
            allu0(16*ij0- 1)= allu0(16*ij0- 1) - rhs_aij(4,3)
            allu0(16*ij0   )= allu0(16*ij0   ) - rhs_aij(4,4)
          endif
 
          if (j.gt.i) then
            ij0= iw2(j)
            aulu0(16*ij0-15)= aulu0(16*ij0-15) - rhs_aij(1,1)
            aulu0(16*ij0-14)= aulu0(16*ij0-14) - rhs_aij(1,2)
            aulu0(16*ij0-13)= aulu0(16*ij0-13) - rhs_aij(1,3)
            aulu0(16*ij0-12)= aulu0(16*ij0-12) - rhs_aij(1,4)
            aulu0(16*ij0-11)= aulu0(16*ij0-11) - rhs_aij(2,1)
            aulu0(16*ij0-10)= aulu0(16*ij0-10) - rhs_aij(2,2)
            aulu0(16*ij0- 9)= aulu0(16*ij0- 9) - rhs_aij(2,3)
            aulu0(16*ij0- 8)= aulu0(16*ij0- 8) - rhs_aij(2,4)
            aulu0(16*ij0- 7)= aulu0(16*ij0- 7) - rhs_aij(3,1)
            aulu0(16*ij0- 6)= aulu0(16*ij0- 6) - rhs_aij(3,2)
            aulu0(16*ij0- 5)= aulu0(16*ij0- 5) - rhs_aij(3,3)
            aulu0(16*ij0- 4)= aulu0(16*ij0- 4) - rhs_aij(3,4)
            aulu0(16*ij0- 3)= aulu0(16*ij0- 3) - rhs_aij(4,1)
            aulu0(16*ij0- 2)= aulu0(16*ij0- 2) - rhs_aij(4,2)
            aulu0(16*ij0- 1)= aulu0(16*ij0- 1) - rhs_aij(4,3)
            aulu0(16*ij0   )= aulu0(16*ij0   ) - rhs_aij(4,4)
          endif
 
        enddo
      enddo
 
      call ilu1a44 (dkinv, &
        dlu0(16*i-15), dlu0(16*i-14), dlu0(16*i-13), dlu0(16*i-12), &
        dlu0(16*i-11), dlu0(16*i-10), dlu0(16*i- 9), dlu0(16*i- 8), &
        dlu0(16*i- 7), dlu0(16*i- 6), dlu0(16*i- 5), dlu0(16*i- 4), &
        dlu0(16*i- 3), dlu0(16*i- 2), dlu0(16*i- 1), dlu0(16*i   ) )
      dlu0(16*i-15)= dkinv(1,1)
      dlu0(16*i-14)= dkinv(1,2)
      dlu0(16*i-13)= dkinv(1,3)
      dlu0(16*i-12)= dkinv(1,4)
      dlu0(16*i-11)= dkinv(2,1)
      dlu0(16*i-10)= dkinv(2,2)
      dlu0(16*i- 9)= dkinv(2,3)
      dlu0(16*i- 8)= dkinv(2,4)
      dlu0(16*i- 7)= dkinv(3,1)
      dlu0(16*i- 6)= dkinv(3,2)
      dlu0(16*i- 5)= dkinv(3,3)
      dlu0(16*i- 4)= dkinv(3,4)
      dlu0(16*i- 3)= dkinv(4,1)
      dlu0(16*i- 2)= dkinv(4,2)
      dlu0(16*i- 1)= dkinv(4,3)
      dlu0(16*i   )= dkinv(4,4)
    enddo
 
    deallocate (iw1, iw2)
  end subroutine form_ilu0_44
 
  !C
  !C***
  !C*** ILU1a44
  !C***
  !C
  !C    computes LU factorization of 4*4 Diagonal Block
  !C
  subroutine ilu1a44 (ALU, D11,D12,D13,D14,D21,D22,D23,D24,D31,D32,D33,D34,D41,D42,D43,D44)
    use hecmw_util
    implicit none
    real(kind=kreal) :: alu(4,4), pw(4)
    real(kind=kreal) :: d11,d12,d13,d14,d21,d22,d23,d24,d31,d32,d33,d34,d41,d42,d43,d44
    integer(kind=kint) :: i,j,k
 
    alu(1,1)= d11
    alu(1,2)= d12
    alu(1,3)= d13
    alu(1,4)= d14
    alu(2,1)= d21
    alu(2,2)= d22
    alu(2,3)= d23
    alu(2,4)= d24
    alu(3,1)= d31
    alu(3,2)= d32
    alu(3,3)= d33
    alu(3,4)= d34
    alu(4,1)= d41
    alu(4,2)= d42
    alu(4,3)= d43
    alu(4,4)= d44
 
    do k= 1, 4
      alu(k,k)= 1.d0/alu(k,k)
      do i= k+1, 4
        alu(i,k)= alu(i,k) * alu(k,k)
        do j= k+1, 4
          pw(j)= alu(i,j) - alu(i,k)*alu(k,j)
        enddo
        do j= k+1, 4
          alu(i,j)= pw(j)
        enddo
      enddo
    enddo
 
    return
  end subroutine ilu1a44
 
  !C
  !C***
  !C*** ILU1b44
  !C***
  !C
  !C    computes L_ik * D_k_INV * U_kj at ILU factorization
  !C    for 4*4 Block Type Matrix
  !C
  subroutine ilu1b44 (RHS_Aij, DkINV, Aik, Akj)
    use hecmw_util
    implicit none
    real(kind=kreal) :: rhs_aij(4,4), dkinv(4,4), aik(4,4), akj(4,4)
    real(kind=kreal) :: x1,x2,x3,x4
 
    !C
    !C-- 1st Col.
    x1= akj(1,1)
    x2= akj(2,1)
    x3= akj(3,1)
    x4= akj(4,1)
 
    x2= x2 - dkinv(2,1)*x1
    x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
    x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
 
    x4= dkinv(4,4)*  x4
    x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
    x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
    x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
 
    rhs_aij(1,1)=  aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
    rhs_aij(2,1)=  aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
    rhs_aij(3,1)=  aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
    rhs_aij(4,1)=  aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
 
    !C
    !C-- 2nd Col.
    x1= akj(1,2)
    x2= akj(2,2)
    x3= akj(3,2)
    x4= akj(4,2)
 
    x2= x2 - dkinv(2,1)*x1
    x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
    x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
 
    x4= dkinv(4,4)*  x4
    x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
    x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
    x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
 
    rhs_aij(1,2)=  aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
    rhs_aij(2,2)=  aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
    rhs_aij(3,2)=  aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
    rhs_aij(4,2)=  aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
 
    !C
    !C-- 3rd Col.
    x1= akj(1,3)
    x2= akj(2,3)
    x3= akj(3,3)
    x4= akj(4,3)
 
    x2= x2 - dkinv(2,1)*x1
    x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
    x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
 
    x4= dkinv(4,4)*  x4
    x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
    x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
    x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
 
    rhs_aij(1,3)=  aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
    rhs_aij(2,3)=  aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
    rhs_aij(3,3)=  aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
    rhs_aij(4,3)=  aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
 
    !C
    !C-- 4th Col.
    x1= akj(1,4)
    x2= akj(2,4)
    x3= akj(3,4)
    x4= akj(4,4)
 
    x2= x2 - dkinv(2,1)*x1
    x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
    x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
 
    x4= dkinv(4,4)*  x4
    x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
    x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
    x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
 
    rhs_aij(1,4)=  aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
    rhs_aij(2,4)=  aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
    rhs_aij(3,4)=  aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
    rhs_aij(4,4)=  aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
 
    return
  end subroutine ilu1b44
 
end module     hecmw_precond_bilu_44