SUBROUTINE lubksb_sp(a,indx,b)
  use mo_nrtype
  use mo_nrutil, ONLY : assert_eq
  IMPLICIT NONE
  REAL(SP), DIMENSION(:,:), INTENT(IN) :: a
  INTEGER(I4), DIMENSION(:), INTENT(IN) :: indx
  REAL(SP), DIMENSION(:), INTENT(INOUT) :: b
  INTEGER(I4) :: i,n,ii,ll
  REAL(SP) :: summ
  n=assert_eq(size(a,1),size(a,2),size(indx),'lubksb')
  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
     end if
     b(i)=summ
  end do
  do i=n,1,-1
     b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
  end do
END SUBROUTINE lubksb_sp
SUBROUTINE lubksb_dp(a,indx,b)
  use mo_nrtype
  use mo_nrutil, ONLY : assert_eq
  IMPLICIT NONE
  REAL(DP), DIMENSION(:,:), INTENT(IN) :: a
  INTEGER(I4), DIMENSION(:), INTENT(IN) :: indx
  REAL(DP), DIMENSION(:), INTENT(INOUT) :: b
  INTEGER(I4) :: i,n,ii,ll
  REAL(DP) :: summ
  n=assert_eq(size(a,1),size(a,2),size(indx),'lubksb')
  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
     end if
     b(i)=summ
  end do
  do i=n,1,-1
     b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
  end do
END SUBROUTINE lubksb_dp
