solve_trl Subroutine

   pure subroutine solve_trl(n, t, ldt, b, job)
   !! Solve systems of the form:
   !!
   !!  `t * x = b  or  trans(t) * x = b`
   !!
   !! where `t` is an upper or lower triangular matrix of order `n`, and the solution `x`
   !! overwrites the RHS `b`.
      ! Adapted from LINPACK subroutine DTRSL.
      ! @note: This is one of the most time-consuming subroutines in ODRPACK (~25% of total).

      use odrpack_kinds, only: zero

      integer, intent(in) :: n
         !! Number of rows and columns of data in array `t`.
      real(wp), intent(in) :: t(ldt, n)
         !! Upper or lower tridiagonal system.
      integer, intent(in) :: ldt
         !! Leading dimension of array `t`.
      real(wp), intent(inout) :: b(:)
         !! On input: the right hand side; On exit: the solution.
      integer, intent(in) :: job
         !! What kind of system is to be solved:
         !!   1: Solve `t * x = b`, where `t` is lower triangular,
         !!   2: Solve `t * x = b`, where `t` is upper triangular,
         !!   3: Solve `trans(t) * x = b`, where `t` is lower triangular,
         !!   4: Solve `trans(t) * x = b`, where `t` is upper triangular.

      ! Local scalars
      real(wp) :: temp
      integer :: j1, j, jn

      ! Variable Definitions (alphabetically)
      !  B:       On input: the right hand side; On exit: the solution
      !  J1:      The first nonzero entry in T.
      !  J:       An indexing variable.
      !  JN:      The last nonzero entry in T.
      !  JOB:     What kind of system is to be solved, where if JOB is
      !           1   Solve T*X=B, T lower triangular,
      !           2   Solve T*X=B, T upper triangular,
      !           3   Solve trans(T)*X=B, T lower triangular,
      !           4   Solve trans(T)*X=B, T upper triangular.
      !  LDT:     The leading dimension of array T.
      !  N:       The number of rows and columns of data in array T.
      !  T:       The upper or lower tridiagonal system.

      !  Find first nonzero diagonal entry in T
      j1 = 0
      do j = 1, n
         if (j1 == 0 .and. t(j, j) /= zero) then
            j1 = j
         elseif (t(j, j) == zero) then
            b(j) = zero
         end if
      end do
      if (j1 == 0) return

      ! Find last nonzero diagonal entry in T
      jn = 0
      do j = n, j1, -1
         if (jn == 0 .and. t(j, j) /= zero) then
            jn = j
         elseif (t(j, j) == zero) then
            b(j) = zero
         end if
      end do

      select case (job)
      case (1)
         ! Solve T*X=B for T lower triangular
         b(j1) = b(j1)/t(j1, j1)
         do j = j1 + 1, jn
            temp = -b(j - 1)
            b(j:jn) = b(j:jn) + temp*t(j:jn, j - 1)
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do

      case (2)
         ! Solve T*X=B for T upper triangular.
         b(jn) = b(jn)/t(jn, jn)
         do j = jn - 1, j1, -1
            temp = -b(j + 1)
            b(1:j) = b(1:j) + temp*t(1:j, j + 1)
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do

      case (3)
         ! Solve trans(T)*X=B for T lower triangular.
         b(jn) = b(jn)/t(jn, jn)
         do j = jn - 1, j1, -1
            b(j) = b(j) - dot_product(t(j + 1:jn + 1, j), b(j + 1:jn + 1))
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do

      case (4)
         ! Solve trans(T)*X=B for T upper triangular
         b(j1) = b(j1)/t(j1, j1)
         do j = j1 + 1, jn
            b(j) = b(j) - dot_product(t(1:j - 1, j), b(1:j - 1))
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do
      case default
         error stop "Invalid value of JOB."
      end select

   end subroutine solve_trl