sgtsl.f90

Go to the documentation of this file.

!** from netlib, tue aug 28 08:28:34 edt 1990 ***
!
subroutine sgtsl(n, c, d, e, b, info)

  use itm_types

  implicit none

  integer(itm_i4), intent(in) :: n
  real(r8), dimension(n + 2), intent(inout) :: c, d, e, b
  integer(itm_i4), intent(out) :: info

!-----------------------------------------------------------------------------
!     sgtsl given a general tridiagonal matrix and a right hand
!     side will find the solution.
!
!     on entry
!
!        n       integer
!                is the order of the tridiagonal matrix.
!
!        c       real(n)
!                is the subdiagonal of the tridiagonal matrix.
!                c(2) through c(n) should contain the subdiagonal.
!                on output c is destroyed.
!
!        d       real(n)
!                is the diagonal of the tridiagonal matrix.
!                on output d is destroyed.
!
!        e       real(n)
!                is the superdiagonal of the tridiagonal matrix.
!                e(1) through e(n-1) should contain the superdiagonal.
!                on output e is destroyed.
!
!        b       real(n)
!                is the right hand side vector.
!
!     on return
!
!        b       is the solution vector.
!
!        info    integer
!                = 0 normal value.
!                = k if the k-th element of the diagonal becomes
!                    exactly zero.  the subroutine returns when
!                    this is detected.
!
!     linpack. this version dated 08/14/78.
!     jack dongarra, argonne national laboratory.
!
!     no externals
!     fortran abs
!
!     internal variables
!----------------------------------------------------------------------------

  integer(itm_i4) :: k, kb, kp1, nm1, nm2
  real(r8) :: t

!-- begin block permitting ...exits to 100
  info = 0
  c(1) = d(1)
  nm1 = n - 1

  if (nm1 >= 1) then
    d(1) = e(1)
    e(1) = 0.0_r8
    e(n) = 0.0_r8
    do k = 1, nm1
      kp1 = k + 1
!-- find the largest of the two rows
      if (abs(c(kp1)) >= abs(c(k))) then
!-- interchange row
        t = c(kp1)
        c(kp1) = c(k)
        c(k) = t
        t = d(kp1)
        d(kp1) = d(k)
        d(k) = t
        t = e(kp1)
        e(kp1) = e(k)
        e(k) = t
        t = b(kp1)
        b(kp1) = b(k)
        b(k) = t
      end if
!-- zero elements
      if (c(k) == 0.0_r8) then
        info = k
!-- exit
        return
      end if
      t = -c(kp1) / c(k)
      c(kp1) = d(kp1) + t * d(k)
      d(kp1) = e(kp1) + t * e(k)
      e(kp1) = 0.0_r8
      b(kp1) = b(kp1) + t * b(k)
    end do

  end if

  if (c(n) == 0._r8) then
    info = n
    return
  end if

!-- back solve

  nm2 = n - 2
  b(n) = b(n) / c(n)
  if (n == 1) return
  b(nm1) = (b(nm1) - d(nm1) * b(n)) / c(nm1)
  if (nm2 < 1) return
  do kb = 1, nm2
    k = nm2 - kb + 1
    b(k) = (b(k) - d(k) * b(k + 1) - e(k) * b(k + 2)) / c(k)
  end do

  return
end subroutine sgtsl