ImsLinearMisc.f90

Go to the documentation of this file.

MODULE imslinearmisc
 
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dzero, done, dem30
 
  private
  public :: ims_misc_thomas
  public :: ims_misc_dvscale
 
CONTAINS
 
  !> @brief Tridiagonal solve using the Thomas algorithm
    !!
    !! Subroutine to solve tridiagonal linear equations using the
    !! Thomas algorithm.
    !!
  !<
  subroutine ims_misc_thomas(n, tl, td, tu, b, x, w)
    implicit none
    ! -- dummy variables
    integer(I4B), intent(in) :: n !< number of matrix rows
    real(dp), dimension(n), intent(in) :: tl !< lower matrix terms
    real(dp), dimension(n), intent(in) :: td !< diagonal matrix terms
    real(dp), dimension(n), intent(in) :: tu !< upper matrix terms
    real(dp), dimension(n), intent(in) :: b !< right-hand side vector
    real(dp), dimension(n), intent(inout) :: x !< solution vector
    real(dp), dimension(n), intent(inout) :: w !< work vector
    ! -- local variables
    integer(I4B) :: j
    real(dp) :: bet
    real(dp) :: beti
    !
    ! -- initialize variables
    w(1) = dzero
    bet = td(1)
    beti = done / bet
    x(1) = b(1) * beti
    !
    ! -- decomposition and forward substitution
    do j = 2, n
      w(j) = tu(j - 1) * beti
      bet = td(j) - tl(j) * w(j)
      beti = done / bet
      x(j) = (b(j) - tl(j) * x(j - 1)) * beti
    end do
    !
    ! -- backsubstitution
    do j = n - 1, 1, -1
      x(j) = x(j) - w(j + 1) * x(j + 1)
    end do
  end subroutine ims_misc_thomas
 
  !
  !> @ brief Scale X and RHS
  !!
  !!  Scale X and B to avoid big or small values. Scaling value is the
  !!  maximum ABS(X).
  !<
  SUBROUTINE ims_misc_dvscale(IOPT, NEQ, DSCALE, X, B)
    ! -- dummy variables
    integer(I4B), INTENT(IN) :: iopt !< flag to scale (0) or unscale the system of equations
    integer(I4B), INTENT(IN) :: neq !< number of equations
    real(dp), INTENT(INOUT) :: dscale !< scaling value
    real(dp), DIMENSION(NEQ), INTENT(INOUT) :: x !< dependent variable
    real(dp), DIMENSION(NEQ), INTENT(INOUT) :: b !< right-hand side
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- SCALE SCALE AMAT, X, AND B
    IF (iopt == 0) THEN
      dscale = abs(dscale)
      if (dscale < dem30) then
        dscale = done
      end if
      DO n = 1, neq
        b(n) = b(n) / dscale
        x(n) = x(n) / dscale
      END DO
      !
      ! -- UNSCALE X, AND B
    ELSE
      DO n = 1, neq
        b(n) = b(n) * dscale
        x(n) = x(n) * dscale
      END DO
    END IF
  end SUBROUTINE ims_misc_dvscale
 
END MODULE imslinearmisc