sgesl Subroutine

  subroutine sgesl(a,lda,n,ipvt,b,job)

    !! Routine to solve the the real system  Ax = b  or  transp(A).x = b
    !! author: Cleve Moler, University of New Mexico, Argonne National Lab.
    !! author: P J Knight, CCFE, Culham Science Centre
    !! a(lda,n) : input real array : output from <A HREF="sgefa.html">sgefa</A>
    !! lda : input integer : leading dimension of the array A
    !! n : input integer : order of the matrix A
    !! ipvt(n) : input integer array : pivot vector from <CODE>sgefa</CODE>
    !! b(n) : input/output real array : RHS vector on input,
    !! solution vector x on output
    !! job : input integer : switch
    !! = 0         to solve  A*x = b ,
    !! = nonzero   to solve  transp(A)*x = b  where
    !! transp(A)  is the transpose
    !! This routine solves the real system  A*x = b  or  transp(A)*x = b
    !! using the factors computed by <A HREF="sgefa.html">sgefa</A>.
    !! <P>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 <CODE>lda</CODE>.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    integer, intent(in) :: lda,n,job
    integer, dimension(n), intent(in) :: ipvt
    real(dp), dimension(lda,n), intent(in) :: a
    real(dp), dimension(n), intent(inout) :: b

    !  Local variables

    integer :: k,kb,l,nm1
    real(dp) :: t

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    nm1 = n - 1

    if (job == 0) then  !  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
             end if
             call saxpy(n-k,t,a(k+1:n,k),1,b(k+1:n),1)
          end do
       end if

       !  Now solve  u*x = y

       do kb = 1, n
          k = n + 1 - kb
          b(k) = b(k)/a(k,k)
          t = -b(k)
          call saxpy(k-1,t,a(1:n,k),1,b(1:n),1)
       end do

    else  !  Solve  transp(A) * x = b

       !  First solve  transp(u)*y = b

       do k = 1, n
          t = sdot(k-1,a(:,k),1,b(:),1)
          b(k) = (b(k) - t)/a(k,k)
       end do

       !  Now solve transp(l)*x = y

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

    end if

  end subroutine sgesl