dsolve Subroutine

public pure subroutine dsolve(n, t, ldt, b, job)

Uses

  • proc~~dsolve~~UsesGraph proc~dsolve dsolve module~blas_interfaces blas_interfaces proc~dsolve->module~blas_interfaces module~odrpack_kinds odrpack_kinds proc~dsolve->module~odrpack_kinds module~blas_interfaces->module~odrpack_kinds iso_fortran_env iso_fortran_env module~odrpack_kinds->iso_fortran_env

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.

References: * Dongarra J.J., Bunch J.R., Moler C.B., Stewart G.W., LINPACK Users Guide, SIAM, 1979.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

The number of rows and columns of data in array t.
real(kind=wp), intent(in) :: t(ldt,n)

The upper or lower tridiagonal system.
integer, intent(in) :: ldt

The leading dimension of array t.
real(kind=wp), intent(inout) :: b(n)

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.

Calls

proc~~dsolve~~CallsGraph proc~dsolve dsolve interface~daxpy daxpy proc~dsolve->interface~daxpy interface~ddot ddot proc~dsolve->interface~ddot

Called by

proc~~dsolve~~CalledByGraph proc~dsolve dsolve proc~dodstp dodstp proc~dodstp->proc~dsolve proc~dvevtr dvevtr proc~dodstp->proc~dvevtr proc~dvevtr->proc~dsolve proc~dodlm dodlm proc~dodlm->proc~dodstp proc~dodvcv dodvcv proc~dodvcv->proc~dodstp proc~dodmn dodmn proc~dodmn->proc~dodlm proc~dodmn->proc~dodvcv proc~doddrv doddrv proc~doddrv->proc~dodmn proc~dodcnt dodcnt proc~dodcnt->proc~doddrv proc~odr odr proc~odr->proc~dodcnt proc~odr_long_c odr_long_c proc~odr_long_c->proc~odr proc~odr_medium_c odr_medium_c proc~odr_medium_c->proc~odr proc~odr_short_c odr_short_c proc~odr_short_c->proc~odr program~example1 example1 program~example1->proc~odr program~example2 example2 program~example2->proc~odr program~example3 example3 program~example3->proc~odr program~example4 example4 program~example4->proc~odr program~example5 example5 program~example5->proc~odr

Variables

Type Visibility Attributes Name Initial
real(kind=wp), public :: temp
integer, public :: j1
integer, public :: j
integer, public :: jn

Source Code

   pure subroutine dsolve(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.
   !!
   !! References:
   !! * Dongarra J.J., Bunch J.R., Moler C.B., Stewart G.W., *LINPACK Users Guide*, SIAM, 1979.
      ! @note: This is one of the most time-consuming subroutines in ODRPACK (~25% of total).
      ! Routines Called  DAXPY, DDOT
      ! Date Written   920220   (YYMMDD)
      ! Revision Date  920619   (YYMMDD)

      use odrpack_kinds, only: zero
      use blas_interfaces, only: daxpy, ddot

      integer, intent(in) :: n
         !! The number of rows and columns of data in array `t`.
      real(wp), intent(in) :: t(ldt, n)
         !! The upper or lower tridiagonal system.
      integer, intent(in) :: ldt
         !! The leading dimension of array `t`.
      real(wp), intent(inout) :: b(n)
         !! 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

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

      elseif (job == 2) then
         ! 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)
            call daxpy(j, temp, t(1, j + 1), 1, b(1), 1)
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do

      elseif (job == 3) then
         ! 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) - ddot(jn - j + 1, t(j + 1, j), 1, b(j + 1), 1)
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do

      elseif (job == 4) then
         ! 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) - ddot(j - 1, t(1, j), 1, b(1), 1)
            if (t(j, j) /= zero) then
               b(j) = b(j)/t(j, j)
            else
               b(j) = zero
            end if
         end do
      end if

   end subroutine dsolve