!--
!----------------------------------------------------------------------
! Copyright (c) 2011 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  lumatrix_omp : 行列の LU 分解による線形連立方程式の解法(Vector版)
!
!履歴  1990/08/31(numaguti)
!      1999/11/10(takepiro) LUSOLM 追加
!      2002/01/20(takepiro) ispack-f90 用に抜きだし
!      2002/06/10(takepiro) ベクトル長問題対応 lusol2 を用意.
!      2008/10/29(takepiro) OpenMP 用に改造
!      2009/08/06(takepiro) LUMAK1 追加
!      2011/02/20 佐々木洋平 Fortran90 化
!      2016/04/11(takepiro) ベクトル長問題対応 lusol3 を用意.
!++
subroutine lumake(alu, kp, jdim, ndim)
  !
  ! * ndim x ndim の行列 jdim 個を一度に計算
  ! * LU 行列は入力に上書きされる
  !
  implicit none
  integer, intent(in)    :: jdim
  integer, intent(in)    :: ndim
  integer, intent(out)   :: kp(jdim, ndim)
  real(8), intent(inout) :: alu(jdim, ndim, ndim)
  integer :: j, k, m, n
  real(8) :: pivot, temp

  do k=1,ndim-1
    do j=1,jdim
!! pivot 選択
      pivot   = alu(j,k,k)
      kp(j,k) = k
      do m=k+1,ndim
        if ( abs(alu(j,m,k)) .gt. abs(pivot)) then
          pivot   = alu(j,m,k)
          kp(j,k) = m
        end if
      end do
      if ( kp(j,k) .ne. k ) then
        do n=1,ndim
          temp             = alu(j,k,n)
          alu(j,k,n)       = alu(j,kp(j,k),n)
          alu(j,kp(j,k),n) = temp
        end do
      end if

!! LU 分解
      do n=k+1,ndim
        alu(j,k,n) = alu(j,k,n)/pivot
        do m=k+1,ndim
          alu(j,m,n) = alu(j,m,n) - alu(j,m,k) * alu(j,k,n)
        end do
      end do
    end do
  end do
end subroutine lumake

subroutine lumak1(alu, kp, ndim)
  !
  ! * ndim x ndim の行列一個を計算
  ! * LU 行列は入力に上書き
  !
  implicit none
  integer, intent(in)    :: ndim
  real(8), intent(inout) :: alu(ndim,ndim)
  integer, intent(out)   :: kp(ndim)
  call lumake(alu, kp, 1, ndim)
end subroutine lumak1

subroutine lusolv(xv, alu, kp, idim, jdim, ndim)
  !
  ! * ndim x ndim 行列を jdim 個並べた連立方程式 AX = B を idim 個 の B
  !   について計算する.
  ! * 解は右辺入力ベクトルに上書きされる
  !
  implicit none
  integer, intent(in) :: idim
  integer, intent(in) :: jdim
  integer, intent(in) :: ndim
  real(8), intent(inout) :: xv(idim, jdim, ndim)
  real(8), intent(in) :: alu(jdim, ndim, ndim)
  integer, intent(in) :: kp(jdim, ndim)
  integer :: i, j, k, n, nn
  real(8) :: temp

  do k=1, ndim-1
    do i=1, idim
      do j=1,jdim
        if ( kp(j,k) .ne. k ) then
          temp            = xv(i,j,k)
          xv(i,j,k)       = xv(i,j,kp(j,k))
          xv(i,j,kp(j,k)) = temp
        end if
      end do
    end do
  end do

  do n=1, ndim
    do i=1, idim
      do j=1, jdim
        xv(i,j,n) = xv(i,j,n)/alu(j,n,n)
      end do
    end do
    do nn=n+1,ndim
      do i=1, idim
        do j=1, jdim
          xv(i,j,nn) = xv(i,j,nn) - xv(i,j,n) * alu(j,nn,n)
        end do
      end do
    end do
  end do
  do k=ndim-1, 1, -1
    do n=k+1,ndim
      do i=1,idim
        do j=1,jdim
          xv(i,j,k) = xv(i,j,k) - xv(i,j,n) * alu(j,k,n)
        end do
      end do
    end do
  end do
end subroutine lusolv

subroutine lusol3(xv, alu, kp, jdim, ndim)
  !
  ! * ndim x ndim 行列を jdim 個並べた連立方程式 AX = B を jdim 個 の B
  !   について計算する.
  ! * 解は右辺入力ベクトルに上書きされる
  !
  implicit none
  integer, intent(in) :: jdim
  integer, intent(in) :: ndim
  real(8), intent(inout) :: xv(jdim, ndim)
  real(8), intent(in) :: alu(jdim, ndim, ndim)
  integer, intent(in) :: kp(jdim, ndim)
  integer :: j, k, n, nn
  real(8) :: temp

  do k=1, ndim-1
     do j=1,jdim
        if ( kp(j,k) .ne. k ) then
           temp            = xv(j,k)
           xv(j,k)       = xv(j,kp(j,k))
           xv(j,kp(j,k)) = temp
        end if
     end do
  end do

  do n=1, ndim
     do j=1, jdim
        xv(j,n) = xv(j,n)/alu(j,n,n)
     end do
     do nn=n+1,ndim
        do j=1, jdim
          xv(j,nn) = xv(j,nn) - xv(j,n) * alu(j,nn,n)
        end do
     end do
  end do
  do k=ndim-1, 1, -1
    do n=k+1,ndim
       do j=1,jdim
          xv(j,k) = xv(j,k) - xv(j,n) * alu(j,k,n)
       end do
    end do
  end do
end subroutine lusol3

subroutine lusol2(xv, alu, kp, idim, ndim)
  !
  ! * ndim x ndim 型行列の連立方程式 A X = B を idim 個の B に対して計算.
  ! * 解は右辺の入力ベクトルに上書きされる.
  ! * LUSOLV の JDIM = 1 に相当. JDIM=1 には際のベクトル長が短くなるため,
  !   このルーチンを用意している
  !
  implicit none
  integer, intent(in) :: idim
  integer, intent(in) :: ndim
  real(8), intent(inout) :: xv(idim, ndim)
  real(8), intent(in) :: alu(ndim,ndim)
  integer, intent(in) :: kp(ndim)
  integer :: i, k, n, nn
  real(8) :: temp

  do k = 1, ndim-1
    do i = 1, idim
      if ( kp ( k ) .ne. k ) then
        temp           = xv ( i,k )
        xv ( i,k )     = xv ( i,kp(k) )
        xv ( i,kp(k) ) = temp
      endif
    end do
  end do

  do n = 1, ndim
    do i = 1, idim
      xv ( i,n ) = xv ( i,n ) / alu ( n,n )
      do nn = n+1, ndim
        xv ( i,nn ) = xv ( i,nn ) - xv ( i,n ) * alu ( nn,n )
      end do
    end do
  end do

  do k = ndim-1, 1, -1
    do n = k+1, ndim
      do i = 1, idim
        xv ( i,k ) = xv ( i,k ) - xv ( i,n ) * alu ( k,n )
      end do
    end do
  end do

end subroutine lusol2

subroutine lusolm(xv, alu, kp, jmtx, idim, jdim, ndim)
  !
  ! * jdim 個の ndim x ndim 型行列について
  !   連立方程式 A X = B を idim 個の B に対して計算.
  ! * 解は右辺の入力ベクトルに上書きされる.
  !
  implicit none
  integer, intent(in)    :: idim
  integer, intent(in)    :: jdim
  integer, intent(in)    :: ndim
  real(8), intent(inout) :: xv(idim, ndim)
  real(8), intent(in)    :: alu(jdim, ndim, ndim)
  integer, intent(in)    :: kp(jdim, ndim)
  integer, intent(in)    :: jmtx(idim)
  integer :: i, k, n, nn
  real(8) :: temp

  do k = 1, ndim-1
    do i = 1, idim
      if ( kp ( jmtx(i),k ) .ne. k ) then
        temp             = xv ( i,k )
        xv ( i,k )       = xv ( i,kp(jmtx(i),k) )
        xv ( i,kp(jmtx(i),k) ) = temp
      endif
    end do
  end do
  do  n = 1, ndim
    do i = 1, idim
      xv ( i,n ) = xv ( i,n ) / alu ( jmtx(i),n,n )
    end do
    do  nn = n+1, ndim
      do i = 1, idim
        xv ( i,nn ) = xv ( i,nn ) - xv ( i,n ) * alu ( jmtx(i),nn,n )
      end do
    end do
  end do
  do  k = ndim-1, 1, -1
    do  n = k+1, ndim
      do  i = 1, idim
        xv ( i,k ) = xv ( i,k ) - xv ( i,n ) * alu ( jmtx(i),k,n )
      end do
    end do
  end do
end subroutine lusolm