subroutine ludecomp_prep_simple_many( mm, nn, mat, ms, me )

    integer , intent(in   ) :: mm, nn
    real(DP), intent(inout) :: mat( mm, nn, nn )
    integer , intent(in   ) :: ms, me


    !
    ! local variables
    !
    real(DP) :: lmat( ms:me, nn, nn )
    real(DP) :: ratio( ms:me )
    integer  :: i, j, k, m


    do j = 1, nn
       do i = 1, nn
          do m = ms, me
             lmat( m, i, j ) = 0.0d0
          end do
       end do
    end do
    do i = 1, nn
       do m = ms, me
          lmat( m, i, i ) = 1.0d0
       end do
    end do

    do k = 1, nn-1
       do i = k+1, nn
          do m = ms, me
             ratio( m ) = mat( m, i, k ) / mat( m, k, k )
          end do

          do j = 1, nn
             do m = ms, me
                mat( m, i, j ) = mat( m, i, j ) - mat( m, k, j ) * ratio( m )
             end do
          end do

          do m = ms, me
             lmat( m, i, k ) = ratio( m )
          end do
       end do

    end do

    !
    ! assemble into 1 matrix
    !
    do j = 1, nn
       do i = j+1, nn
          do m = ms, me
             mat( m, i, j ) = lmat( m, i, j )
          end do
       end do
    end do


  end subroutine ludecomp_prep_simple_many

  subroutine ludecomp_solve_simple_many( mm, nn, mat, vec, ms, me )

    integer , intent(in   ) :: mm, nn
    real(DP), intent(in   ) :: mat( mm, nn, nn )
    real(DP), intent(inout) :: vec( mm, nn )
    integer , intent(in   ) :: ms, me


    !
    ! local variables
    !
    real(DP) :: tmp( ms:me )
    integer  :: i, j, m


    ! solve matrix

    do i = 1+1, nn
       do m = ms, me
          tmp( m ) = 0.0d0
       end do
       do j = 1, i-1
          do m = ms, me
             tmp( m ) = tmp( m ) + mat( m, i, j ) * vec( m, j )
          end do
       end do
       do m = ms, me
          vec( m, i ) = vec( m, i ) - tmp( m )
       end do
    end do


    do i = nn, 1, -1
       do m = ms, me
          tmp( m ) = 0.0d0
       end do
       do j = i+1, nn
          do m = ms, me
             tmp( m ) = tmp( m ) + mat( m, i, j ) * vec( m, j )
          end do
       end do
       do m = ms, me
          vec( m, i ) = ( vec( m, i ) - tmp( m ) ) / mat( m, i, i )
       end do
    end do



  end subroutine ludecomp_solve_simple_many