lumatrix_lapack.f90

subroutine ludecomp21(alu,kp)
  !
  ! ALU(NDIM,NDIM), KP(NDIM)
  ! NDIM x NDIM の行列を LU 分解.
  ! ＬＵ行列は 入力行列に上書きされる.
  !
  use dc_message

  real(8), intent(inout) :: alu(:,:)                  ! 入力／ＬＵ行列
  integer, intent(out)   :: kp(size(alu,2))           ! ピボット
  integer  :: icon
  integer  :: kk, nn

  kk = size(alu,1) 
    nn = size(alu,2)

  if ( kk /= nn ) then
     call MessageNotify('E','ludecomp21', 'The first dimension differs from the second')
  endif

  !" 行列のＬＵ分解
  call dgetrf( nn, nn, alu, nn, kp, icon )
  if ( icon /= 0 ) call MessageNotify('E','ludecomp21','LU decompostion not succeeded.')

end subroutine ludecomp21

subroutine ludecomp32(alu,kp)
  !
  ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM)
  ! NDIM x NDIM の行列 JDIM 個を一度に LU 分解.
  ! ＬＵ行列は 入力行列に上書きされる.
  !
  use dc_message

  real(8), intent(inout) :: alu(:,:,:)                  ! 入力／ＬＵ行列
  integer, intent(out)   :: kp(size(alu,1),size(alu,3)) ! ピボット

  integer :: icon
  integer :: jj, kk, nn
  integer :: j

  jj = size(alu,1) 
   kk = size(alu,2) 
    nn = size(alu,3)

  if ( kk /= nn ) then
     call MessageNotify('E','ludecomp32', 'The second dimension differs from the third')
  endif

 !" 行列のＬＵ分解（クラウト法）
 !$omp parallel do private(icon)
  do j=1,jj
     call dgetrf( nn, nn, alu(j,:,:), nn, kp(j,:), icon )
     if ( icon /= 0 ) call MessageNotify('E','ludecomp32','LU decompostion not succeeded.')
  end do
!$omp end parallel do

end subroutine ludecomp32

function lusolve211(alu,kp,b)
  !
  ! ALU(NDIM,NDIM), KP(NDIM), B(NDIM)
  ! NDIM x NDIM 型行列の連立方程式
  ! A X = B を 1 個の B に対して計算する.
  !
  use dc_message

  real(8), intent(in)  :: alu(:,:)              ! 入力／ＬＵ行列
  integer, intent(in)  :: kp(:)                 ! ピボット
  real(8), intent(in)  :: b(:)                  ! 右辺ベクトル

  real(8) :: lusolve211(size(b))                ! 解
  integer :: kk, nn
  integer :: icon

  kk = size(alu,1) 
    nn = size(alu,2)

  if ( kk /= nn ) then
     call MessageNotify('E','ludsolve211', 'The first dimension is greater than the second')
  elseif ( nn /= size(b)) then
     call MessageNotify('E','lusolve211', 'The input vector length differs from the second dimension of the matrix')
  endif

  lusolve211 = b
  call dgetrs( 'N', nn, 1, alu, nn, kp, lusolve211, nn, icon )
  if ( icon /= 0 ) call MessageNotify('E','lusplve211','LU decompostion not succeeded.')

end function lusolve211

function lusolve212(alu,kp,b)
  !
  ! ALU(NDIM,NDIM), KP(NDIM), B(JDIM,NDIM)
  ! NDIM x NDIM 型行列の連立方程式
  ! A X = B を JDIM 個の B に対して計算する.
  !
  use dc_message

  real(8), intent(in)  :: alu(:,:)              ! 入力／ＬＵ行列
  integer, intent(in)  :: kp(:)                 ! ピボット
  real(8), intent(in)  :: b(:,:)                ! 右辺ベクトル

  real(8) :: lusolve212(size(b,1),size(b,2))       ! 解

  integer :: kk, nn
  integer :: j
  integer :: icon

  kk = size(alu,1) 
    nn = size(alu,2)

  if ( kk /= nn ) then
     call MessageNotify('E','ludsolve212', 'The first dimension differs from the second')
  elseif ( nn /= size(b,2)) then
     call MessageNotify('E','lusolvep212', 'The vector length differs from the second dimension of the matrix')
  endif

  lusolve212 = b

!$omp parallel do private(icon)
  do j=1,size(b,1)
     call dgetrs( 'N', nn, 1, alu, nn, kp, lusolve212(j,:), nn, icon )
     if ( icon /= 0 ) call MessageNotify('E','lusplve212','LU decompostion not succeeded.')
  enddo
!$omp end parallel do

end function lusolve212

function lusolve322(alu,kp,b)
  !
  ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(JDIM,NDIM)
  ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
  ! A X = B をひとつの B の並びに対して計算する.
  !
  use dc_message

  real(8), intent(in)  :: alu(:,:,:)                   ! 入力／ＬＵ行列
  integer, intent(in)  :: kp(:,:)                      ! ピボット
  real(8), intent(in)  :: b(:,:)                       ! 右辺ベクトル

  real(8) :: lusolve322(size(b,1),size(b,2))             ! 解

  integer  :: jj, kk, nn
  integer  :: j
  integer :: icon

  jj = size(alu,1) 
   kk = size(alu,2) 
    nn = size(alu,3)

  if ( kk /= nn ) then
     call MessageNotify('E','ludsolve322', 'The second dimension differs from the third')
  elseif ( nn /= size(b,2)) then
     call MessageNotify('E','lusolvep322', 'The vector length differs from the second dimension of the matrix')
  endif

  lusolve322 = b

!$omp parallel do private(icon)
  do j=1,jj
     call dgetrs( 'N', nn, 1, alu(j,:,:), nn, kp(j,:), lusolve322(j,:), nn, icon )
     if ( icon /= 0 ) call MessageNotify('E','lusplve322','LU decompostion not succeeded.')
  end do
!$omp end parallel do

end function lusolve322

function lusolve323(alu,kp,b)
  !
  ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(IDIM,JDIM,NDIM)
  ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
  ! A X = B を IDIM 個の B に対して計算する.
  !
  use dc_message

  real(8), intent(in)  :: alu(:,:,:)                   ! 入力／ＬＵ行列
  integer, intent(in)  :: kp(:,:)                      ! ピボット
  real(8), intent(in)  :: b(:,:,:)                     ! 右辺ベクトル

  real(8) :: lusolve323(size(b,1),size(b,2),size(b,3)) ! 解

  integer  :: jj, kk, nn
  integer  :: i, j
  integer :: icon

  jj = size(alu,1) 
   kk = size(alu,2) 
    nn = size(alu,3)


  if ( kk /= nn ) then
     call MessageNotify('E','ludsolve323', 'The first dimension differs from the second')
  elseif ( nn /= size(b,3)) then
     call MessageNotify('E','lusolvep323', 'The vector length differs from the second dimension of the matrix')
  endif

  lusolve323 = b

!$omp parallel do private(icon)
  do i=1,size(b,1)
     do j=1,jj
        call dgetrs( 'N', nn, 1, alu(j,:,:), nn, kp(j,:), lusolve323(i,j,:), nn, icon )
     if ( icon /= 0 ) call MessageNotify('E','lusplve323','LU decompostion not succeeded.')
     enddo
  enddo
!$omp end parallel do

end function lusolve323