dgbsl Subroutine

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

solve the real band system a*x=b or trans(a)*x=b using the factors computed by dgbco or dgbfa.

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) go to ...
  do j = 1, p
     call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
  end do

Error condition

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 rcond > 0.0 or dgbfa has set info == 0.

Author

  • moler, c. b., (u. of new mexico)

Revision history

  • 780814 date written
  • 890531 changed all specific intrinsics to generic. (wrb)
  • 890831 modified array declarations. (wrb)
  • 890831 revision date from version 3.2
  • 891214 prologue converted to version 4.0 format. (bab)
  • 900326 removed duplicate information from description section. (wrb)
  • 920501 reformatted the references section. (wrb)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: abd(lda,*)

the output from dgbco or dgbfa.
integer, intent(in) :: lda

the leading dimension of the array abd.
integer, intent(in) :: n

the order of the original matrix.
integer, intent(in) :: ml

number of diagonals below the main diagonal.
integer, intent(in) :: mu

number of diagonals above the main diagonal.
integer, intent(in) :: ipvt(*)

the pivot vector from dgbco or dgbfa.
real(kind=wp), intent(inout) :: b(*)
  • input: the right hand side vector.
  • output: the solution vector x.
integer, intent(in) :: job
  • 0 to solve a*x = b ,
  • nonzero to solve trans(a)*x = b , where trans(a) is the transpose.

Calls

proc~~dgbsl~~CallsGraph proc~dgbsl dvode_linpack_module::dgbsl proc~daxpy dvode_blas_module::daxpy proc~dgbsl->proc~daxpy proc~ddot dvode_blas_module::ddot proc~dgbsl->proc~ddot

Called by

proc~~dgbsl~~CalledByGraph proc~dgbsl dvode_linpack_module::dgbsl proc~dvsol dvode_module::dvode_t%dvsol proc~dvsol->proc~dgbsl proc~dvnlsd dvode_module::dvode_t%dvnlsd proc~dvnlsd->proc~dvsol proc~dvstep dvode_module::dvode_t%dvstep proc~dvstep->proc~dvnlsd proc~dvode dvode_module::dvode_t%dvode proc~dvode->proc~dvstep

Source Code

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

       integer,intent(in) :: lda !! the leading dimension of the array  `abd`.
       integer,intent(in) :: n !! the order of the original matrix.
       integer,intent(in) :: ml !! number of diagonals below the main diagonal.
       integer,intent(in) :: mu !! number of diagonals above the main diagonal.
       integer,intent(in) :: ipvt(*) !! the pivot vector from [[dgbco]] or [[dgbfa]].
       integer,intent(in) :: job !! * 0         to solve  `a*x = b` ,
                                 !! * nonzero   to solve  `trans(a)*x = b` , where
                                 !!             `trans(a)`  is the transpose.
       real(wp),intent(in) :: abd(lda,*) !! the output from [[dgbco]] or [[dgbfa]].
       real(wp),intent(inout) :: b(*) !! * **input:** the right hand side vector.
                                      !! * **output:** the solution vector `x`.

       real(wp) :: t
       integer :: k , kb , l , la , lb , lm , m , nm1

#ifdef HAS_BLAS
       ! user is linking against an external BLAS library
       double precision,external :: ddot 
#endif

       m = mu + ml + 1
       nm1 = n - 1
       if ( job/=0 ) then

         ! job = nonzero, solve  trans(a) * x = b
         ! first solve  trans(u)*y = b

          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)
          enddo

         ! now solve trans(l)*x = y

          if ( ml/=0 ) then
             if ( nm1>=1 ) then
                do kb = 1 , nm1
                   k = n - kb
                   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
                   endif
                enddo
             endif
          endif
       else

         ! job = 0 , solve  a * x = b
         ! first solve l*y = b

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

         ! now solve  u*x = y

          do kb = 1 , n
             k = n + 1 - kb
             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)
          enddo
       endif

   end subroutine dgbsl