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-1
       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
    d(ke) = ( iv(ke) - ld(ke) * d(ke-1) ) / ( md(ke) - ld(ke) * c(ke-1) )
 
    ! 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, KS, KE, &
       IA, IS, IE, &
       JA, JS, JE, &
       ud,         &
       md,         &
       ld,         &
       iv,         &
       ov,         &
       mask        )
    implicit none
 
    integer,  intent(in)  :: KA, KS, KE   ! 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
 
    logical,  intent(in), optional :: mask(IA,JA)
 
    integer :: i, j
    !---------------------------------------------------------------------------
 
    if ( present(mask) ) then
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JS,JE,IS,IE,KA,KS,KE,ud,md,ld,iv,ov,mask) &
       !$omp private(i,j)
       do j = js, je
       do i = is, ie
          if ( mask(i,j) ) then
             call matrix_solver_tridiagonal_1d( ka, ks, ke, &
                                                ud(:,i,j), md(:,i,j), ld(:,i,j), & ! (in)
                                                iv(:,i,j),                       & ! (in)
                                                ov(:,i,j)                        ) ! (out)
          end if
       enddo
       enddo
    else
       !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
       !$omp shared(JS,JE,IS,IE,KA,KS,KE,ud,md,ld,iv,ov) &
       !$omp private(i,j)
       do j = js, je
       do i = is, ie
          call matrix_solver_tridiagonal_1d( ka, ks, ke, &
                                             ud(:,i,j), md(:,i,j), ld(:,i,j), & ! (in)
                                             iv(:,i,j),                       & ! (in)
                                             ov(:,i,j)                        ) ! (out)
       enddo
       enddo
    end if
 
    return
  end subroutine matrix_solver_tridiagonal_3d
 
end module scale_matrix