【ガウスの消去法(Gaussian elimination method)】Rubyによるガウスの消去法アルゴリズム(構造化プログラミング/ オブジェクト指向プログラミング)

#**はじめに**
先日、[こちらの記事](https://qiita.com/takuyanin/items/3addb4d6e4ff85673e3b)に`Fortran`による**ガウスの消去法アルゴリズム**をまとめました。
今回は、そのプログラムを`ruby`で書き換えたものを掲載するとします。

注：**ガウスの消去法**に関しては上記リンク先に簡単にまとめています。

#**プログラム**
Fortranで記述したプログラムのように**トップダウン式**にプログラムを組む**1. 構造化プログラミング（全体の流れを考えてから、個々の機能を細分化していく方法）**と**ボトムアップ式**にプログラムを組む**2. オブジェクト指向プログラミング（必要な部品をオブジェクトという概念で捉えて、そのオブジェクトをモジュールとして組み合わせて全体を構成する方法）**のそれぞれで組んだプログラムを下記に載せています。


##**1. (構造化プログラミングによる)プログラム**
###**メインプログラム**
メインプログラムには手続き（処理）のみを記述しています。

```main.rb
require 'csv'
require_relative 'gauss_elim'
require_relative 'disp_matrixes'

=begin
変数の意味
@a: 行列A
@b: 行列B
@n_a: 行列Aのサイズ
@n_b_col: 行列Bの列数

変数@a, @bの初期値
@a = [
  [-1, 2, 2],
  [2, 3, -1],
  [1, 1, -1]
]
@b = [
  [1, 2],
  [-3, 2],
  [-2, 2]
]
=end

# 1. 変数の設定
if ARGV[0] && ARGV[1]
  @a = CSV.readlines("./#{ARGV[0]}")
  @b = CSV.readlines("./#{ARGV[1]}")
else
  # 初期値
  @a = CSV.readlines("./a.csv")
  @b = CSV.readlines("./b.csv")
end

@n_a = @a.size
@n_b_col = @b[0].size

# 2. 行列A, Bを表示
display_matrixA_and_matrixB(@a, @b, @n_a, @n_b_col)

# 3. ガウスの消去法を適用
gauss_elim(@a, @b, @n_a, @n_b_col)

# 4. 行列X(解)を表示
display_matrixX(@b, @n_a, @n_b_col)
```

###**サブルーチン1**
サブルーチン1に行列A、B、Xを表示する関数を定義しています。

```disp_matrixes.rb
def display_matrixA_and_matrixB(a, b, n_a, n_b_col)
  # 行列Aを出力
  puts "A ="
  n_a.times do |i|
    n_a.times {|j| printf("%+5f\t", a[i][j])}
    puts ""
  end

  # 行列Bを出力
  puts "B ="
  n_a.times do |i|
    n_b_col.times {|j| printf("%+5f\t", b[i][j])}
    puts ""
  end
end

def display_matrixX(b, n_a, n_b_col)
  # 行列X(解)を出力
  puts "X (satisfying AX=B) ="
  n_a.times do |i|
    n_b_col.times {|j| printf("%+5f\t", b[i][j])}
    puts ""
  end
end
```

###**サブルーチン2**
サブルーチン2にガウスの消去法関数を定義しています。

```gauss_elim.rb
require_relative 'ppivot'

# ガウスの消去法
def gauss_elim(a, b, n_a, n_b_col)
  # 前進消去
  (n_a - 1).times do |j|
    partial_pivoting(a, b, n_a, n_b_col, j)
    # ピボット列の掃き出し
    (j + 1).upto(n_a - 1) do |i|
      r = a[i][j].to_f / a[j][j].to_f
      # 各行減算
      # 1.行列Aについて
      j.upto(n_a - 1) do |k|
        a[i][k] = a[i][k].to_f - a[j][k].to_f * r
      end
      # 2.行列Bについて
      0.upto(n_b_col - 1) do |k|
        b[i][k] = b[i][k].to_f - b[j][k].to_f * r
      end
    end
  end

  # 後退代入
  (n_b_col).times do |j|
    (n_a - 1).downto(0) do |i|
      r = b[i][j].to_f
      # シグマ演算
      (i + 1).upto(n_a - 1) do |k|
        # ここでのrはx
        r = r - a[i][k].to_f * b[k][j].to_f
      end
      b[i][j] = r / a[i][i].to_f
    end
  end
end
```

###**サブルーチン3**
サブルーチン3に部分ピボット選択関数を定義しています。

```ppivot.rb
# 部分ピボット選択
def partial_pivoting(a, b, n_a, n_b_col, j)
  #最大値を記録
  amax = a[j][j].to_f.abs
  imax = j

  (j + 1).upto(n_a - 1) do |i|
    r = a[i][j].to_f.abs
    if r > amax
    # 最大値を変更
      amax = r
      imax = i
    end
  end

  return if j == imax
  # ピボット行と最大値行が異なれば交換
  # 1.行列Aについて
  n_a.times do |k|
    r = a[imax][k]
    a[imax][k] = a[j][k]
    a[j][k] = r
  end
  # 2.行列Bについて
  n_b_col.times do |k|
    r = b[imax][k]
    b[imax][k] = b[j][k]
    b[j][k] = r
  end
end
```

###**実行**
```
$ ruby main.rb
```
行列を定義したCSVファイル(例： a.csv, b.csv)を指定することもできます。

```
$ ruby main.rb a.csv b.csv
```

一応、CSVファイルの例としてa.csvとb.csvの中身を下記に載せておきます。

###**CSVファイル**
```a.csv
-1, 2, 2
2, 3, -1
1, 1, -1
```

```b.csv
1, 2
-3, 2
-2, 2
```

###**結果**
```
$ ruby main.rb a.csv b.csv
A =
-1.000000	+2.000000	+2.000000
+2.000000	+3.000000	-1.000000
+1.000000	+1.000000	-1.000000
B =
+1.000000	+2.000000
-3.000000	+2.000000
-2.000000	+2.000000
X (satisfying AX=B) =
+1.000000	-6.000000
-1.000000	+3.000000
+2.000000	-5.000000
```

##**2. （オブジェクト指向プログラミングによる）プログラム**

###**メインプログラム**
計算を行うメインプログラムファイルを準備します。
このプログラムに必要モジュールをインポートします。

```calc.rb
require './gauss_elimination.rb'

if __FILE__ == $0
  begin
    # 計算クラスインスタンス化
    obj = GaussElimination.new
    # ガウスの消去法で連立方程式を解く
    obj.display_matrixA_and_matrixB
    obj.calculate
    obj.display_matrixX
  rescue => e
    $stderr.puts "[#{e.class}]"
    $stderr.puts "#{e.message}"
    e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
  end
end
```

###**モジュール1**
**ガウスの消去法クラス**を設定し、その中に必要メソッドを定義します。

注：今回はモジュールは一つのみですが、複数モジュールを準備できるということを示すために、便宜上モジュール1としています。

```gauss_elimination.rb
require 'csv'

class GaussElimination
=begin
  A. インスタンス変数
    @a: 行列A
    @b: 行列B
    @n_a: 行列Aの次元
    @n_b_col: 行列Bの列数

  B. メソッド
    initialize: コンストラクタ
    display_matrixA_and_matrixB: 行列A, Bを出力
    display_matrixX: 行列X(解)を出力
    calculate: ガウスの消去法で計算を行う
    partial_pivoting: 部分ピボット選択
=end

  def initialize
    # 初期値
    @a = [
      [-1, 2, 2],
      [2, 3, -1],
      [1, 1, -1]
    ]
    @b = [
      [1, 2],
      [-3, 2],
      [-2, 2]
    ]

    if ARGV[0] && ARGV[1]
      @a = CSV.readlines("./#{ARGV[0]}")
      @b = CSV.readlines("./#{ARGV[1]}")
    end

    @n_a = @a.size
    @n_b_col = @b[0].size
  end

  def display_matrixA_and_matrixB
    # 行列Aを出力
    puts "A ="
    @n_a.times do |i|
      @n_a.times {|j| printf("%+.5f\t", @a[i][j])}
      puts ""
    end

    # 行列Bを出力
    puts "B ="
    @n_a.times do |i|
      @n_b_col.times {|j| printf("%+.5f\t", @b[i][j])}
      puts ""
    end
  end

  def display_matrixX
    # 行列X(解)を出力
    puts "X (satisfying AX=B) ="
    @n_a.times do |i|
      @n_b_col.times {|j| printf("%+.5f\t", @b[i][j])}
      puts ""
    end
  end

  def calculate
    # 前進消去
    (@n_a - 1).times do |j|
      partial_pivoting(j)
      # ピボット列の掃き出し
      (j + 1).upto(@n_a - 1) do |i|
        r = @a[i][j].to_f / @a[j][j].to_f
        # 各行減算
        # 1.行列Aについて
        j.upto(@n_a - 1) do |k|
          @a[i][k] = @a[i][k].to_f - @a[j][k].to_f * r
        end
        # 2.行列Bについて
        0.upto(@n_b_col - 1) do |k|
          @b[i][k] = @b[i][k].to_f - @b[j][k].to_f * r
        end
      end
    end

    # 後退代入
    (@n_b_col).times do |j|
      (@n_a - 1).downto(0) do |i|
        r = @b[i][j].to_f
        # シグマ演算
        (i + 1).upto(@n_a - 1) do |k|
          # ここでのrはx
          r = r - @a[i][k].to_f * @b[k][j].to_f
        end
        @b[i][j] = r / @a[i][i].to_f
      end
    end
  end

  private
  def partial_pivoting(j)
    print @a
    puts ""
    #最大値を記録
    amax = @a[j][j].to_f.abs
    imax = j

    (j + 1).upto(@n_a - 1) do |i|
      r = @a[i][j].to_f.abs
      if r > amax
        # 最大値を変更
        amax = r
        imax = i
      end
    end

    return if j == imax
    # ピボット行と最大値行が異なれば交換
    # 1.行列Aについて
    @n_a.times do |k|
      r = @a[imax][k]
      @a[imax][k] = @a[j][k]
      @a[j][k] = r
    end
    # 2.行列Bについて
    @n_b_col.times do |k|
      r = @b[imax][k]
      @b[imax][k] = @b[j][k]
      @b[j][k] = r
    end
  end
end
```

###**実行**
```
$ ruby calc.rb
```
行列を定義したCSVファイル(例： c.csv, d.csv)を指定することもできます。

```
$ ruby calc.rb c.csv d.csv
```

一応、CSVファイルの例としてc.csvとd.csvの中身を下記に載せておきます。

###**CSVファイル**
```c.csv
2, 4, 2
2, 0, 1
1, 1, 1
```

```d.csv
6, 1, 0, 0
2, 0, 1, 0
1, 0, 0, 1
```

###**結果**
```
$ ruby calc.rb c.csv d.csv
A =
+2.000000	+4.000000	+2.000000
+2.000000	+0.000000	+1.000000
+1.000000	+1.000000	+1.000000
B =
+6.000000	+1.000000	+0.000000	+0.000000
+2.000000	+0.000000	+1.000000	+0.000000
+1.000000	+0.000000	+0.000000	+1.000000
X (satisfying AX=B) =
+3.000000	+0.500000	+1.000000	-2.000000
+2.000000	+0.500000	-0.000000	-1.000000
-4.000000	-1.000000	-1.000000	+4.000000
```

#**おわりに**
今回、構造化プログラミンング、オブジェクト指向プログラミングそれぞれでガウスの消去法プラグラムを組んでみましたが、**臨機応変にどちらの手法を採用すべきか**考えていきたいものです。
いずれにしても、メソッド（関数）、クラスの中身はブラックボックスとして捉えることができるので、メインプログラムを見たときに、**何をしているのか直感的に理解できるように記述することが大切**です。

こちらの記事が役にたったという方は、いいね、よろしくお願いします(^^)