連立一次方程式と行列 (その1: 行列計算用変数タイプの設計)

目次

今回は行列計算のためのライブラリーの作成が主になります．前回の授業で，ポインタ，構造体は完全に理解できていないと，作成は大変困難を極めます．あまり自身がない人は再度ポインタの使い方を読み返して，その役割を理解しておきましょう．

連立一次方程式

数値解析における問題の多くが，最終的に連立一次方程式を解くことに帰着します．そのため，連立方程式の解法は，数値解析の中でも重要な要素の一つです． 連立一次方程式を計算機で扱う場合，方程式の要素を行列を用いて表します．例えば，

\[\begin{cases} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n = b_1 \\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n = b_2 \\ \qquad \vdots \\ a_{m1} x_1 + a_{m2} x_2 + \cdots + a_{mn} x_n = b_m \\ \end{cases}\]

で表される \(m\) 元連立方程式は，

\[\mathbf{A} = \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \\ \end{pmatrix}, \quad \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}, \quad \mathbf{b} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{pmatrix}\]

とおくことにより，

\[\mathbf{A x} = \mathbf{b}\]

と，行列の形で表すことができます．

行列用変数タイプの設計

行列には，行と列の2次元の要素があるため，行列は本来2次元の配列です．C言語では，2次元配列を扱うこともできますが，2次元配列の動的確保は面倒なことが多いので，ここでは，行列を1次元配列を使って表現します．
(2次元配列を用いた定義に興味がある人はこちらを参考にしてみてください）

配列の添え字について

C言語における配列の添字は，”0”から始まります．一方，数学における行列の列数・行数は，通常 “1”から数えます． よって，行列をC言語で扱う場合，添字（or 行数・列数）を0から始めるか，1から始めるかをまず決めないといけません．

これに関してはいろいろな考え方がありますが，現代の比較的新しいプログラミング言語においては，配列の添え字は0から始まるのが一般的であるため，それに慣れる意味も含め，本講義では0から始まる添え字を使ってプログラムを作成していきます．一方で数学的な説明については，前回から引き続き1から始まる添え字を用います．

このような書き換えは，実際に自分で本を読みながらプログラミングをする際にも役立つものですので，ぜひこの表現の違いに慣れてください．

補足として，一部のプログラミング言語 (FortranやMATLABなど) は配列が1から始まる添え字を採用しています．

1次元配列による行列要素の表現

\(m\) 行 \(\times\) \(n\) 列 の行列を1次元の配列で表現するには，次のようにします．

  1. 行列の要素数は全部で \(m \times n\) 個になるので， \(m \times n\) 個の要素を持った1次元配列を作る．
  2. 行列の \(i\) 行 \(j\) 列の要素を，1次元配列の \(n \times i + j\) 番目の要素に割り当てる．

行列が持つ情報には，各要素の値のほかに，行数，列数の情報があります．行数，列数がわからないと，各要素にアクセスすることができないので，これらの情報を別途，何らかの変数で管理する必要があります．よって行列を効率良く扱うためには，「要素」「行数」「列数」をまとめて管理できる変数の型があると大変便利です．C言語では，新しい変数型を自分で設計できますので，「要素」「行数」「列数」の３つの情報を持つ変数型を設計してみましょう．

本日の課題の進め方

本日の課題を作成するにあたって，以下のZIPファイルをダウンロードしてください．

このZIPファイルの中には matrix.c と check_matrix.c という2つのファイルが含まれています．まずはZIPファイルを展開し，得られる matrix_template フォルダを matrix に名前変更して，自分の作業ディレクトリに配置しましょう (例えば ~/prog/matrix)．

配置ができたら，ターミナルからフォルダの中に移動し，以下のコマンドでテストプログラムをコンパイルしましょう．

gcc -o check_matrix check_matrix.c -D_NOT_USE_HEADER

こちらのコマンドでは _NOT_USE_HEADERというマクロによって matrix.c の中身を check_matrix.c の中に含めるようにしてあるので，gccに matrix.c を与える必要がないことに注意してください (なお，次回の課題では必要になります)．

すると check_matrix という名前で実行ファイルが作成されるので，これをシェルから

./check_matrix

として実行します．

するとターミナル上に以下のようなステータスがいくつか表示されます．

[ RUN ] mat_alloc_and_free
[ ERROR ] expected = true, actual = false (Line 116)
[ FAILED ] mat_alloc_and_free
=====

(...中略...)

=====
[ RUN ] mat_inverse
[ ERROR ] expected = true, actual = false (Line 532)
[ FAILED ] mat_inverse
=====
12 / 12 tests failed,
[ STATUS ] FAILED

このステータスは，どの関数が正しく実装されているかを確認していて，例えば mat_alloc と mat_free が正しく実装されると [ ERROR ] と書かれた行が消えて [ FAILED ] という表示が [ PASSED ] に変わります．その上で，一番下に表示される失敗テストの個数が1つ減ります．

[ RUN ] mat_alloc_and_free
[ PASSED ] mat_alloc_and_free
====

(...中略...)

=====
[ RUN ] mat_inverse
[ ERROR ] expected = true, actual = false (Line 532)
[ FAILED ] mat_inverse
=====
11 / 12 tests failed,
[ STATUS ] FAILED

今回と次回の講義では，1つずつ matrix.c 内の関数を実装していき，最終的に全てのテストが通るようにします．このとき [ ERROR ]と書かれた行の一番最後にある (Line 116) というのは check_matrix.c 内の行番号です．ここを見てみると，

ASSERT_TRUE(mat_alloc(&A, 4, 5));

のように書かれています．これは mat_alloc の戻り値が true であることを assert (断言する，つまり保証するということ) するという意味です．

ここでテストが失敗しているということは mat_alloc が true を返していないということです．もちろん，単に true を返すようにするだけでは他の部分でテストが失敗してしまいますので，まずは関数を自分なりに実装してみて，どのくらいのテストに合格するかを見てみるのが良いと思います．

課題1 (配点: 2点)

行列用変数タイプ matrix を設計せよ．すでに matrix.c には matrix 型の構造体の定義と，行列の要素にアクセスするためのマクロ mat_elem が与えられている (行列の要素番号が 0 から始まることに注意)．

なお，以下のプログラムでは，各処理が成功したときにはブール値のtrueを，そうでないときにはfalseを返すようにしている．このようなブール値はC99から導入されたもので，使用するには stdbool.h をインクルードする必要があることに注意 (C++では標準型になったので必要ない)．

このテンプレートファイルを元にして，行列をメモリ上に確保する関数(mat_alloc)と，それを解放する関数(mat_free)を設計せよ．設計がうまくいくと12個あるテストのうち1つが合格するようになる．

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) {
    return false;
}

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

課題2 (配点: 6点)

matrix.c に与えられた以下のようなプロトタイプ宣言に従い，行列を扱う関数群をプログラムせよ (引数の順序に注意)．これらが全て完了すると，12個あるテストのうち9個が合格するようになる．

注意点

mat_mul と mat_trans に関しては mat_mul(&a, a, b) というように，結果と被演算変数に同じ変数を指定してもうまく動くように作業用の行列を関数内で用意すること．

matrix 型は内部にポインタを使っているため，意図せずに要素を書き換えてしまう場合がある．特に，複数のmatrix型変数のelemsが同じメモリ領域を指している場合は注意が必要である．どのメモリ領域を書き換えているのかを意識しながらプログラムを書くことが重要．

matrix.c (一部)

// mat_copy: srcの中身を*dstにコピーする
bool mat_copy(matrix *dst, matrix src) {
    return false;
}

// mat_add: mat2+mat3を*mat1に代入する
bool mat_add(matrix *res, matrix mat1, matrix mat2) {
    return false;
}

// mat_sub: mat2-mat3を*mat1に代入する
bool mat_sub(matrix *res, matrix mat1, matrix mat2) {
    return false;
}

// mat_mul: mat2とmat3の行列積を*mat1に代入する
bool mat_mul(matrix *res, matrix mat1, matrix mat2) {
    return false;
}


// mat_muls: mat2をc倍（スカラー倍）した結果を*mat1に代入する
bool mat_muls(matrix *res, matrix mat, double c) {
    return false;
}

// mat_trans: mat2の転置行列を*mat1に代入する
bool mat_trans(matrix *res, matrix mat) {
    return false;
}

// mat_ident: 単位行列を与える
bool mat_ident(matrix *mat) {
    return false;
}

// mat_equal: mat1とmat2が等しければtrueを返す
bool mat_equal(matrix mat1, matrix mat2) {
    return false;
}

課題3 (配点: 2点)

上記のテストに合格(12個中9個だけ合格すればよい)したら，実際にmatrix型を使ってプログラムを書いてみましょう．現在 matrix.cをコピーしたファイルmatrix2.cを作成し，そのファイルにmain関数を追加して，以下の計算を行ってください．なお matrix.cをコピーせずに直接main関数を付け加えると，上記の check_matrix.c のコンパイルが通らなくなるため，必ず，matrix2.cのような別のソースコードのファイルを作成して，課題に取り組んでください．

(1)     \(\mathbf{A} = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}\) に対して \(3 \mathbf{A}^T\)

(2)     \(\mathbf{A} = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} -1 & 2 & -1 \\ 2 & -1 & 2 \\ -1 & 2 & -1 \end{pmatrix}, \quad \mathbf{C} = \begin{pmatrix} 3 & 1 & 4 \\ 1 & 5 & 9 \\ 2 & 6 & 5 \end{pmatrix}\) に対して \(\mathbf{C}^T (5 \mathbf{A} - \mathbf{B} / 2)\)

なお，それぞれの結果は以下のようになる．

(1)     \(\begin{pmatrix} 3 & 12 & 21 \\ 6 & 15 & 24 \\ 9 & 18 & 27 \\ \end{pmatrix}\)

(2)     \(\begin{pmatrix} 106.5 & 130.5 & 166.5 \\ 313.5 & 370.5 & 433.5 \\ 370.5 & 460.5 & 550.5 \\ \end{pmatrix}\)