計算機言語 II 第 2 回 線形代数とプログラム
http://www.math.u-ryukyu.ac.jp/~suga/gengo/2019-2/02.pdf
1 線形代数をプログラムする
1年次の授業で線形代数学を学習しましたが,その内容は,とても広く応用されています. 特に,「連立1次 方程式を解く」は, 応用科学の分野では重要で,もちろんコンピュータを用いてなされます. さらに,「線形変 換」は,シミュレーションやコンピュータグラフィックスの分野では重要で, これが行列の積で記述できると いうのは,実は線型代数学の最も重要な結果といっても過言ではありません. これ以外にも固有値の計算や,逆 行列の計算なども重要です.
連立方程式の解を求めたり,逆行列を求める計算は, 線形代数で学習した「行列の行基本変形を用いた掃き 出し」が, 1つの基本的な計算手法です(固有値, 固有ベクトルの計算は大変.). これらの計算をコンピュータ に効率よく計算させる手法を述べます.
注意 「正則行列は, 行基本変形で逆行列が求められる」というのは, 線型代数学の重要な結果ですが,
「そんなことをしなくても逆行列がわかる」こともあります. 例えば, グラフィックスのプログラムで,
「原点を中心に回転する」という操作をプログラムする際には,「直交行列による線形変換」をプログラ ムしますが,直交行列の逆行列は転置行列なので,行列要素の入れ替えだけで,逆行列が求まります. 前 期の授業でも述べましたが,何事も「臨機応変」な態度が重要ですし, このような基本的な事を理解し ている事も重要です.
行基本変形は,次の3つの操作からなります. 1. ある行を(0でない)定数倍する.
2. 2つの行を入れ替える.
3. ある行に別の行の定数倍を加える.
この中で, 2番目はうまくプログラムを組むことで,動作の高速化ができます. 行列を単純な2次元配列とす るのではなく,行ベクトルへのポインタの配列にするのです. すなわち, 列ベクトルの大きさ(行ベクトルの個 数 を持つポインタ配列 を準備し
という風にします. こうすると,i行とj行の交換は,m[i] とm[j]の交換で済みます. (素朴に行ベクトルの 値の交換は, 行ベクトルの個数だけ実行しなければならない.) また, このようにデータを用意しておくと,i, j 成分はm[i][j]でアクセスでき, 行列計算(和や積)のプログラムも容易です.
2 動的なメモリの確保(malloc())
行列計算や連立方程式を解く等のプログラムを少し汎用的に書こうとすると,行列の大きさや連立方程式の 本数も変化する量としてプログラムを書く必要が出てきます.
上のように配列を利用するのですが, 可変長配列というのがC11(2011年制定の規格)で使えるようになり ました. しかし,それは局所的な使い方(ある関数の内部だけでしか使えない)しかできまので, これから書く ようなプログラムには向きません.
前期の授業で少し述べましたが, Cではmalloc()というライブラリ関数があり, これはシステムに引数で 与えらえる量のメモリの確保を要求するものです. 返り値は,以前解説したvoid型のポインタで,これを必要 な型にキャストして使います. malloc()で確保したメモリは, 不要になればfree()という関数で, 再び確保 可能な領域に戻すことができます. ファイルの読み書きに対するfopen()とfclose()と似たような関係で す. ファイルのときと同様に,プログラム終了時には,malloc()で確保したメモリは自動的にfree()されま す. 従って,すぐに終了するようなプログラムではfree()はあまり意味を持ちませんが,巨大なデータを扱っ たり(最近流行りのbig data),ずっと動き続けるようなプログラムでは,free()を忘れてはいけません.
malloc()は, メモリ確保に失敗すると, NULL ポインタを返します. 以前に述べましたが, NULLポイン タは,「どこも指していない」ことが明確になっているポインタです.
これから,何回かに渡って行列計算のプログラムを書く予定です. このような際には, ディレクトリ(フォル ダ)を作って, そこで作業するのが通常です. 適当な場所に, 例えばmatrixというディレクトリを作り,そこ
にworkingディレクトリを移動します.
Bash3-2$ mkdir matrix Bash3-2$ cd matrix
上で述べたことを踏まえて, 下の場所にあるプログラムの解説をします. (指定した教科書のプログラムは, 添字をずらすなどの操作が入っており,読みづらい)
https://github.com/okumuralab/algo-c/blob/master/src/matutil.c
このソースを上で作った matrixというフォルダにダウンロードして(あるいは, Gedit を用いてプログラ ム内容をコピペして),コンパイル,実行してみてください. main()の前後にある#if 0 と#endifを取り除 くかコメントアウトすると,実行可能ファイルが作れるソースになります.
3 分割コンパイル
Cではいくつかに分割したソースコードを使って, ひとつのプログラムを作成することができます. コンパ イルの方法は単純で,file1.c, file2.c, ...から1つのプログラムを作るには,
cc file1.c file2.c ...
と順番にファイル名を並べるだけです.
この分割コンパイルの仕組みがあったため, C で多くのプログラムが書かれました. 大きなプログラム(例
えば, Linuxのカーネル)では, 多くの人が機能単位でプログラムを別ファイルに記述して行き, それを最後に
まとめあげて最終成果物を作り上げます.
これまで利用してきたライブラリ関数は,多くの場合, 1つのファイルに1つの関数だけが書かれていて,そ れをまとめたものです.
#include <stdio.h>のようなヘッダファイルの読み込みは,分割コンパイルの別の仕組みです. これは, コンパイル時に別のファイルを取り込むという操作をしています. 上のコマンドで複数ファイルを指定した場 合は,それぞれのファイルをコンパイルして機械語ファイルを作った後, 最後にそれをまとめるという動作を gccコマンドが面倒を見ます. (最後のまとめる部分は, 実際にはld(リンクエディタとかリンカと呼ばれる) というコマンドが実行しますが,それは gccコマンドが面倒を見ます.)
これらの仕組みの使い分けですが,基本的に, ヘッダファイルにはcppが扱うマクロ置換や,様々なコンパ イル時の定義および関数のプロトタイプ宣言を書き, Cのソースにはプログラム本体を書きます.
4 matutil.c をライブラリ的に使う
先程のmatutil.cのmain()以外の関数は, 線形代数(行列)をプログラムする上で便利な関数が多く記述さ
れています. これをライブラリ的に使うようにします.
このフォルダに, matutil.hというヘッダファイルと,matutil.cというCのソースを作ります.
matutil.hは, matutil.cで利用するマクロ置換と記述する関数のプロトタイプ宣言を書きます. 次のよ うな内容になります.https://github.com/okumuralab/algo-c/blob/master/src/matutil.cをダウン ロードするかコピペして適当に編集するのが早道です.
/***********************************************************
matutil.h: Header file for matutil.c
***********************************************************/
#ifndef MATUTIL
#define MATUTIL
#endif
#ifndef SCALAR
#define SCALAR double
#endif
typedef SCALAR *vector, **matrix;
/* ここからが関数のプロトタイプ宣言 */
void error(char *message);
vector newvec(int n);
matrix newmat(int nrow, int ncol);
vector new_vector(int n);
matrix new_matrix(int nrow, int ncol);
void free_vector(vector v);
void free_matrix(matrix a);
double innerproduct(int n, vector u, vector v);
void vecprint(vector v, int n, int perline, char *format);
void matprint(matrix a, int ncol, int perline, char *format);
/* とりあえずここまで */
次にmatutil.cを書きます. 上のヘッダファイルがあるので, matutil.cは次ページのようになります. やはり,上でダウンロードしたファイルを編集(一部の追加と不要な部分の削除)します.
/***********************************************************
matutil.c -- 行列
***********************************************************/
/* 行列操作の小道具集 */
#include <stdio.h>
#include <stdlib.h>
#include "matutil.h" /* この行が付け加わる */
void error(char *message) {
fprintf(stderr, "\n%s\n", message); exit(1);
}
vector newvec(int n) {
return malloc(sizeof(SCALAR) * n);
}
... 中略 ...
長いので途中は省略するが, 最初の関数 error から下の関数 matprint まではそのまま利用して, それ以降は全部削除する.
... 中略 ...
void matprint(matrix a, int ncol, int perline, char *format) {
int i;
for (i = 0; a[i] != NULL; i++) {
vecprint(a[i], ncol, perline, format);
if (ncol > perline){
printf("\n");
} }
このようにしておくと,matutil.cで記述された関数を利用するソースを書くには,ソースコードの初めの 方にに#include "matutil.h"で,matutil.hを読み込んでおくだけです. #includeでファイルを読み込む には,ファイルまでのパス名を指定します. < >の中にあるファイルは,コンパイラ処理系が標準的にに指定し たところにあるファイル, Linuxだと/usr/includeにあるファイルを読み込みます.
matutil.oにある関数を用いたプログラムをコンパイルには, 次のようにタイプすれば大丈夫です.
Bash3-2$ gcc hohohoge.c matutil.o
分割コンパイルを実行してみる
最後に,分割コンパイルを試してみます. matutil.cのmain部分をコピーして,hogehoge.c を作ります. ファイルの先頭部分にのヘッダファイルの取り込みは,次のようにします.
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "matutil.h"
上のコマンドを実行して,実行可能ファイルができることを確かめて下さい.
レポート問題(締め切り 10月23日(水))
行列を読み込んで, その転置行列を出力するプログラムを書いて下さい. 上で作った, matutil.oは必ず利 用するようにし,次の形のプログラムにして下さい. (メールの件名は, Enshu 2.1)
• 入力行列と転置行列の両方を保存するメモリを確保し, 転置行列には入力行列の値からコピーをして転 置行列を出力する. 入力は,行列のサイズ,各成分を入力するものとする(下を参考にせよ).
行の数を入力してください -> 1 列の数を入力してください -> 2 0,0 成分を入力して下さい -> 3 0,1 成分を入力して下さい -> 4 転置行列->
3.00000 4.00000
このPDFソースの置き場所: ftp://ftp.math.u-ryukyu.ac.jp/pub/gengo/2019-2/
