module subroutine の場合いずれも intent(out) を明示し、 配列に値を渡して戻してもよいが、 渡す際にも答えを得る配列を用意しておかなくてはならない。 inner subroutine で答えを得たい場合には、module function の方が便利である。
以前の FORTRAN 77 ではスカラ関数のみであったが、 配列を返すことができるようになった。 また、部分配列の組み込み演算も可能である。
function 関数(仮引数) result(返す変数) 型, intent(in) :: 仮引数 型 内部変数 返す変数への代入文 end function 関数
module op
implicit none
contains
function unit(mat) result(e)
real(8), intent(in) :: mat(0:,0:)
real(8), allocatable :: e(:,:)
integer j, l, k
k = size(mat,1)
l = size(mat,2)
allocate(e(0:k-1,0:l-1))
e(0:k-1,0:l-1) = 0.0d0
do j = 0, k-1
e(j,j) = 1.0d0
end do
end function unit
function pci(c,i,mat) result(p)
real(8), intent(in) :: c, mat(0:,0:)
real(8), allocatable :: p(:,:)
integer, intent(in) :: i
integer s, t, k, l
k = size(mat,1)
l = size(mat,2)
allocate(p(0:k-1,0:l-1))
p = unit(p)
p(i,i) = c * p(i,i)
p = matmul(transpose(p),mat)
end function pci
function qcij(c,i,j,mat) result(q)
real(8), intent(in) :: c, mat(0:,0:)
real(8), allocatable :: q(:,:)
integer, intent(in) :: i,j
integer s, t, k, l
k = size(mat,1)
l = size(mat,2)
allocate(q(0:k-1,0:l-1))
q = unit(q)
q(i,j) = c * 1.0d0
q = matmul(transpose(q),mat)
end function qcij
function rij(i,j,mat) result(r)
real(8), intent(in) :: mat(0:,0:)
real(8), allocatable :: r(:,:), vec(:)
integer, intent(in) :: i,j
integer s, t, k, l
k = size(mat,1)
l = size(mat,2)
allocate(r(0:k-1,0:l-1))
allocate(vec(0:l-1))
r = unit(r)
vec(0:l-1) = r(i,:)
r(i,0:l-1) = r(j,0:l-1)
r(j,0:l-1) = vec(0:l-1)
r = matmul(transpose(r),mat)
end function rij
end module op
program gauss
use op
implicit none
integer :: row = 4, col = 5
integer m
integer i, j, p, q
real(8) c
real(8), allocatable :: mat(:,:), cal(:,:)
real(8), allocatable :: b(:), x(:)
allocate(mat(0:row-1,0:col-1))
allocate(cal(0:row-1,0:col-1))
allocate(b(0:row-1))
allocate(x(0:row-1))
mat(0,0:col-1) = (/ 2.0d0, 3.0d0, 1.0d0, -3.0d0, 1.0d0/)
mat(1,0:col-1) = (/ -1.0d0, 2.0d0, 2.0d0, 4.0d0, 6.0d0/)
mat(2,0:col-1) = (/ 4.0d0, 1.0d0, -3.0d0, 5.0d0, 3.0d0/)
mat(3,0:col-1) = (/ 5.0d0, -4.0d0, -4.0d0, 1.0d0, 3.0d0/)
do m = 0, row - 1
write(*,'(100e12.4)') mat(m, 0:col-1)
end do
write(*,*)
cal = mat
do i = 0, row - 1
if (cal(i,i) == 1.0d0) stop 'coeff is nan'
do j = i + 1, row - 1
if (abs(cal(i,j)) > cal(i,i)) then
cal = rij(i,j,cal)
end if
end do
c = 1.0d0 / cal(i,i)
cal = pci(c, i, cal)
do j = i+1, row -1
cal = qcij(-1.0d0 * cal(j,i),i,j,cal)
end do
end do
write(*,*)
do m = 0, row - 1
write(*,'(100e12.4)') cal(m, 0:col-1)
end do
write(*,*)
if (row + 1 == col) then
do m = 0, row - 1
p = col - 2 - m
b(p) = cal(p,col - 1)
x(p) = b(p)
do q = p + 1, col - 2
x(p) = x(p) - 1.0d0 * cal(p,q) * x(q)
end do
end do
write(*,'(100e12.4)') x
end if
end program gauss
pivot の代入が本文内で短く済むため、わざわざ呼び出さなかったりしているが、 組み込み関数がいくつか定義されていたり、 Ruby で先に作った Gauss の消去法によるプログラムを焼き直したものである。