scale_matrix.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_matrix
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  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
  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine matrix_solver_tridiagonal_1d( &
       KA, KS, KE, &
       ud, md, ld, &
       iv,         &
       ov          )
    implicit none
    integer,  intent(in)  :: ka, ks, ke

    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(ks) = ud(ks) / md(ks)
    d(ks) = iv(ks) / md(ks)
    do k = ks+1, ke
       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(ke) = d(ke)
    do k = ke-1, ks, -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