dbndsl Subroutine

    subroutine dbndsl(Mode, g, Mdg, Nb, Ip, Ir, x, n, Rnorm)

      integer,intent(in) :: Mode !! Set by the user to one of the values 1, 2, or
                                 !! 3. These values respectively indicate that
                                 !! the solution of `AX = B`, `YR = H` or `RZ = W` is
                                 !! required.
      integer,intent(in) :: Mdg !! The number of rows in the working array
                                !! `G(*,*)`.  The value of `MDG` should be `>= MU`.
                                !! The value of `MU` is defined in the abstract
                                !! of these subprograms.
                                !!
                                !! This argument has the same meaning and
                                !! contents as following the last call to [[DBNDAC]].
      real(wp),intent(in) :: g(Mdg, *) !! `G(MDG,NB+1)`
                                       !!
                                       !! This argument has the same meaning and
                                       !! contents as following the last call to [[DBNDAC]].
      integer,intent(in) :: Nb !! This argument has the same meaning and
                               !! contents as following the last call to [[DBNDAC]].
      integer,intent(in) :: Ip !! This argument has the same meaning and
                               !! contents as following the last call to [[DBNDAC]].
      integer,intent(in) :: Ir !! This argument has the same meaning and
                               !! contents as following the last call to [[DBNDAC]].
      real(wp),intent(inout) :: x(*) !! `X(N)`
                                     !!
                                     !! *Input* With mode=2 or 3 this array contains,
                                     !! respectively, the right-side vectors H or W of
                                     !! the systems YR = H or RZ = W.
                                     !!
                                     !! *Output* This array contains the solution vectors `X`,
                                     !! `Y` or `Z` of the systems `AX = B`, `YR = H` or
                                     !! `RZ = W` depending on the value of `MODE`=1,
                                     !! 2 or 3.
      integer,intent(in) :: n !! The number of variables in the solution
                              !! vector.  If any of the `N` diagonal terms are
                              !! zero the subroutine [[DBNDSL]] prints an
                              !! appropriate message.  This condition is
                              !! considered an error.
      real(wp),intent(out) :: Rnorm !! If `MODE=1`, `RNORM` is the Euclidean length of the
                                    !! residual vector `AX-B`.  When `MODE=2` or `3` RNORM`
                                    !! is set to zero.

      real(wp) :: rsq, s
      integer :: i, i1, i2, ie, ii, iopt, irm1, ix, j, &
                 jg, l, nerr, np1

      real(wp), parameter :: zero = 0.0_wp

      main: block

         Rnorm = zero
         select case (Mode)
         case (1)
            ! ALG. STEP 26
            do j = 1, n
               x(j) = g(j, Nb + 1)
            end do
            rsq = zero
            np1 = n + 1
            irm1 = Ir - 1
            if (np1 <= irm1) then
               do j = np1, irm1
                  rsq = rsq + g(j, Nb + 1)**2
               end do
               Rnorm = sqrt(rsq)
            end if
         case (2)
            do j = 1, n
               s = zero
               if (j /= 1) then
                  i1 = max(1, j - Nb + 1)
                  i2 = j - 1
                  do i = i1, i2
                     l = j - i + 1 + max(0, i - Ip)
                     s = s + x(i)*g(i, l)
                  end do
               end if
               l = max(0, j - Ip)
               if (g(j, l + 1) == 0) exit main
               x(j) = (x(j) - s)/g(j, l + 1)
            end do
            return
         end select

         ! MODE = 3
         do ii = 1, n
            i = n + 1 - ii
            s = zero
            l = max(0, i - Ip)
            if (i /= n) then
               ie = min(n + 1 - i, Nb)
               do j = 2, ie
                  jg = j + l
                  ix = i - 1 + j
                  s = s + g(i, jg)*x(ix)
               end do
            end if
            if (g(i, l + 1) == 0) exit main
            x(i) = (x(i) - s)/g(i, l + 1)
         end do

         return

      end block main

      ! error handling
      nerr = 1
      iopt = 2
      write (*, *) 'A zero diagonal term is in the n by n upper triangular matrix.'

   end subroutine dbndsl