lubksb Subroutine

  subroutine lubksb(a,indx,b)
    !+ Does LU back substitution
    real(double), dimension(:,:),intent(IN)   :: a
    integer,dimension(:),  intent(IN)         :: indx
    real(double), dimension(:),  intent(INOUT):: b
    integer       :: i,n,ii,ll
    real(double)  :: summ
    n=size(indx)
    ii=0
    do i=1,n
       ll=indx(i)
       summ=b(ll)
       b(ll)=b(i)
       if (ii /= 0) then
          summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1))
       else if (summ /= 0.0) then
          ii=i
       endif
       b(i)=summ
    enddo
    do i=n,1,-1
       b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
    enddo
  end subroutine lubksb