lumatrix_ssl2.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  :: is, icon
  real(8)  :: vw(size(alu,2))
  integer  :: kk, nn
  
  kk = size(alu,1) 
    nn = size(alu,2)

  if ( kk > nn ) then
     call MessageNotify('W','ludecomp21', 'The first dimension is greater than the second')
  elseif( kk < nn ) then
     call MessageNotify('E','ludecomp21', 'The first dimension is less than the second')
  endif
  
  !" 行列のＬＵ分解（クラウト法）
  call dalu( alu, kk, nn, 0.0D0, kp, is, vw, 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  :: is, icon
  real(8)  :: vw(size(alu,3))
  integer  :: jj, kk, nn
  integer  :: j
  
  jj = size(alu,1) 
   kk = size(alu,2) 
    nn = size(alu,3)
  
  if ( kk > nn ) then
     call MessageNotify('W','ludecomp32', 'The second dimension is greater than the third')
  elseif( kk < nn ) then
     call MessageNotify('E','ludecomp32', 'The second dimension is less than the third')
  endif
 
 !" 行列のＬＵ分解（クラウト法）
!$omp parallel do private(vw,is,icon)
  do j=1,jj
     call dalu( alu(j,:,:), kk, nn, 0.0D0, kp(j,:), is, vw, 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 :: icon
  integer :: kk, nn

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

  if ( kk > nn ) then
     call MessageNotify('W','ludsolve211', 'The first dimension is greater than the second')
  elseif( kk < nn ) then
     call MessageNotify('E','lusolve211', 'The first dimension is less 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 dlux( lusolve211, alu, kk, nn, 1, kp, icon )
  if ( icon /= 0 ) call MessageNotify('E','lusolve211','LU solve 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 :: icon
  integer :: kk, nn
  integer :: j
  
  kk = size(alu,1) 
    nn = size(alu,2)

  if ( kk > nn ) then
     call MessageNotify('W','ludsolve212', 'The first dimension is greater than the second')
  elseif( kk < nn ) then
     call MessageNotify('E','lusolve212', 'The first dimension is less than 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 dlux( lusolve212(j,:), alu, kk, nn, 1, kp, icon )
     if ( icon /= 0 ) call MessageNotify('E','lusolve211','LU solve not succeeded.')
  end do
!$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  :: 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('W','lusovel322', 'The second dimension is greater than the third')
  elseif( kk < nn ) then
     call MessageNotify('E','lusolve322', 'The second dimension is less than the third')
  elseif( jj  /= size(b,1) .OR. nn /= size(b,2) ) then
     call MessageNotify('E','lusovle322', 'The input vector dimension does not match the matrix dimension')
  endif

  lusolve322 = b
 
!$omp parallel do private(icon)
  do j=1,jj
     call dlux( lusolve322(j,:), alu(j,:,:), kk, nn, 1, kp(j,:), icon )
     if ( icon /= 0 ) call MessageNotify('E','lusolve322','LU solve 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  :: icon
  integer  :: jj, kk, nn
  integer  :: i, j

  jj = size(alu,1) 
   kk = size(alu,2) 
    nn = size(alu,3)
  
  if ( kk > nn ) then
     call MessageNotify('W','lusovel322', 'The second dimension is greater than the third')
  elseif( kk < nn ) then
     call MessageNotify('E','lusolve322', 'The second dimension is less than the third')
  elseif( jj  /= size(b,2) .OR. nn /= size(b,3) ) then
     call MessageNotify('E','lusovle322', 'The input vector dimension does not match the matrix dimension')
  endif
  
  lusolve323 = b

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