連立一次方程式と行列 (その2: 連立一次方程式の解法)

目次

• 分割コンパイル
  • ソースコードの分割
  • ヘッダファイルの利用
  • 課題1 (配点: 2点)
• GNU Make と Makefile
  • 課題2 (配点: 2点)
• ガウスの消去法
  • 後退代入
  • 課題3 (配点: 2点)
• 逆行列の計算
  • 課題4 (配点: 2点)
• ピボット選択
  • 課題5 (配点: 2点)
  • 課題が終わらなかった場合
• 発展: 完全ピボット選択
• 発展: より発展的なMakefileの利用
  • マクロ
  • allターゲット
  • 動的マクロ
  • パターンルール
  • ダミーターゲット
  • さらに実践的な例

今回の講義では，前回の講義で作成した行列タイプを用いて連立一次方程式の解法，ならびに逆行列の計算法について学びます．また，今回の講義で作成する行列タイプは「ライブラリ」として他のプログラムでも利用できるように分割コンパイルとヘッダファイルを利用した再利用の仕組みについても紹介します．

注意

今回からは，プログラムの書き方が汚い（インデントがされていないなど）場合は，採点しません． 前回の述べた見やすいコードの書き方 を参考に，きれいなプログラムを作成してください．

分割コンパイル

今回の本題はガウスの消去法による連立一次方程式の解法ですが，そこに入る前に複数のソースコードからなるプログラムをコンパイルするための分割コンパイルの仕組みについて触れておきます．

前回，行列を扱う関数をいくつか作成しましたが，できれば，これらの関数群を何度でも再利用できるようにしておきたいところです．過去に作成した関数を再利用する場合，初心者は，しばしばコピー＆ペーストに頼りがちです．しかし，コピー＆ペーストで大量にコピーを作っていくと，何かプログラムミスを発見した際の修正が大変です（コピーした数だけ修正しないといけません！）．

分割コンパイルという仕組みを利用すると，作成した関数群を，コピー＆ペーストすることなく簡単に再利用することが可能になります．分割コンパイルの説明をするために，まずは，普通のコンパイルについて復習します．これまで，コンパイルは

gcc source.c

などとしてきましたが，実はこのコマンドを実行すると，内部的には，コンパイル と リンク の二つの作業が行われます．

コンパイルというのは，拡張子が .c のソースファイルを，オブジェクトファイル (通常は拡張子が .o) という中間ファイルに変換する作業であり，リンクというのは，複数のオブジェクトファイルを結合し一つの実行形式ファイル (Windowsの場合であれば拡張子 .exe ですが，CygwinをはじめとするUNIX環境では拡張子無しが一般的です) を作成するための作業です．

オブジェクトファイルを生成する最初のコンパイルの段階では，各関数の中身だけを変換し，関数呼び出しなどの関数間相互の関連や変数参照などは未解決のまま放置されます．そのため，例えば次のプログラムのように，プログラム中で呼び出している関数（func）の実体が存在しないような場合でも，リンクをせずにコンパイルするだけならエラーになりません．

source.c

double func(double x, double y);

int main() {
    double a, b, c;
    a = 3;
    b = 2:

    c = func(a, b);
    return;
}

リンクはせずに，コンパイルだけをするには，

gcc -c source.c

というように -c オプションをつけて gcc を起動すればよく，これにより source.o というオブジェクトファイルが作成されます．また，より明示的に

gcc -c source.c -o source.o

のように出力ファイル名を指定することもできます (出力ファイル名はsource.oでなくてもいい)．

ソースコードの分割

分割コンパイルのおかげで，C言語では，プログラムを複数のファイルに分けて記述することが可能です．たとえば，下記の例のように source1.c と source2.c の二つのファイルに分けてプログラムを記述したとすると，

gcc -c source1.c
gcc -c source2.c

とすることで，それぞれのオブジェクトファイル source1.o と source2.o が作成され，

gcc source1.o source2.o

とすると，両者の間の，関数や変数の外部参照の解決が行われ，リンクされます (いきなり gcc source1.c source2.c としてもOK)．

source1.c

double func(double a, double b);

int main() {
    double a, b, c;
    a = 3.0;
    b = 2.0;
    c = func(a, b);
}

source2.c

double func(double a, double b) {
    return a + b;
}

なお上記の source2.c にはプロトタイプ宣言がありません．もともと，プロトタイプ宣言は main 関数やその他の関数が別の関数を参照する際に，関数名と引数の妥当性をチェックする目的でソースコード上部に書かれるものでした． source2.c においては，関数 func を参照している関数が存在しないので，このようにプロトタイプ宣言を省いて関数を定義することができているのです．

次に，前回の課題3で作成したプログラムを思い出してみましょう．前回は matrix.c に直接 main関数を追加して，以下のようなソースコードを作成していました．

matrix.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

// 要素を交換するマクロ
#define swap(a, b) do { \
    double t = (a);     \
    a = b;              \
    b = t;              \
} while( 0 );

/*
 * 行列用構造体
 * rows: 行数
 * cols: 列数
 * elems: 行列要素を入れた一次元配列
 */
typedef struct {
	int rows;
	int cols;
	double *elems;
} matrix;

// 行列要素を取得するマクロ
#define mat_elem(m, i, j) (m).elems[(i) * (m).cols + (j)]

// ----------------------------------------------------------------------------
// 行列演算用関数群
// ----------------------------------------------------------------------------

// (... 中略 ...)

int main() {
    // 行列の計算
    ...
}

これを，二つのファイルに分割することを考えてみましょう．

分割する場合，分割したファイルの一方を後から再利用できるように，再利用可能な部分と，そうで無い部分の二つに分割します．つまり今回の例では，「行列演算のための関数」を集めたファイル (matrix.c) と，それを利用する main 関数 収めたファイル (main.c) とに分割します．

なお，前回配布した matrix.c は関数のプロトタイプ宣言を省略していましたが，ファイルを分割する場合には main.c をコンパイルするときに関数の情報が必要であるため，以下のようにプロトタイプ宣言を追加します．

main.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

// 要素を交換するマクロ
#define swap(a, b) do { \
    double t = (a);     \
    a = b;              \
    b = t;              \
} while( 0 );

/*
 * 行列用構造体
 * rows: 行数
 * cols: 列数
 * elems: 行列要素を入れた一次元配列
 */
typedef struct {
	int rows;
	int cols;
	double *elems;
} matrix;

// 行列要素を取得するマクロ
#define mat_elem(m, i, j) (m).elems[(i) * (m).cols + (j)]

// ----------------------------------------------------------------------------
// プロトタイプ宣言
// ----------------------------------------------------------------------------

// 行列のメモリ確保 (成否に合わせてブール値を返す)
bool mat_alloc(matrix *mat, int row, int col);
// 行列のメモリ開放
void mat_free(matrix *mat);
// 行列のコピー (サイズが異なるなどコピー不可ならfalseを返す)
bool mat_copy(matrix *dst, matrix src);
// 行列の加算 (サイズが異なるなど計算が不可ならfalseを返す)
bool mat_add(matrix *dst, matrix mat1, matrix mat2);
// 行列の減算 (サイズが異なるなど計算が不可ならfalseを返す)
bool mat_sub(matrix *dst, matrix mat1, matrix mat2);
// 行列の積 (サイズが不適合など計算が不可ならfalseを返す)
bool mat_mul(matrix *dst, matrix mat1, matrix mat2);
// 行列のスカラ倍 (サイズが異なるなど計算が不可ならfalseを返す)
bool mat_muls(matrix *dst, matrix mat, double s);
// 行列の転置 (サイズが不適合ならfalseを返す)
bool mat_trans(matrix *dst, matrix src);
// 行列を単位行列にする (行列が正方でなければfalseを返す)
bool mat_ident(matrix *mat);
// 二つの行列mat1とmat2が等しいかを判定
bool mat_equal(matrix mat1, matrix mat2);
// 行列の中身を標準出力する
void mat_print(matrix mat);

int main() {
    // 行列の計算
    ...
}

一方で，先ほどの source1.c と source2.c の例からも分かる通り matrix.c 側には関数の実態がありますので，こちらはプロトタイプ宣言は不要です．

matrix.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

// 要素を交換するマクロ
#define swap(a, b) do { \
    double t = (a);     \
    a = b;              \
    b = t;              \
} while( 0 );

/*
 * 行列用構造体
 * rows: 行数
 * cols: 列数
 * elems: 行列要素を入れた一次元配列
 */
typedef struct {
	int rows;
	int cols;
	double *elems;
} matrix;

// 行列要素を取得するマクロ
#define mat_elem(m, i, j) (m).elems[(i) * (m).cols + (j)]

// ----------------------------------------------------------------------------
// 行列演算用関数群
// ----------------------------------------------------------------------------

// mat_alloc:行列要素用のメモリを確保する
bool mat_alloc(matrix *mat, int rows, int cols)
{
    ...
}

// mat_free:使い終わった行列のメモリを解放する
void mat_free(matrix *mat)
{
    ...
}

// (... 以下省略 ...)

ここで注意すべき点は，#defineマクロや，matrix構造体を定義している部分は，両方のファイルに共通して書き込む必要があるということです (それぞれのファイルをコンパイルする際に必要とされるため)．

ヘッダファイルの利用

しかし，これらの共通部分を，いちいち複数のファイルに書きこむのは面倒なので，通常は，この共通部分を独立したヘッダファイル (例えば matrix.h) に収めておきmain.c と matrix.c の両者で #include "matrix.h" としてインクルードします．

matrix.h

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

// 要素を交換するマクロ
#define swap(a, b) do { \
    double t = (a);     \
    a = b;              \
    b = t;              \
} while( 0 );

/*
 * 行列用構造体
 * rows: 行数
 * cols: 列数
 * elems: 行列要素を入れた一次元配列
 */
typedef struct {
	int rows;
	int cols;
	double *elems;
} matrix;

// 行列要素を取得するマクロ
#define mat_elem(m, i, j) (m).elems[(i) * (m).cols + (j)]

// ----------------------------------------------------------------------------
// プロトタイプ宣言
// ----------------------------------------------------------------------------

// (... 以下省略 ...)

ヘッダファイルを用いると，以下のように main.c ならびに matrix.c をスッキリとさせることが可能です．

main.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include "matrix.h"

int main() {
    // 中身は省略
    ...
}

matrix.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

#include "matrix.h"

// ----------------------------------------------------------------------------
// 行列演算用関数群
// ----------------------------------------------------------------------------

// mat_alloc:行列要素用のメモリを確保する
bool mat_alloc(matrix *mat, int rows, int cols)
{
    ...
}

// mat_free:使い終わった行列のメモリを解放する
void mat_free(matrix *mat)
{
    ...
}

// (... 以下省略 ...)

なお，この例ではmatrix.h において，プロトタイプ宣言を含めています．これらは (C言語の場合には) 必ずしも matrix.c と main.c の両者に必要とされるわけではありませんが，この例のように matrix.h 中に含めてしまう場合が多いです．

また，標準ヘッダファイルである stdio.h は matrix.h とそれぞれのソースファイルで二重にインクルードされていますが，あくまでそのファイル内で用いる標準ヘッダファイルは別のヘッダファイル (この例では matrix.h) に含めずに，その都度インクルードするほうがバグが少なく，また改変しやすいプログラムとなります．

課題1 (配点: 2点)

前回の課題3で作成したプログラムを，main.c，matrix.c，matrix.h の3つのファイルに分割して，分割コンパイルを実行せよ．

GNU Make と Makefile

分割コンパイルをするときには make コマンド (GNU Make) を使うと便利です．

例えば，上記の例では， main.c と matrix.c という二つのファイルを作りました． すでに一度コンパイル，リンクをした後に matrix.c に間違いを見つけ，matirx.c を修正した場合を考えてみましょう．この場合 main.c には変更が無いので，main.c に関してはコンパイルをする必要がなく matrix.c だけをコンパイルし直して，元からある main.o と，新たにコンパイルした matrix.o とをリンクすればOKです．

このようなときに，変更もしていない main.c を，いちいちコンパイルしなおしていては，時間の無駄です（プログラムのサイズが大きくなってくると，コンパイルに数十秒〜数分かかることもめずらしくないので，この点は重要です）．プログラムの規模が小さい内は，修正したファイルのみを自分で判断して再コンパイルすれば良いですが，プログラムの規模が大きくなり分割ファイル数が増えてくると，このメンテ作業は極めて面倒なものとなります．

make というコマンドを使うと，このようなメンテ作業を自動的に行ってくれます． make は，ファイルのタイムスタンプ（最終更新日の記録）をもとに，ファイルが修正されたかどうかを自動的に判断し，必要なファイルのみを再コンパイルします．

make コマンドを使うには Makefile (先頭のMは大文字．拡張子は無し) という名前のファイルをソースファイルと同じフォルダに作成し，そこにファイルのコンパイル方法，依存関係などを記入していきます．

Makefileの具体例を以下に示します．#で始まる行はコメント行です．この例では program という名前の実行形式ファイルが出力されることを想定しています．

# program という出力ファイルを得るための方法と，必要となる（practice が依存する）ファイル群を定義する．
# 2行目が具体的なコンパイルコマンドを表す．２行目の行頭は，タブ文字である（スペースだとエラーとなる）．
program: main.o matrix.o
	gcc -o program main.o matrix.o

main.o: main.c matrix.h
	gcc -c main.c

matrix.o: matrix.c matrix.h
	gcc -c matrix.c

上記のMakefileでは，まず最初の2行で

「 program は，main.o と matrix.o から生成され，その生成に必要なコマンドは gcc -o program main.o matrix.o です」

というルールを宣言しています．最初の program: ... と書かれた部分でコロンの左側にある program をターゲットと呼びます．一方，右側には matrix.o と main.o が書かれていますが，これらを 生成源 と呼び，もしまだ作成されていなければ Makefile の中に生成ルールが存在するかを探しにいきます．

すると，続く次の2行で main.o をターゲットにした生成ルールが見つかります．このルールは

「main.o は main.c と matrix.h から生成され，生成に必要なコマンドは gcc -c main.c です」

というルールを宣言しており，以下，同様に必要な要素の作成方法を記述しています．なお，ここで 生成源 には matrix.h が指定されているのにGCCには matrix.h が渡されていません．これは matrix.h がプリコンパイラによって暗黙的に main.c の中に展開されるためです．逆に考えると，たとえ main.c が変更されていなくても matrix.h が変更されていれば main.c を再コンパイルする必要があるため，生成源には matrix.h を加えておく必要があります．

GNU Makeはターミナル上で

make

とタイプしてEnterキーを押すと，そのディレクトリにある Makefile を探索して，そこに書かれた一番上のターゲットを実行します．一方で，

make main.o

のようにして main.o のみを最新にアップデートすることもできます．

ところで，生成ルールを示す行 (例えば gcc -c main.c などと書かれた行) の最初に入っている空白は 必ずタブでなければなりません (スペースになっているとエラーになる) ので注意してください．

課題2 (配点: 2点)

課題1で作成したプログラムを分割コンパイルするための Makefile を作成し， make の動作を確認せよ（上に書かれた Makefile をそのまま書き写すだけで良いが，必ず自分で書くこと）．

ガウスの消去法

ここから本題のガウスの消去法について解説します．

今回は，連立一次方程式の解を求める（≒ 逆行列を求める）ための最も一般的な方法である ガウスの消去法 ( 掃き出し法 とも言う) について学びます．ガウスの消去法は，行列の基本変形（通常は行に関する基本変形）を用いて連立方程式の解を求める（あるいは，逆行列を求める）方法です．なお，今回扱うのは変数の数と等式の数が同じであるようなケースのみとします．

\(i\) 行目の等式 (\(i = 2 \ldots n\)) に対して，

の計算をすると，連立方程式は，

となり，2行目以降の1番目の係数が全て消去されます（行変形された式番号や係数に対して，元の値と違うという意味でプライムをつけています）．同様にして，3行目以降の方程式に対し

という計算をすると，3行目以降において2番目の係数も，全て消去されます．同様の操作を順次繰り返していくと，最終的に次のような上三角行列（の形に相当する方程式）が得られます．この上三角行列に対して対角成分全ての積を取ると元々の係数行列の行列式を計算することができます．この行列式が0かどうかを判定することで，行列が正則なのか特異なのか (逆行列を持つか持たないか) を判定することができます．行列が特異な場合には連立方程式は単一の解を持たず不定解となります．

ここまでの操作をガウスの消去法と呼んでいます（広義には，次の後退代入も含めてガウスの消去法と呼ぶ場合もあります）．

後退代入

次に，上記とは反対向きの操作で，対角成分$a_{ii}$ ($i=1, \ldots n$以外の要素を除去していきましょう．
（以下，式が煩雑となるのをさけるため，式番号や係数の ‘ を取り除いて，ガウスの消去法後の式を新たに，$(1)$-$(n)$，また，係数を $a_{ij}$ と表記しなおします．）

まず \(n\) 行目の式をみると \(a_{nn} x_n = b_{n}\) となっているため，これより \(x_n = b_n ~/~ a_{nn}\) という解が得られます．

次に \(n-1\) 行目は \(a_{n-1, n-1} x_{n-1} + a_{n-1, n} x_n = b_{n-1}\) となっているはずなので， \(x_{n-1} = \left( b_{n-1} - a_{n-1, n} x_n \right) ~/~ a_{n}\) とすることで \(x_{n-1}\) の解が求まります．

この処理を一般化すると \(x_i\) の解は，

を，下の行から順に計算していくことで，全ての解が求まることがわかります．

課題3 (配点: 2点)

ガウスの消去法＋後退代入により連立１次方程式を解く関数を作成し，前回のレポート課題のプログラムファイルに付け加えよ． 関数のプロトタイプ宣言は matrix.c に記述がある通り，

bool mat_solve(matrix *x, matrix A, matrix b);

とし \(\mathbf{Ax} = \mathbf{b}\) で与えられる方程式の解を matrix *x に代入することとする．また係数行列の \(\mathbf{a}\) が正則行列でない (行列式が0になる) 場合には方程式が解けないので，そのような場合には false を，正確に方程式が解けた場合には true を返すようにすること．また，関数呼び出しの前後で，引数 A と b の要素の値が書き変わることの無いよう注意せよ．

作成に当たっては以下のスケルトン・プログラムを参考にして良い．

bool mat_solve(matrix *x, matrix A, matrix b) {
    int i, j, k;
    matrix A2, b2;

    // 行列のサイズチェック．この例では，bとxを列ベクトルに限定しているが，
    // この部分は後ほど改良する
    if (A.rows != A.cols) return false;
    if ((b.cols != 1) || (A.cols != b.rows)) return false;
    if ((x->rows != b.rows) || (x->cols != 1)) return false;

    // 行列aと行列bの値を書き換えないよう，別の行列を用意する 
    if (!mat_alloc(&A2, A.rows, A.cols)) return false;
    if (!mat_alloc(&b2, b.rows, b.cols)) return false;

    // 用意した行列にAとbの値をコピーする．以下，これらの行列を用いて計算する
    mat_copy(&A2, a);
    mat_copy(&b2, b);

    /*
     * ガウスの消去法：
     * 普通に作れば10行程度． forループを3つ使う？
     * 行列式がゼロかどうかの判定も忘れないこと
     */

    /*
     *  後退代入：
     *  普通に作れば5-7行程度． forループを2つ使う？
     */

    // 結果を x にコピー
    mat_copy(x, b2);

    // 終わったらメモリを解放
    mat_free(&A2);
    mat_free(&b2);

    return true;
}

プログラムが完成したら前回同様に check_matrix.c を用いてテストを行う．前回は，

gcc -o check_matrix matrix.c check_matrix.c -D_NOT_USE_HEADER

としていたが，今回の講義の最初でヘッダファイルを用意したので，これ以降は

gcc -o check_matrix matrix.c check_matrix.c

とコンパイルすることができる (ヘッダが正しく書けていないとエラーになる)．正しくガウスの消去法が実装できていれば12個あるテストのうち10個が合格するようになる．

逆行列の計算

上記の手法では，連立方程式左辺の係数さえ同じであれば，右辺のみが異なる複数の連立方程式を同時に解くことができます．
例えば，連立方程式を行列の形に書き換えて，

と表される三つの連立方程式があるとすると，これらを一つにまとめて，

ただし \([\mathbf{x}_1 ~~ \mathbf{x}_2 ~~ \mathbf{x}_3 ]\) は，三つの列ベクトルを並べた３列の行列であるとします．

このように書き換えた状態で，この行列で表された方程式に対し，上記と同様にガウスの消去法＋後退代入の基本変形を施していけば，三つの連立方程式を一度に解くことができます．

さて，ある \(n\) 次の正方行列 \(\mathbf{A}\) に対する逆行列というのは，\(\mathbf{A}\mathbf{S} = \mathbf{I}\) (\(\mathbf{I}\) は \(n\) 次正方行列) を満たす \(\mathbf{S}\) のことでした．これは上記の式 \((\alpha)\) と同様の形をしているため，ガウスの消去法と後退代入を用いることで逆行列も求められることが分かります．

課題4 (配点: 2点)

ガウスの消去法＋後退代入により逆行列を計算する関数を作成したい．そこで，課題3で製作した mat_solve 関数を改良して，引数 x, b が列ベクトルで無い場合（すなわち2列以上ある場合）でも，正しく動作するように変更せよ．

作成する関数 mat_inverse のプロトタイプ宣言は

bool mat_inverse(matrix *invA, matrix A);

となっている．この関数の中では内部的に mat_solve を呼び出すことで逆行列を求めるようにプログラムせよ (mat_solveと同じ処理をコピペしたりせず，mat_inverse の中から mat_solve を呼び出すこと)

戻り値については mat_solve と同様に逆行列が得られれば true を，そうでなければ false を返すようにせよ．

なお，関数の作成にあたっては，

mat_inverse(&A, A);

のように同じ行列を入力と出力の両方に指定しても動くように実装し，仮に別々の行列が入出力に与えられ，

mat_inverse(&B, A);

のように呼び出されたときは A の内容が書き変わらないようにせよ．

プログラムが完成したら前回同様に check_matrix.c を用いてテストを行う．正しく書けると12個のテストのうち11個が合格するようになる．

ピボット選択

ガウスの消去法における \(a_{11}\) や \(a'_{22}\), \(a''_{33}\) などのように解を求める過程で除算の分母となっている係数をピボット (pivot) と呼びます．計算の課程においてピボットが 0 となると，0による除算が発生するため，計算を続けることができなくなります．ところが，その場合でも，行基本変形の一つである行の入れ替えを行うと計算が可能となる場合があります．このように，行の入れ替えによって，より性質の良いピボットを選ぶ操作のことを ピボット入れ替え，あるいはピボッティングと呼びます．

具体例を見てみましょう．方程式 \(\mathbf{Ax} = \mathbf{b}\) において，

となる場合を考えます．

まず， 2行目から(1行目×2)を，3行目から(1行目×0.5)，4行目から(1行目×1)をそれぞれ引くと，

となり，次のピボットが 0 となるため，この先の計算の続行が不可能になります．

このような場合は，行の交換を行うことで，計算が可能になります．例えば，2 行目と 4 行目を入れ替えると，

となり，計算を続行できます．

このようにして，ピボットを選択し直す作業を ピボット選択 と呼びます．この操作は連立方程式において，等式の順番を入れ替える操作に対応するので，この操作によって解が変化することはありません．

ピボット選択においては，具体的にどの行と入れ替えを行うかが重要です．一般に，絶対値が小さい値で割り算を行うと，以後の計算において桁落ち誤差を引き起こしやすくなるため，ピボット選択は絶対値が最も大きくなるものを選ぶように行の入れ替え行うと誤差の少ない結果が得られます．

さて，通常，ピボット選択をプログラムに実装する際には，ピボットが 0 となる場合だけではなく，全ての場合に対して，ピボットが最大となるように，ピボット選択を行います．この例をを以下に示します．

まず，ピボット選択を行わない場合には，1 行目の実数倍を，2〜4 行目より引くことで，2〜4 行目の第 1 列を消去しますが，ピボット選択を行う場合には，ピボットの絶対値が最大となるように行入れ替えを行った後に，2 行目以降の第 1 列を消去します．この例では，1 列目を眺めると，2 行目の 4 という数字が最大であることがわかるので，2 行目と 1 行目を入れ替えることで，ピボットの絶対値を最大とすることができます．よって，2 行目と 1 行目を入れ替えて，

として，2 行目以降の第 1 列を消去します．すると，

が得られます．

次に 2 行目を用いて，3 行目以降の第 2 列を消去しますが，その場合にも，ピボット（ここでは2行2列目の値）の絶対値が最大となるように，2 行目と 4 行目を入れ替えて，

としてから，計算を行います (以下，省略)．

課題5 (配点: 2点)

課題4までで作成した mat_solve 関数を，ピボット選択を行うように改良せよ（ピボットが0となる場合以外でも，常に最大のピボットを選ぶようにすること）．

プログラムが完成したら前回同様に check_matrix.c を用いてテストを行う．ここまで完成すると12個全てのテストに合格するようになる．

課題が終わらなかった場合

課題3, 4, 5で作成中の

• main.c
• matrix.h
• matrix.c
• Makefile

の4つのファイルをメールで提出して下さい．加えて，できたところまでのcheck_matrixの実行結果をスクリーンショットとしてメールにて添付してください．

発展: 完全ピボット選択

今回作成したガウスの消去法のプログラムでは連立方程式の行だけを入れ替えていました．本文中にも示した通り，この行の入れ替えは連立方程式に含まれるスカラー係数の一次式の順序を入れ替えるだけなので，連立方程式の解が変化することはありません．

これと同様に， \(x_1, x_2, \ldots, x_n\) とその係数の順序を入れ替える操作は各一次方程式の解を変更しません．当然ながら，

という二式の意味は変わらないわけです．したがって，変数 \(x_1, x_2, \ldots, x_n\)の順序が変化することに注意すると，係数行列の列の順序を入れ替える操作も可能となります．このようにして，最大のピボットを行と列，両方の入れ替えによって選択する操作を完全ピボット選択と言います．

ピボット選択は，特に係数行列に含まれる要素の絶対値の値に大きな大小差があるときに有効なので，完全ピボット選択を用いると，より精度が向上するように思えます．しかしながら，大雑把に考察してみると，完全ピボット選択も部分ピボット選択も最終的に使われるピボットの相対的な大きさにあまり大きな差はなく，実際に精度はそれほど向上しません．

しかも，完全ピボット選択の場合には，列の入れ替えにより解の順序が変わることを，どこかで記録しておかなければならないために，実装が複雑 (実際には少し変えるだけですが)となり，また計算のコストも列入れ替えの分だけ上昇します．そのため，多くの実装においては，計算コストを重視して，部分ピボット選択を用いるのが一般的です．

もし部分ピボット選択で精度が足りないとすれば，LU分解やQR分解を用いたり，ヤコビ法や共役勾配法などの反復法を用いるなど，ガウスの消去法以外の解法を使うのがより有効でしょう．

また，実用的にはピボットの値の絶対値が一定の値 \(\epsilon\) 以下になった時に初めて部分ピボット選択を行い，それでも \(\epsilon\) を上回らなければ完全ピボット選択をするのが効率と精度のバランスに優れてた実装となるようです．

発展: より発展的なMakefileの利用

本文で紹介したMakefileの使い方は本来の使い方からすると，かなり機能を制限しており，一般的な使い方とは言いづらい部分もあります．

ここでは，一般的と言えるレベルでの最低限の使い方を以下のMakefileを元にして解説したいと思います．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

CC		= gcc
SRC		= main.c matrix.c
HDR     = matrix.h
CFLAGS		= -O2 -Wall
LDFLAGS		= -lm
PROGRAM		= main

all: $(PROGRAM)

$(PROGRAM): main.o matrix.o
	$(CC) $(LDFLAGS) $^ -o $@

%.o: %.c
	$(CC) -c $(CFLAGS) $< -o $@

.PHONY: run
run: $(PROGRAM)
	./$(PROGRAM)

.PHONY: clean
clean:
	$(RM) -f *.o $(PROGRAM)

マクロ

まずは上記のMakefileの1行目から5行目を見てください．ここでは，使用するコンパイラの種類や，ソースコードをマクロとして変数に代入しています．この5行はそれぞれ，以下のような意味を表す標準的なマクロ変数名です．

このようにして定義したマクロ変数は $(CC) のような形でターゲットの作成時に利用することができます．

現在 CFLAGS に指定されている2つのオプションのうち -O2 は最適化オプションと呼ばれ，数字の2は0から3の値を取ることができます．-O0 は最適化なしで，数字が大きくなるごとに，よりプログラムが速く動くようになります．一般的には -O2 を利用することが多く，特に速く動かしたいときにだけ -O3 を使います (ただし-O3 は-O2と比べて常に速くなるわけではなく，コンパイル時間は増加します)．2つめの-WallはWarningを全て(all)表示する，というオプションで，バグを生む可能性のあるコードに対して注意を出力してくれます．ときに煩わしいこともありますが，注意の内容からバグが解消することも多いので，特別な理由がなければつけておくことをおすすめします．

次に LDFLAGS についてですが，こちらは外部のライブラリをリンクする際に，ライブラリのある場所やライブラリの名前を指定するのに使います．現在指定されている -lm は特にLinux環境で math.h を使用した場合に必要になります．この講義では，あまりライブラリをリンクする機会は多くありませんが，並列計算の回やグラフィックスプログラミングの回などから使用するようになりますので，頭の片隅にとどめておいてください．

allターゲット

7行目を見てください．Makefileにおいては一番上に all というターゲットを用意して，そのターゲットの生成元として，必要な全てのプログラムを列挙するのが一般的です．

今回は作成するプログラムが一つだけなので，マクロ変数を展開して $(PROGRAM) だけを指定しています．

動的マクロ

10行目を見てください．ここでは，マクロを使ってリンクを行なっていますが，この命令には $^ や $@ といった見慣れない記号が含まれています．これらの記号は動的マクロといって，ターゲット名や生成元を，その直前になされた宣言に基づいて代入します．

動的マクロには多くの種類がありますが，ここではよく使う3種類だけを紹介します．

これに従うと10行目は次のように書いてあるのと同義になります．

gcc main.o matrix.o -o main

今は LDFLAGS が空文字なので少しややこしいですが， $^ は生成元として指定された main.o と matrix.o になり，$@ はターゲット名である main になっています．

パターンルール

12行目を見てください．これはパターンルールといって，生成元のパターンにマッチするファイルを make を実行したディレクトリにおいて探索し，ヒットしたファイルを1つずつ処理していきます．

生成元の部分には %.c と指定されているので，現在のディレクトリにある main.c と matrix.c がそれぞれ，ここに記載された命令で処理されます．

ターゲットの部分にある %.o は生成元である %.c に対応する形で % の部分が置き換えられます．例えば main.c が %.c にマッチしている時には %.o は main.o がターゲットとなります．

ダミーターゲット

19行目を見てください．ここには .PHONY と書かれたあとにコロンで区切って clean が指定されています．また，この clean は直下にターゲットとして定義されています． phony というのは英語で「偽りの」という意味で，この .PHONY の後に書かれたターゲットがダミーターゲットであることを表しています．

ダミーターゲットとは，ターゲットに指定された名前に対応するもの (プログラムはオブジェクトファイルなど) を作成せずに，命令だけを実行するターゲットです．例えば clean ターゲットでは clean という名前のプログラムや clean.o という名前のオブジェクトファイルが作られるわけではありません．

実際には .PHONY の行は省略することができるのですが，たまたま make を実行したディレクトリに clean という名前のファイルがあると .PHONY なしでは clean ターゲットが既に作成済みと判断されて命令が実行されなくなってしまいます．つまり，カレントディクトリの状況に関係なく，必ず命令を実行したい際には，この .PHONY が必要となるわけです．

また，このダミーターゲットを用いると，プログラムを ./a.out のように実行する代わりに make run (ダミーターゲット名は必ずしも run でなくとも良い) というコマンドでプログラムを起動できます．これを行なっているのが 16行目 です．ここでは17行目の命令でプログラムを実行しているのですが，ダミーターゲットの生成元に $(PROGRAM) を指定することによって，もしプログラムがまだ作成されていなかったり，ソースファイルが更新されていたりする場合には，プログラムの作成や更新を行なった上でプログラムを実行してくれます．こうすることで，「コンパイル・リンク」→「実行」というコマンドを make run だけで実行できるようになります．

さらに実践的な例

以上が一般的なMakefileの使い方になります．ここまで覚えると，だいぶプログラム作成やソースコードの修正が楽になると思いますので，ぜひ積極的に活用してください．

ちなみに上記の Makefile でも，まだいくつかの問題があります．例えば，

• ソースファイルが増えた時 $(PROGRAM) の生成源に *.o のファイルを手動で列挙しないといけない
• ヘッダファイル matrix.h が更新されても make がそれを検出してくれない

などの問題が残っています．これらの問題をここで説明するのは本題からそれますので，そこれらの問題を解決した Makefile を掲載するに留めます．詳細は，各自で調べてみてください．

CC		:= gcc
SRC		:= main.c matrix.c
OBJS		:= $(SRC:%.c=%.o)
DEPS		:= $(SRC:%.c=%.d)
HDR		:= matrix.h
CFLAGS		:= -O2 -Wall -MP -MMD
LDFLAGS		:= -lm
PROGRAM		:= main

all: $(PROGRAM)

-include $(DEPS)

$(PROGRAM): $(OBJS)
	$(CC) $(LDFLAGS) $^ -o $@

%.o: %.c
	$(CC) -c $(CFLAGS) $< -o $@

.PHONY: run
run: $(PROGRAM)
	./$(PROGRAM)

.PHONY: clean
clean:
	$(RM) -f $(OBJS) $(DEPS) $(PROGRAM)