Fortran - LU 分解を用いた連立1次方程式の解法!

更新日時:


連立1次方程式を LU 分解を用いて解くアルゴリズムを Fortran 95 で実装してみました。
(使用する LU 分解法は「外積形式ガウス法(outer-product form)」)

前回 Ruby で同じことをしました。

0. 前提条件

  • LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
  • GCC 6.3.0 (GFortran 6.3.0) でのコンパイルを想定。

1. LU 分解(外積形式ガウス法)について

当ブログ過去記事を参照。

Ruby - LU 分解(外積形式ガウス法(outer-product form))!

2. Fortran ソースコードの作成

  • 本来、 L と U の2つの行列に分けるものだが1つの行列にまとめている。(実際に LU 分解を使用する際に L と U を意識して取り扱えばよいだけなので)
    また、 行列 L の対角成分が 1 であることを想定。

File: sle_lu.f95

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
!************************************************************
! 連立1次方程式の解法 ( LU 分解(外積形式ガウス法) )
!
!   date          name            version
!   2019.03.13    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!************************************************************
!
module const
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const

module lu
  use const
  implicit none
  private
  public :: decompose, solve

contains
  ! LU 分解
  ! * 外積形式ガウス法(outer-product form)
  !
  ! :param(inout) real(8) a(:,:): 行列
  subroutine decompose(a)
    implicit none
    real(DP), intent(inout) :: a(:, :)
    integer(SP) :: i, j, k, n
    real(DP)    :: tmp

    n = int(sqrt(real(size(a))))
    do k = 1, n
      if (a(1, 1) == 0.0_DP) then
        print *, "Can't divide by 0 ..."
        stop
      end if
      tmp = 1.0_DP / a(k, k)
      a(k+1:n, k) = a(k+1:n, k) * tmp
      do j = k + 1, n
        tmp = a(k, j)
        a(k+1:n, j) = a(k+1:n, j) - a(k+1:n, k) * tmp
      end do
    end do
  end subroutine decompose

  ! 連立方程式を解く
  !
  ! :param(in)  real(8) a(n,n): 行列 A
  ! :param(in)  real(8)   b(n): ベクトル b
  ! :param(out) real(8)   x(n): 解ベクトル x
  subroutine solve(a, b, x)
    implicit none
    real(DP),    intent(in)  :: a(:, :), b(:)
    real(DP),    intent(out) :: x(:)
    integer(SP) :: n, i, j
    real(DP)    :: s
    real(DP), allocatable :: y(:)

    ! 元数
    n = size(x)

    ! 配列 y のメモリ確保
    allocate(y(n))

    ! 前進代入
    ! * Ly = b から y を計算
    do i = 1, n
      s = sum((/(a(i, j) * y(j), j=1,i)/))
      y(i) = b(i) - s
    end do

    ! 後退代入
    ! * Ux = y から x を計算
    do i = n, 1, -1
      s = sum((/(a(i, j) * x(j), j=i+1,n)/))
      x(i) = (y(i) - s) / a(i, i)
    end do

    ! 配列 y のメモリ解放
    deallocate(y)
  end subroutine solve
end module lu

program sle_lu
  use const
  use lu
  implicit none
  character(9), parameter :: F_INP = "input.txt"
  integer(SP),  parameter :: UID   = 10
  integer(SP)   :: n, i                  ! 元数、ループインデックス
  real(DP), allocatable :: a(:,:), b(:)  ! 係数配列
  real(DP), allocatable :: x(:)          ! 解配列

  ! IN ファイル OPEN
  open(UID, file = F_INP, status = "old")

  ! 元数取得
  read(UID, *) n
  if (n < 1) stop
  print '("n = ", I0)', n

  ! 配列メモリ確保
  allocate(a(n, n))
  allocate(b(n))
  allocate(x(n))

  ! 行列 A 取得
  do i = 1, n
    read(UID, *) a(i,:)
  end do
  print '(A)', "A ="
  call display_mtx(a)

  ! ベクトル B 取得
  read(UID, *) b(:)
  print '(A)', "b ="
  call display_vec(b)

  ! IN ファイル CLOSE
  close (UID)

  ! 行列 A の LU 分解
  call decompose(a)
  !print '(A)', "(LU) ="
  !call display_mtx(a)

  ! 連立方程式を解く
  call solve(a, b, x)
  print '(A)', "x ="
  call display_vec(x)

  ! 配列メモリ解放
  deallocate(a)
  deallocate(b)
  deallocate(x)

contains
  subroutine display_mtx(m)
    implicit none
    real(DP), intent(in) :: m(:, :)
    integer(SP) :: n
    character(20) :: f  ! 書式文字列

    n = size(m(1, :))
    write (f, '("(", I0, "(F10.2)", ")")') n + 1
    do i = 1, n
      print f, m(i,:)
    end do
  end subroutine display_mtx

  subroutine display_vec(v)
    implicit none
    real(DP), intent(in) :: v(:)
    integer(SP) :: n
    character(20) :: f  ! 書式文字列

    n = size(v)
    write (f, '("(", I0, "(F10.2)", ")")') n
    print f, v(:)
  end subroutine display_vec
end program sle_lu

3. ソースコードのコンパイル

$ gfortran -o sle_lu sle_lu.f95

4. 動作確認

まず、以下のような入力ファイルを用意する。
(先頭行:元の数、2行目〜最終行の前の行:行列 A、最終行:ベクトル b)

File: input.txt

1
2
3
4
5
6
4
3.0  1.0  1.0  1.0
1.0  2.0  1.0 -1.0
1.0  1.0 -4.0  1.0
1.0 -1.0  1.0  1.0
0.0  4.0 -4.0  2.0

そして、実行。

$ ./sle_lu
n = 4
A =
      3.00      1.00      1.00      1.00
      1.00      2.00      1.00     -1.00
      1.00      1.00     -4.00      1.00
      1.00     -1.00      1.00      1.00
b =
      0.00      4.00     -4.00      2.00
x =
     27.00    -28.00    -10.00    -43.00

検算してみると、正しいことが分かる。


以上。

 

Sponsored Link

 

コメントする