LDU 分解するプログラムをつくろう。
program ldu
implicit none
integer :: fr = 15
integer a(0:2,0:2), b(0:2), i, j
real(8) r(0:2,0:2), br(0:2)
real(8) :: d(0:2,0:2) = 0.0d0 , l(0:2,0:2) = 0.0d0 ,
& u(0:2,0:2) = 0.0d0
open(fr, file = 'test.dat')
do i = 0, 2
read(fr, *) a(i,0:2), b(i)
end do
close(fr)
do i = 0, 2
br(i) = dble(b(i))
do j = 0, 2
r(i,j) = dble(a(i,j))
end do
end do
do i = 0, 2
do j = 0, 2
if (j == i) then
d(i,j) = r(i,j)
else if (i > j) then
l(i,j) = r(i,j)
else
u(i,j) = r(i,j)
end if
end do
end do
do i = 0, 2
write(*,*) d(i,:)
end do
do i = 0, 2
write(*,*) l(i,:)
end do
do i = 0, 2
write(*,*) u(i,:)
end do
end program ldu
長い宣言文になったため、6 列目以前に & をつけ、 前の行からの続きであることを示している。
b≠ 0 かつ aii≠ 0、ならば、 T = (D - L)-1 U から作られる行列が収束行列、すなわち k が無限大で T は限りなく零行列に近づく。 このとき、行列 T の 固有値の最大値を max|λi| とすると、
ρ(T) = max|λi| > 1 (i > i > n)
が成り立つ。A と A からなる対角行列が正定値行列の場合は、収束解を得ることができる。
module matmod
implicit none
contains
function sddm(ra) result(sd)
real(8), intent(in) :: ra(0:,0:)
real(8) sd(0:size(ra,1)-1,0:size(ra,2)-1), v(0:size(ra,1)-1)
integer i, j, k, l
k = size(ra(1,:))
v(0:k-1) = 0.0d0
sd(:,:) = ra(:,:)
do i = 0, k - 1
do j = 0, i - 1
v(i) = v(i) + abs(sd(i, j))
end do
do j = i + 1, k
v(i) = v(i) + sd(i,j)
end do
if (abs(ra(i,i)) > v(i)) then
sd(0:k - 1,0:k - 1) = ra(0:k - 1,0:k - 1)
else
sd(0:k - 1,0:k - 1) = 0.0d0
end if
end do
end function sddm
function trace(ra) result(tr)
real(8), intent(in) :: ra(0:,0:)
real(8) tr
integer i, k
k = size(ra,1)
tr = 0.0d0
do i = 0, k
if (ra(i,i) == 0.0d0) then
stop 'diag i is 0'
end if
tr = tr + ra(i,i)
end do
end function trace
function prev(ra) result(pr)
real(8), intent(in) :: ra(0: ,0:)
real(8) pr(0 : size(ra, 1) - 1, 0:size(ra, 2) - 1)
real(8) pv(0 : size(ra, 2) - 1), pv2(0 : size(ra, 2) - 1)
integer i, j, k, l
l = size(ra, 2)
k = size(ra, 1)
pr = ra
do i = 0, k - 2
j = i + 1
if (ra(i, i) == 0.0d0) then
pv(:) = ra(i, :)
pv2(:) = ra(j, :)
pr(i, :) = pv2(:)
pr(i + 1, :) = pv(:)
end if
end do
end function prev
subroutine gauss_seidel(ra, rb, rx, maxit, eps)
real(8), intent(in) :: ra(0: ,0:), rb(0:), eps
integer, intent(in) :: maxit
real(8), intent(out) :: rx(0: size(ra, 1) - 1)
real(8) ry(0:size(ra,1) - 1)
integer i, j, k, ct
k = size(ra,2)
rx(:) = 0.0d0
do ct = 1, maxit
do i = 0, k - 1
rx(i) = -1.0d0 * dot_product(ra(i, 0: i - 1), rx(0 : i - 1))
rx(i) = rx(i)
& -1.0d0 *
& dot_product(ra(i, i + 1 : k -1), rx(i + 1 : k -1))
rx(i) = rx(i) + rb(i)
rx(i) = rx(i) / ra(i, i)
ry(:) = rb(:) - matmul(ra,rx)
end do
if (dot_product(ry,ry) <= eps) then
write(*,*) ct
exit
end if
end do
end subroutine gauss_seidel
end module matmod
program gs
use matmod
implicit none
integer :: n = 0, fr = 15, maxit = 10000
real(8) :: eps = 1.0d-6
integer l
integer io
real(8), allocatable :: a(:,:), b(:), x(:), sum(:)
integer i, j
open(fr, file = 'test.dat', action = 'read')
do
read(fr, *, iostat = io) l
if (io < 0) then
exit
else
n = n + 1
end if
end do
close(fr)
allocate (a(0:n-1,0:n-1))
allocate (b(0:n-1))
allocate (x(0:n-1))
allocate (sum(0:n-1))
open(fr, file = 'test.dat', action = 'read')
do i = 0, n-1
read(fr, *, iostat = io) a(i, 0 : n - 1), b(i)
a(i,0:n-1) = dble(a(i, 0: n - 1))
b(i) = dble(b(i))
end do
close(fr)
a = sddm(a)
do i = 0, n-1
write(*,'(100e12.4)') a(i, 0 : n - 1), b(i)
end do
a = prev(a)
write(*,*) trace(a)
do i = 0, n-1
write(*,'(100e12.4)') a(i, 0 : n - 1), b(i)
end do
call gauss_seidel(a, b, x, maxit, eps)
do i = 0, n-1
write(*,'(100e12.4)') x(i)
end do
end program gs
Gauss--Seidel 法のプログラムのおよその形はこのようになる。