scale_matrix.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_matrix
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: matrix_solver_tridiagonal

  interface matrix_solver_tridiagonal
     module procedure matrix_solver_tridiagonal_1d
     module procedure matrix_solver_tridiagonal_3d
  end interface matrix_solver_tridiagonal

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine matrix_solver_tridiagonal_1d( &
       KA, &
       ud, &
       md, &
       ld, &
       iv, &
       ov  )
    implicit none

    integer,  intent(in)  :: KA     ! array size
    real(RP), intent(in)  :: ud(KA) ! upper  diagonal
    real(RP), intent(in)  :: md(KA) ! middle diagonal
    real(RP), intent(in)  :: ld(KA) ! lower  diagonal
    real(RP), intent(in)  :: iv(KA) ! input  vector
    real(RP), intent(out) :: ov(KA) ! output vector

    real(RP) :: c(KA)
    real(RP) :: d(KA)

    integer :: k
    !---------------------------------------------------------------------------

    ! foward reduction
    c(1) = ud(1) / md(1)
    d(1) = iv(1) / md(1)
    do k = 2, ka
       c(k) =           ud(k)            / ( md(k) - ld(k) * c(k-1) )
       d(k) = ( iv(k) - ld(k) * d(k-1) ) / ( md(k) - ld(k) * c(k-1) )
    enddo

    ! backward substitution
    ov(ka) = d(ka)
    do k = ka-1, 1, -1
       ov(k) = d(k) - c(k) * ov(k+1)
    enddo

    return
  end subroutine matrix_solver_tridiagonal_1d

  !-----------------------------------------------------------------------------
  subroutine matrix_solver_tridiagonal_3d( &
       KA,         &
       IA, IS, IE, &
       JA, JS, JE, &
       ud,         &
       md,         &
       ld,         &
       iv,         &
       ov          )
    implicit none

    integer,  intent(in)  :: KA           ! array size
    integer,  intent(in)  :: IA, IS, IE   ! array size
    integer,  intent(in)  :: JA, JS, JE   ! array size
    real(RP), intent(in)  :: ud(KA,IA,JA) ! upper  diagonal
    real(RP), intent(in)  :: md(KA,IA,JA) ! middle diagonal
    real(RP), intent(in)  :: ld(KA,IA,JA) ! lower  diagonal
    real(RP), intent(in)  :: iv(KA,IA,JA) ! input  vector
    real(RP), intent(out) :: ov(KA,IA,JA) ! output vector

    real(RP) :: c(KA,IA,JA)
    real(RP) :: d(KA,IA,JA)

    integer :: k, i, j
    !---------------------------------------------------------------------------

    ! foward reduction
    do j = js, je
    do i = is, ie
       c(1,i,j) = ud(1,i,j) / md(1,i,j)
       d(1,i,j) = iv(1,i,j) / md(1,i,j)
    enddo
    enddo

    !$omp parallel do default(none) &
    !$omp shared(JS,JE,IS,IE,KA,ud,md,ld,iv,c,d) &
    !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
    do j = js, je
    do i = is, ie
    do k = 2,  ka
       c(k,i,j) =               ud(k,i,j)                &
                / ( md(k,i,j) - ld(k,i,j) * c(k-1,i,j) )
       d(k,i,j) = ( iv(k,i,j) - ld(k,i,j) * d(k-1,i,j) ) &
                / ( md(k,i,j) - ld(k,i,j) * c(k-1,i,j) )
    enddo
    enddo
    enddo

    ! backward substitution
    do j = js, je
    do i = is, ie
       ov(ka,i,j) = d(ka,i,j)
    enddo
    enddo

    do j = js, je
    do i = is, ie
    do k = ka-1, 1, -1
       ov(k,i,j) = d(k,i,j) - c(k,i,j) * ov(k+1,i,j)
    enddo
    enddo
    enddo

    return
  end subroutine matrix_solver_tridiagonal_3d

end module scale_matrix