Path: | libsrc/lumatrix/lumatrix_vec.f90 |
Last Update: | Mon Aug 19 17:49:30 +0900 2013 |
Subroutine : | |
alu(ndim,ndim) : | real(8), intent(inout) |
kp(ndim) : | integer, intent(out) |
ndim : | integer, intent(in) |
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 : | |
alu(jdim, ndim, ndim) : | real(8), intent(inout) |
kp(jdim, ndim) : | integer, intent(out) |
jdim : | integer, intent(in) |
ndim : | integer, intent(in) |
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 : | |
xv(idim, ndim) : | real(8), intent(inout) |
alu(ndim,ndim) : | real(8), intent(in) |
kp(ndim) : | integer, intent(in) |
idim : | integer, intent(in) |
ndim : | integer, intent(in) |
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 : | |
xv(idim, ndim) : | real(8), intent(inout) |
alu(jdim, ndim, ndim) : | real(8), intent(in) |
kp(jdim, ndim) : | integer, intent(in) |
jmtx(idim) : | integer, intent(in) |
idim : | integer, intent(in) |
jdim : | integer, intent(in) |
ndim : | integer, intent(in) |
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
Subroutine : | |
xv(idim, jdim, ndim) : | real(8), intent(inout) |
alu(jdim, ndim, ndim) : | real(8), intent(in) |
kp(jdim, ndim) : | integer, intent(in) |
idim : | integer, intent(in) |
jdim : | integer, intent(in) |
ndim : | integer, intent(in) |
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