!******************************************************************************
! subroutine ludecomp
!******************************************************************************
!     Arrangement of array elements
!----------------------------------------------------------------------
!  A x = b
!          a(1,1)  a(1,2)  a(1,3)  a(1,4) ...  a(1,n)  |   x(1)  |   b(1)
!          a(2,1)  a(2,2)  a(2,3)  a(2,4) ...  a(2,n)  |   x(2)  |   b(2)
!          a(3,1)  a(3,2)  a(3,3)  a(3,4) ...  a(3,n)  |   x(3)  |   b(3)
!          a(4,1)  a(4,2)  a(4,3)  a(4,4) ...  a(4,n)  |   x(4)  |   b(4)
!          a(5,1)  a(5,2)  a(5,3)  a(5,4) ...  a(5,n)  |   x(5)  |   b(5)
!          ...     ...     ...     ...    ...  ...     |   ...   |   ...
!          a(n,1)  a(n,2)  a(n,3)  a(n,4) ...  a(n,n)  |   x(n)  |   b(n)
!----------------------------------------------------------------------
!  A = L * U
!          a(1,1)  a(1,2)  a(1,3)  a(1,4) ...  a(1,n)
!          a(2,1)  a(2,2)  a(2,3)  a(2,4) ...  a(2,n)
!          a(3,1)  a(3,2)  a(3,3)  a(3,4) ...  a(3,n)
!          a(4,1)  a(4,2)  a(4,3)  a(4,4) ...  a(4,n)
!          a(5,1)  a(5,2)  a(5,3)  a(5,4) ...  a(5,n)
!          ...     ...     ...     ...    ...  ...
!          a(n,1)  a(n,2)  a(n,3)  a(n,4) ...  a(n,n)
!----------------------------------------------------------------------
!  matrix L
!          1       0       0       0      ...  0
!          a(2,1)  1       0       0      ...  0
!          a(3,1)  a(3,2)  1       0      ...  0
!          a(4,1)  a(4,2)  a(4,3)  1      ...  0
!          a(5,1)  a(5,2)  a(5,3)  a(5,4) ...  0
!          ...     ...     ...     ...    ...  ...
!          a(n,1)  a(n,2)  a(n,3)  a(n,4) ...  a(n,n)
!----------------------------------------------------------------------
!  matrix L
!          a(1,1)  a(1,2)  a(1,3)  a(1,4) ...  a(1,n)
!          0       a(2,2)  a(2,3)  a(2,4) ...  a(2,n)
!          0       0       a(3,3)  a(3,4) ...  a(3,n)
!          0       0       0       a(4,4) ...  a(4,n)
!          0       0       0       0      ...  a(5,n)
!          ...     ...     ...     ...    ...  ...
!          0       0       0       0      ...  a(n,n)
!******************************************************************************

  subroutine ludecomp_prep_simple( nn, mat )

    use vtype_module

    implicit none

    integer(i4b), intent(in   ) :: nn
    real(dp)    , intent(inout) :: mat( nn, nn )


    !
    ! local variables
    !
    real(dp)     :: lmat( nn, nn )
!!$    real(dp)     :: mat_prod( nn, nn )
    real(dp)     :: ratio
    integer(i4b) :: i, j, k


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

    do k = 1, nn-1
       do i = k+1, nn
          ratio = mat( i, k ) / mat( k, k )

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

          lmat( i, k ) = ratio
       end do


!!$       write( 6, * ) '-----', k
!!$    do i = 1, nn
!!$       write( 6, * ) i, '|', ( real( lmat(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )
!!$    do i = 1, nn
!!$       write( 6, * ) i, '|', ( real(  mat(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )

    end do


!!$    write( 6, * ) 'ludecomp_prep : LU'
!!$    do i = 1, nn
!!$       write( 6, * ) i, '|', ( real( lmat(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )
!!$    do i = 1, nn
!!$       write( 6, * ) i, '|', ( real(  mat(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )
!!$
!!$
!!$    do j = 1, nn
!!$       do i = 1, nn
!!$          mat_prod( i, j ) = 0.0d0
!!$       end do
!!$    end do
!!$    do j = 1, nn
!!$       do i = 1, nn
!!$          do k = 1, nn
!!$             mat_prod( i, j ) = mat_prod( i, j ) + lmat( i, k ) * mat( k, j )
!!$          end do
!!$       end do
!!$    end do
!!$
!!$    write( 6, * ) 'ludecomp_prep : product'
!!$    do i = 1, nn
!!$       write( 6, * ) i, ( real( mat_prod(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )


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


  end subroutine ludecomp_prep_simple

  !**************************************************************************************

  subroutine ludecomp_solve_simple( nn, mat, vec )

    use vtype_module

    implicit none

    integer(i4b), intent(in   ) :: nn
    real(dp)    , intent(inout) :: mat( nn, nn )
    real(dp)    , intent(inout) :: vec( nn )


    !
    ! local variables
    !
    real(dp)     :: tmp
    real(dp)     :: lmat( nn, nn )
!!$    real(dp)     :: mat_prod( nn, nn )
    integer(i4b) :: i, j


    ! decomposition into 2 matrices, L and U

    do j = 1, nn
       do i = 1, j-1
          lmat( i, j ) = 0.0d0
       end do
       do i = j, j
          lmat( i, j ) = 1.0d0
       end do
       do i = j+1, nn
          lmat( i, j ) = mat( i, j )
       end do
    end do
    do j = 1, nn
       do i = j+1, nn
          mat( i, j ) = 0.0d0
       end do
    end do


!!$    write( 6, * ) 'ludecomp_solve : LU'
!!$    do i = 1, nn
!!$       write( 6, * ) i, '|', ( real( lmat(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )
!!$    do i = 1, nn
!!$       write( 6, * ) i, '|', ( real(  mat(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )
!!$
!!$
!!$    do j = 1, nn
!!$       do i = 1, nn
!!$          mat_prod( i, j ) = 0.0d0
!!$       end do
!!$    end do
!!$    do j = 1, nn
!!$       do i = 1, nn
!!$          do k = 1, nn
!!$             mat_prod( i, j ) = mat_prod( i, j ) + lmat( i, k ) * mat( k, j )
!!$          end do
!!$       end do
!!$    end do
!!$
!!$    write( 6, * ) 'ludecomp_solve : product'
!!$    do i = 1, nn
!!$       write( 6, * ) i, ( real( mat_prod(i,j) ), j=1, nn )
!!$    end do
!!$    write( 6, * )


    ! solve matrix

    do i = 1+1, nn
       tmp = 0.0d0
       do j = 1, i-1
          tmp = tmp + lmat( i, j ) * vec( j )
       end do
       vec( i ) = vec( i ) - tmp
    end do


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



  end subroutine ludecomp_solve_simple

  !**************************************************************************************

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

    use vtype_module

    implicit none

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


    !
    ! local variables
    !
    real(dp)     :: lmat( ms:me, nn, nn )
    real(dp)     :: ratio( ms:me )
    integer(i4b) :: 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 )

    use vtype_module

    implicit none

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


    !
    ! local variables
    !
    real(dp)     :: tmp( ms:me )
    integer(i4b) :: 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

  !**************************************************************************************
