dgbsl Subroutine

public pure subroutine dgbsl(abd, lda, n, ml, mu, ipvt, b, job)

Uses

  • proc~~dgbsl~~UsesGraph proc~dgbsl dgbsl module~blas_interfaces blas_interfaces proc~dgbsl->module~blas_interfaces module~daskr_kinds daskr_kinds module~blas_interfaces->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

DGBSL solves a real banded system factored by DGBCO or DGBFA.

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: abd(lda,n)
integer, intent(in) :: lda
integer, intent(in) :: n
integer, intent(in) :: ml
integer, intent(in) :: mu
integer, intent(in) :: ipvt(n)
real(kind=rk), intent(inout) :: b(n)
integer, intent(in) :: job

Calls

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

Called by

proc~~dgbsl~~CalledByGraph proc~dgbsl dgbsl proc~dslvd dslvd proc~dslvd->proc~dgbsl proc~psol_banpre psol_banpre proc~psol_banpre->proc~dgbsl

Variables

Type Visibility Attributes Name Initial
integer, public :: k
integer, public :: l
integer, public :: la
integer, public :: lb
integer, public :: lm
integer, public :: m
real(kind=rk), public :: t

Source Code

   pure subroutine dgbsl(abd, lda, n, ml, mu, ipvt, b, job)
   !! DGBSL solves a real banded system factored by DGBCO or DGBFA.
   !
   !  Discussion:
   !
   !    DGBSL can solve either A * X = B  or  A' * X = B.
   !
   !    A division by zero will occur if the input factor contains a
   !    zero on the diagonal.  Technically this indicates singularity
   !    but it is often caused by improper arguments or improper
   !    setting of LDA.  It will not occur if the subroutines are
   !    called correctly and if DGBCO has set 0.0 < RCOND
   !    or DGBFA has set INFO == 0.
   !
   !    To compute inverse(A) * C  where C is a matrix with P columns:
   !
   !      call dgbco ( abd, lda, n, ml, mu, ipvt, rcond, z )
   !
   !      if ( rcond is too small ) then
   !        exit
   !      end if
   !
   !      do j = 1, p
   !        call dgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 )
   !      end do
   !
   !  Licensing:
   !
   !    This code is distributed under the GNU LGPL license.
   !
   !  Modified:
   !
   !    17 May 2005
   !
   !  Author:
   !
   !    FORTRAN90 version by John Burkardt.
   !
   !  Reference:
   !
   !    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
   !    LINPACK User's Guide,
   !    SIAM, 1979,
   !    ISBN13: 978-0-898711-72-1,
   !    LC: QA214.L56.
   !
   !  Parameters:
   !
   !    Input, real ( kind = 8 ) ABD(LDA,N), the output from DGBCO or DGBFA.
   !
   !    Input, integer ( kind = 4 ) LDA, the leading dimension of the array ABD.
   !
   !    Input, integer ( kind = 4 ) N, the order of the matrix.
   !
   !    Input, integer ( kind = 4 ) ML, MU, the number of diagonals below and above
   !    the main diagonal.  0 <= ML < N, 0 <= MU < N.
   !
   !    Input, integer ( kind = 4 ) IPVT(N), the pivot vector from DGBCO or DGBFA.
   !
   !    Input/output, real ( kind = 8 ) B(N).  On input, the right hand side.
   !    On output, the solution.
   !
   !    Input, integer ( kind = 4 ) JOB, job choice.
   !    0, solve A*X=B.
   !    nonzero, solve A'*X=B.

      use blas_interfaces, only: daxpy, ddot

      real(rk), intent(in) :: abd(lda, n)
      integer, intent(in) :: lda
      integer, intent(in) :: n
      integer, intent(in) :: ml
      integer, intent(in) :: mu
      integer, intent(in) :: ipvt(n)
      real(rk), intent(inout) :: b(n)
      integer, intent(in) :: job

      integer :: k
      integer :: l
      integer :: la
      integer :: lb
      integer :: lm
      integer :: m
      real(rk) :: t

      m = mu + ml + 1
   !
   !  JOB = 0, Solve A * x = b.
   !
   !  First solve L * y = b.
   !
      if (job == 0) then

         if (0 < ml) then

            do k = 1, n - 1
               lm = min(ml, n - k)
               l = ipvt(k)
               t = b(l)
               if (l /= k) then
                  b(l) = b(k)
                  b(k) = t
               end if
               call daxpy(lm, t, abd(m + 1, k), 1, b(k + 1), 1)
            end do

         end if
   !
   !  Now solve U * x = y.
   !
         do k = n, 1, -1
            b(k) = b(k)/abd(m, k)
            lm = min(k, m) - 1
            la = m - lm
            lb = k - lm
            t = -b(k)
            call daxpy(lm, t, abd(la, k), 1, b(lb), 1)
         end do
   !
   !  JOB nonzero, solve A' * x = b.
   !
   !  First solve U' * y = b.
   !
      else

         do k = 1, n
            lm = min(k, m) - 1
            la = m - lm
            lb = k - lm
            t = ddot(lm, abd(la, k), 1, b(lb), 1)
            b(k) = (b(k) - t)/abd(m, k)
         end do
   !
   !  Now solve L' * x = y.
   !
         if (0 < ml) then

            do k = n - 1, 1, -1
               lm = min(ml, n - k)
               b(k) = b(k) + ddot(lm, abd(m + 1, k), 1, b(k + 1), 1)
               l = ipvt(k)
               if (l /= k) then
                  t = b(l)
                  b(l) = b(k)
                  b(k) = t
               end if
            end do

         end if

      end if

   end subroutine dgbsl