\documentclass[a4paper,11pt]{jsarticle}
\usepackage[dvips]{graphicx}
\newcommand{\mathbd}[1]{\mbox{\boldmath $#1$}}
\title{計算機工学概論課題}
\author{計数工学科 4年 30500 岡山友昭}
\begin{document}
\maketitle
\section{簡易ベクトル・行列演算ライブラリ}
講義で配られたCプログラム、 matrixvector.c の骨組みを元に、改善や
新たな実装を行った。主な部分を以下に示す。

\subsection{行列の記憶領域確保と解放(AllocMat, FreeMat)}
行列の記憶領域確保(AllocMat)は、例えば A が m×n 行列の場合には
A[m] に1次元配列として m $\ast$ n 個分のデータを格納し、
A[i] (i = 0, 1, ... m-1) は A[m][i $\ast$ n] を指すポインタとなっている。
例えば 4×4 行列の場合は図\ref{4x4} のようになっている。

\begin{figure}[h]
\begin{center}
\includegraphics[clip]{matrix.eps}
\end{center}
\caption{4×4行列のデータ格納の様子}
\label{4x4}
\end{figure}

初めは A[0] にデータを格納していたのとは違い、 A[m] にデータを格納する
ようにしたため、この記憶領域を解放する FreeMat 関数には列のサイズを
渡さなければならなくなったことに注意。

\subsection{行列の行の入れ換え(SwapRowsInMat)}
ガウスの消去法や LU 分解の際に必須な行の入れ換えをする関数、
SwapRowsInMat は、初めは律義に全ての要素をとりかえていたが、せっかく
データの格納を先ほどのような形にしたので、その利点を生かしてポインタの
入れ換えをするだけにした。これでここの操作のオーダーが m から 1 になった。
これとは違い、列の入れ換えをする関数 SwapColsInMat は、まともに全ての
要素をとりかえる操作をする(遅い)。

\subsection{行列の無限大ノルム(InfNormMat)}
行列の無限大ノルムを求める InfNormMat を実装した。ここで、
$n \times n$ の実数を要素とする行列 $A = \{ a_{ij} \}$ に対する
最大値作用素ノルム $\Vert A \Vert_{\infty}$ は
$$
\Vert A \Vert_{\infty} :=
 \max_{\Vert \mathbd{x} \Vert_{\infty} \leq 1}\Vert A\mathbd{x}\Vert_{\infty}
 =
 \max_{1\leq i\leq n}\sum_{j=1}^{n}|a_{ij}|
$$
で与えられることを用いている。

\subsection{LU分解(LUdecomp)}
逆行列を求めるときや連立1次方程式を解く際にとても有効な手法、LU分解の
メソッド LUdecomp を実装した。この際に参考にしたサイトを最後(参考文献)に
挙げておく。これらのプログラムに比べて私の実装したものの良い点を以下に示す。
\begin{enumerate}
 \item ピボット選択を行っていること(というかメソッド中ほとんどがこのコード)
 \item ピボット選択の際、行のスケーリング(重みづけ)を[きちんと]していること
 \item 行列式も求められること。ただし計算量は全く変わらない。
 \item 行の入れ換え操作はオーダー1でできること(これは上で述べた)
 \item 与えられた行列に対して破壊的な操作を行わないこと
 \item LU分解した情報を保持できる形で返していること
\end{enumerate}

5. は、与えられた行列をコピーした新しい行列を用意してそれに操作すれば
いいものを、与えられた行列そのものに対して操作をしてしまうという
プログラムがあったことを言っている。

6. については、LU分解した情報を何度も使いまわせるということがLU分解の
利点であるはずなのに、その情報を捨ててしまうものが多かったということで
ある。逆行列はそれでもガウス・ジョルダンの方法と計算量は変わらないが、
連立1次方程式を解く際にはLU分解をしなければいけない分、計算量がガウスの
消去法よりも増えてしまう。1回連立1次方程式を解くだけなら
ガウスの消去法で解く SolvLinEqGauss メソッドを使えば十分である。
同じ行列を使って何回も問題を解く場合は、LUdecomp と以下の SolvLinEqLU を
使うとよい。

\subsection{LU分解を利用して逆行列を求める(InvMatLU)}
LU分解した情報を利用して、逆行列を求めるメソッド InvMatLU を
実装した。これを行う前に必ず LUdecomp を実行して、LU分解済みの
行列とピボット選択情報を InvMatLU に渡す必要がある。

\subsection{LU分解を利用して連立1次方程式を解く(SolvLinEqLU)}
LU分解した情報を利用して、連立1次方程式を解くメソッド SolvLinEqLU を
実装した。これを行う前に必ず LUdecomp を実行して、LU分解済みの
行列とピボット選択情報を SolvLinEqLU に渡す必要がある。

\section{連立1次方程式の精度保証つき計算}
実装した簡易ベクトル・行列演算ライブラリ、matrixvector.c を使って、
連立1次方程式を精度保証つきで計算するプログラム accuracy.c を実装した。
理論については適当な文献参照のこと、ということで逃げて、メソッドに
ついての簡単な説明にとどめる。ちなみに、ノルムは全て無限大ノルムで
実装した。

\subsection{誤差評価に必要なパラメータの値(coefficient)}
今回の場合、計算機で求めた $A$ の逆行列を $\widetilde{A}^{-1}$ と
すると、最大誤差評価をするためには
$$
\Vert \widetilde{A}^{-1}A - I \Vert < 1
$$
であることが条件になるので、$\Vert \widetilde{A}^{-1}A - I \Vert$
の値を調べる coefficient を実装した。実際は
$\widetilde{A}^{-1}$ を使う必要はなく、
$\Vert B^{-1}A - I \Vert < 1$
を満たす正則行列 $B$ ならば何でもよいのだが、適当な $B$ を
探すことが難しいので今回は $\widetilde{A}^{-1}$ を用いた
(ただし計算量はそれなりに必要)。

\subsection{計算結果の最大誤差を求める(MaxError)}
連立1次方程式 $A \mathbd{x} = \mathbd{b}$ を計算機で求めた解を
$\widetilde{\mathbd{x}}$ とすると、
$\Vert B^{-1}A - I \Vert < 1$ を満たす適当な正則行列 $B$ に対して
その計算結果の誤差
$\Vert \mathbd{x} - \widetilde{\mathbd{x}} \Vert$ は
$$
\Vert \mathbd{x} - \widetilde{\mathbd{x}} \Vert \leq
\frac
{\Vert B^{-1}(A \widetilde{\mathbd{x}} - \mathbd{b})\Vert}
{1- \Vert I - B^{-1}A \Vert}
$$
と評価されるので、その最大値を求める MaxError を実装した。
$B^{-1}$ は実際には上と同じく $\widetilde{A}^{-1}$ を用いた。

\section{テスト計算プログラム(calculate.c)}
乱数で作った連立1次方程式を精度保証つきで計算するテストプログラム、
calculate.c を実装した。乱数の発生には、C標準関数 rand() よりも様々な
良い点があるとされる Mersenne Twister\cite{MersenneTwister}
を用いている。使い方は make したあと ./calculate を実行するだけである。

1回目の問題は乱数で発生させた行列を用いるが、
2回目は性質の悪い行列にするために意図的な操作を加えた行列を用いている。
1回目と2回目を比べると精度がだいぶ違うことがわかると思われる。
calculate.c の行列のサイズ(DIM)を1000などの大きな数にして試す場合は、
出力命令(prinf)の部分はコメントアウトすることをおすすめする。

\section{アーカイブ一覧}
matrix.tar.gz の中身を記す。
\begin{itemize}
 \item Makefile
 \item accuracy.c
 \item accuracy.h
 \item calculate.c
 \item matrix.pdf
 \item matrixvector.c
 \item matrixvector.h
 \item mt19937ar.c
 \item mt19937ar.h
\end{itemize}
\begin{thebibliography}{4}
 \bibitem{SimpleEquation} 連立一次方程式 \\
 http://www.fuka.info.waseda.ac.jp/\~{}kozo/suuchi/simple\_equation/simple\_equation.html
 \bibitem{Algorithm} C言語によるアルゴリズム(コメント付き) \\
 http://www.sra.co.jp/people/miyata/algorithm/
 \bibitem{LUdecomposition} Gauss の消去法と LU 分解 \\
 http://daisy.math.sci.ehime-u.ac.jp/users/tsuchiya/math/linear/gauss/index.html
 \bibitem{MersenneTwister} Mersenne Twister Home Page \\
 http://www.math.sci.hiroshima-u.ac.jp/\~{}m-mat/MT/mt.html
\end{thebibliography}
\end{document}
