dgesl Subroutine

public subroutine dgesl(a, lda, n, ipvt, b, job)

solve the real system a*x=b or trans(a)*x=b using the factors computed by dgeco or dgefa.

to compute inverse(a) * c where c is a matrix with p columns:

  call dgeco(a,lda,n,ipvt,rcond,z)
  if (rcond is too small) go to ...
  do j = 1, p
     call dgesl(a,lda,n,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 dgeco has set rcond > 0.0 or dgefa has set info == 0 .

Author

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

Revision history

  • 780814 date written
  • 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) :: a(lda,*)

(lda, n) the output from dgeco or dgefa.
integer, intent(in) :: lda

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

the order of the matrix a.
integer, intent(in) :: ipvt(*)

the pivot vector from dgeco or dgefa.
real(kind=wp), intent(out) :: b(*)

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~~dgesl~~CallsGraph proc~dgesl dvode_linpack_module::dgesl proc~daxpy dvode_blas_module::daxpy proc~dgesl->proc~daxpy proc~ddot dvode_blas_module::ddot proc~dgesl->proc~ddot

Called by

proc~~dgesl~~CalledByGraph proc~dgesl dvode_linpack_module::dgesl proc~dvsol dvode_module::dvode_t%dvsol proc~dvsol->proc~dgesl 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 dgesl(a,lda,n,ipvt,b,job)

       integer,intent(in) :: lda !! the leading dimension of the array `a`.
       integer,intent(in) :: n !! the order of the matrix `a`.
       integer,intent(in) :: ipvt(*) !! the pivot vector from [[dgeco]] or [[dgefa]].
       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) :: a(lda,*) !! `(lda, n)` the output from [[dgeco]] or [[dgefa]].
       real(wp),intent(out) :: b(*) !! the solution vector `x`.

       real(wp) :: t
       integer :: k , kb , l , nm1
       
#ifdef HAS_BLAS
       ! user is linking against an external BLAS library
       double precision,external :: ddot 
#endif

       nm1 = n - 1
       if ( job/=0 ) then

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

          do k = 1 , n
             t = ddot(k-1,a(1,k),1,b(1),1)
             b(k) = (b(k)-t)/a(k,k)
          enddo

         ! now solve trans(l)*x = y

          if ( nm1>=1 ) then
             do kb = 1 , nm1
                k = n - kb
                b(k) = b(k) + ddot(n-k,a(k+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
       else

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

          if ( nm1>=1 ) then
             do k = 1 , nm1
                l = ipvt(k)
                t = b(l)
                if ( l/=k ) then
                   b(l) = b(k)
                   b(k) = t
                endif
                call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
             enddo
          endif

         ! now solve  u*x = y

          do kb = 1 , n
             k = n + 1 - kb
             b(k) = b(k)/a(k,k)
             t = -b(k)
             call daxpy(k-1,t,a(1,k),1,b(1),1)
          enddo
       endif

   end subroutine dgesl