16. 疎行列の線形計算

16.1 概要

行列の要素のうち多くが 0 であるような行列を疎行列といいます. そのような行列は偏微分方程式を数値的に解くための連立一次方程式の係数などとしてよく現われます.

疎行列を通常の二次元配列で表すと, 0 を格納するために多くのメモリが消費され, 計算時には 0 との演算に多くの CPU 時間が使われてしまいます. それを避けるために疎行列専用の格納形式が使われます.

疎行列の線形計算では, 格納形式が異なる他に, 適用される問題の性質上大規模 (例えば, 数万元以上) であることが多いため, 適するアルゴリズムも異なってきます. そのため, 従来の線形計算プログラムはそのまま使うことはできません. 具体的には, 連立一次方程式では反復法 (クリロフ部分空間法), 固有値・固有ベクトルではランチョス法やアーノルディ法が使われることが多くなります.

以下, 疎行列計算の例題, XLPack で使われている疎行列の格納形式, 基本演算プログラム, 疎行列のチェック, 連立一次方程式の解法, ファイル入出力などについて説明します.

16.2 疎行列計算の例題

5 点差分近似を使ってラプラス方程式の数値解を求める例について説明します.

\(x = 0 \sim 1, y = 0 \sim 1\) の正方形領域において, 次の偏微分方程式を解くものとします.
\[
-d^2u/dx^2 – d^2u/dy^2 = 0
\] ただし, 境界条件は
\[
u(x, 0) = 0, u(0, y) = 0, u(x, 1) = x, u(1, y) = y \space (ディリクレ条件)
\] とします.

このような偏微分方程式は, 例えば熱平衡状態にある板の温度分布の記述などに適用されます.

差分近似を行うために \(x, y\) 両方向に領域を \(N\) 等分します. 分割幅を \(h = 1/N\) として, 格子点 \(x = ih, y = jh (i = 0 \sim N, j = 0 \sim N)\) における \(u\) を \(u[i,j]\) で表すことにします. そして, 微分を次のように差分により近似します.
\[
d^2u/dx^2 = (u[i-1,j] – 2u[i,j] + u[i+1,j])/h^2 \\
d^2u/dy^2 = (u[i,j-1] – 2u[i,j] + u[i,j+1])/h^2 \\
\] これより, ラプラス方程式は次式により近似されます.
\[
-u[i,j-1] – u[i-1,j] + 4u[i,j] – u[i+1,j] – u[i,j+1] = 0
\] 境界条件 (ディリクレ条件) は次のように境界上の格子点の値を直接与えます.
\[
u[i,0] = 0 (i = 0 \sim N) \\
u[0,j] = 0 (j = 0 \sim N) \\
u[i,N] = ih (i = 0 \sim N) \\
u[N,j] = jh (j = 0 \sim N) \\
\] すなわち, 正方形領域のうち辺の部分は既知となり, 内側の部分だけ求めればよいことになります. そこで, \(m = N – 1\) とおいて, \(u_1 = u[1,1], u_2 = u[2,1], \dots, u_m = u[m,1], u_{m+1} = u[1,2], \dots, u_{m*m} = u[m,m]\) と記号をつけ直すことにします.

\(N = 7 (m = 6)\) の場合を図示すると次のようになります.

この場合, 64 格子点のうち 28 点は境界条件より既知で, 残り 36 点における \(\boldsymbol{u}\) の値が未知ということになり, 次の連立一次方程式に帰着されます.
\[
\boldsymbol{A}\boldsymbol{u} = \boldsymbol{b}
\] ただし, \(\boldsymbol{A}\) は \(36 \times 36\) 正方行列, \(\boldsymbol{u}\) は \(u_1, u_2, \dots, u_{36}\) を要素とするベクトル, \(\boldsymbol{b}\) は右辺ベクトルです.

上のラプラス方程式の近似式において, \(i = 1 \sim 6, j = 1 \sim 6\) とした \(36\) 本の式より \(\boldsymbol{A}\) と \(\boldsymbol{b}\) を定めることができます.
\[
\boldsymbol{A} =
\begin{pmatrix}
4&-1& & & & &-1 \\
-1& 4&-1& & & & &-1 \\
&-1& 4&-1& & & & &-1 \\
& &-1& 4&-1& & & & &-1 \\
& & &-1& 4&-1& & & & &-1 \\
& & & &-1& 4& & & & & &-1 \\
-1& & & & & & 4&-1& & & & &-1 \\
&-1& & & & &-1& 4&-1& & & & &-1 \\
& &-1& & & & &-1& 4&-1& & & & &-1 \\
& & &-1& & & & &-1& 4&-1& & & & &-1 \\
& & & &-1& & & & &-1& 4&-1& & & & &-1 \\
& & & & &-1& & & & &-1& 4& & & & & &-1 \\
& & & & & & & & & & & & & & & & & & \\
& & & & & & &\ddots& & & & &&\ddots& & & & & \\
& & & & & & & & & & & & & & & & & & \\
& & & & & & & & & &-1& & & & & & 4&-1 \\
& & & & & & & & & & &-1& & & & &-1& 4&-1 \\
& & & & & & & & & & & &-1& & & & &-1& 4&-1 \\
& & & & & & & & & & & & &-1& & & & &-1& 4&-1 \\
& & & & & & & & & & & & & &-1& & & & &-1& 4&-1 \\
& & & & & & & & & & & & & & &-1& & & & &-1& 4 \\
\end{pmatrix}
\] \[
\boldsymbol{b} =
\begin{pmatrix}
u[1,0] + u[0,1] \\
u[2,0] \\
u[3,0] \\
u[4,0] \\
u[5,0] \\
u[6,0] + u[6,1] \\
u[0,2] \\
\vdots \\
u[5,7] \\
u[6,7] + u[7,6] \\
\end{pmatrix}
=
\begin{pmatrix}
0 \\
0 \\
0 \\
0 \\
0 \\
h \\
0 \\
\vdots \\
5h\\
12h \\
\end{pmatrix}
\] このように \(\boldsymbol{A}\) は 1 行の中に高々 5 個の非ゼロ要素しかない対称疎行列になります.

\(N = 7\) の例を説明しましたが, より精密に解を求めたければ \(N\) を大きくしなければなりません. 例えば, \(N = 100\) とすれば, 方程式は約 1 万元になります. さらに, 3 次元の問題であれば約 100 万元になります. そうなると, 通常の線形計算プログラムは使えなくなり, 大規模疎行列用のプログラムが必要になります.

16.3 疎行列の格納形式

XLPack では次の3つの格納形式を使用しています.

(1) 圧縮行格納 (Compressed Sparse Row (CSR)) 形式
(2) 圧縮列格納 (Compressed Sparse Column (CSC)) 形式
(3) 座標 (Coordinate (COO)) 形式

(1) と (2) は計算に使われ, (3) は入出力の際に使用されます. XLPack 基本機能では (1) と (3) のみを使用します.

例として次の 3 x 4 行列を考えます.
\[
\begin{pmatrix}
1 & 0 & -2 & 4 \\
-1 & 2 & 0 & -9 \\
0 & 5 & 0 & -8 \\
\end{pmatrix}
\]

(1) 圧縮行格納 (Compressed Sparse Row (CSR)) 形式

1 本の浮動小数配列 val と 2 本の整数配列 rowptr および colind に非ゼロ要素だけを格納します. val には非ゼロ要素の値を行方向順に格納します. rowptr は各行の先頭要素の val 配列中の位置を表します. colind は対応する val の要素の列番号です. rowptr の最後には val 配列の最後 + 1 を入れておきます.

行番号と列番号は, Fortran プログラムの場合には 1 から, C プログラムの場合には 0 から始まることが多いので互換性に注意が必要です.

 
  C言語形式
    val    { 1, -2, 4, -1, 2, -9, 5, -8 }
    rowptr { 0, 3, 6, 8 }
    colind { 0, 2, 3, 0, 1, 3, 1, 3 }

  Fortran形式
    val    { 1, -2, 4, -1, 2, -9, 5, -8 }
    rowptr { 1, 4, 7, 9 }
    colind { 1, 3, 4, 1, 2, 4, 2, 4 }

 
対称行列の場合には, 対角要素とその左側のみ, あるいは, 対角要素とその右側のみを格納して使用メモリを減らすことがあります. ここでは, そのような行列を SSR 形式と表記します (一般的な呼び方ではありません).

上の例では各行の中で列番号の昇順に並んでいます. プログラムによっては行内の並び順は任意でよい場合もありますが, XLPack では計算効率のため昇順に並んでいなければなりません.

(2) 圧縮列格納 (Compressed Sparse Column (CSC)) 形式

XLPack (基本機能) では使用しないため説明は省略しますが, 非ゼロ要素の値を列方向順に格納したものです.

(3) 座標 (Coordinate (COO)) 形式

1 本の浮動小数配列 val と 2 本の整数配列 rowind および colind に非ゼロ要素だけを格納します. val には非ゼロ要素の値を格納します. rowind, colind には対応する行番号, 列番号を入れます.

行番号と列番号は, Fortran プログラムの場合には 1 から, C プログラムの場合には 0 から始まることが多いので互換性に注意が必要です.

 
  C言語形式
    val    { 1, -2, 4, -1, 2, -9, 5, -8 }
    rowind { 0, 0, 0, 1, 1, 1, 2, 2 }
    colind { 0, 2, 3, 0, 1, 3, 1, 3 }

  Fortran形式
    val    { 1, -2, 4, -1, 2, -9, 5, -8 }
    rowind { 1, 1, 1, 2, 2, 2, 3, 3 }
    colind { 1, 3, 4, 1, 2, 4, 2, 4 }

 
上の例では行番号, 列番号の昇順に並んでいますが, プログラムによっては格納の順番は任意でよい場合があります.

16.4 基本演算プログラム

XLPack (基本機能)では, CSR 形式の行列に対して次のような基本演算プログラムが用意されています.

 
  CsrDusmv  y <- αAx + βy または y <- αA^Tx + βy
  CsrDussv  Ax = b または A^Tx = b の解 (三角行列)
  CsrDusmm  C <- αAB + βC または C <- αA^TB + βC
  CsrDussm  AX = B または A^TX = B の解 (三角行列)
  SsrDusmv  y <- αAx + βy (対称行列)
  CsrTrans  疎行列の転置

 
また, 次の変換プログラムが用意されています.

 
  CooCsr    COO -> CSR
  CsrCoo    CSR -> COO
  SsrCsr    SSR (CSR 対称行列形式) -> CSR (対称なフル行列)
  CsrSsr    CSR (対称なフル行列) -> SSR (CSR 対称行列形式)
  CsrDense  CSR -> 密行列
  DenseCsr  密行列 -> CSR

16.5 疎行列のチェック

CSR 形式の行列を使用する際に注意が必要なのはポインタとインデックスの誤りです. これらは変数のアドレスに等しいので, 間違いがあると予期せぬアドレスの読み書きをしてしまいアクセスエラーを引き起こします. 通常 VBA ではポインタは扱えず, 配列の範囲も常にチェックされるのでクラッシュすることはありませんが, CSR 形式行列のポインタやインデックスの誤りでは Excel がクラッシュすることもあります. プログラム作成時には細心の注意が必要になります.

このようなエラーを少しでも減らすために疎行列のチェックプログラム CsrCheck が提供されています. エラーが完全に防げるわけではありませんが, 活用するとよいでしょう. 下に使用例を示します.

Sub Start_2()
    Const M = 3, N = 4
    Dim Val(M * N) As Double, Ptr(M) As Long, Ind(M * N) As Long
    Dim Result(9) As Long, Info As Long
    '-- Check (1)
    Debug.Print "** Check 1 **"
    Val(0) = 1: Val(1) = -2: Val(2) = 4: Val(3) = -1: Val(4) = 2
    Val(5) = -9: Val(6) = 5: Val(7) = -8
    Ptr(0) = 0: Ptr(1) = 3: Ptr(2) = 6: Ptr(3) = 8
    Ind(0) = 0: Ind(1) = 2: Ind(2) = 3: Ind(3) = 0: Ind(4) = 1
    Ind(5) = 3: Ind(6) = 1: Ind(7) = 3
    Call CsrCheck(M, N, Val(), Ptr(), Ind(), Result(), Info)
    Call PrintResult(Info, Result())
    '-- Check (2)
    Debug.Print "** Check 2 **"
    Val(0) = 0: Val(1) = 4: Val(2) = -2: Val(3) = -1: Val(4) = 2
    Val(5) = -9: Val(6) = 5: Val(7) = -8
    Ptr(0) = 0: Ptr(1) = 3: Ptr(2) = 6: Ptr(3) = 8
    Ind(0) = 0: Ind(1) = 3: Ind(2) = 2: Ind(3) = 0: Ind(4) = 1
    Ind(5) = 3: Ind(6) = 1: Ind(7) = 5
    Call CsrCheck(M, N, Val(), Ptr(), Ind(), Result(), Info)
    Call PrintResult(Info, Result())
End Sub

Sub PrintResult(Info As Long, Result() As Long)
    Debug.Print "Info =" + Str(Info)
    Debug.Print "  Sym =" + Str(Result(0))
    Debug.Print "  Nnz =" + Str(Result(1));
    Debug.Print "  Nnz(L) =" + Str(Result(2));
    Debug.Print "  Nnz(U) =" + Str(Result(3));
    Debug.Print "  Nnz(D) =" + Str(Result(4))
    Debug.Print "  Zero =" + Str(Result(5));
    Debug.Print "  Zero(D) =" + Str(Result(6))
    Debug.Print "  Empty rows =" + Str(Result(7));
    Debug.Print "  Unsorted rows =" + Str(Result(8))
    Debug.Print "  Invalid inds =" + Str(Result(9))
End Sub

出力は次のようになりました.

 
** Check 1 **
Info = 0
  Sym = 0
  Nnz = 8  Nnz(L) = 2  Nnz(U) = 4  Nnz(D) = 2
  Zero = 0  Zero(D) = 0
  Empty rows = 0  Unsorted rows = 0
  Invalid inds = 0
** Check 2 **
Info = 14
  Sym = 0
  Nnz = 8  Nnz(L) = 2  Nnz(U) = 4  Nnz(D) = 2
  Zero = 1  Zero(D) = 1
  Empty rows = 0  Unsorted rows = 1
  Invalid inds = 1

 
Check 2 では, 値が 0 の要素 1 (実害はない), 昇順になっていない行 1 (基本演算プログラムなど多くのプログラムが誤動作する), 不正なインデックス 1 (アクセスエラーを起こす) が検出されました.

16.6 連立一次方程式

16.6.1 連立一次方程式の解法

疎行列の線形計算では連立一次方程式の解法としては主に反復法 (クリロフ部分空間法) が使用されます.

一般の連立一次方程式の解法として使われる LU 分解やコレスキー分解では, 係数行列にゼロ要素があっても計算を進める過程で 0 でなくなっていきます. これに対して反復法では, 係数行列のゼロ要素は 0 のままで計算を進めることができるという特長があり, 疎行列の行列格納形式に向いています.

16.6.1.1 共役勾配 (CG) 法

係数行列 \(\boldsymbol{A}\) が正定値対称の場合には共役勾配 (CG) 法が用いられます.

\(\boldsymbol{A}\) が正定値対称行列のとき, \(\boldsymbol{x}\) が方程式
\[
\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}
\] の解であることと, \(\boldsymbol{x}\) が関数
\[
f(\boldsymbol{x}) = (1/2)(\boldsymbol{x}, \boldsymbol{Ax}) – (\boldsymbol{x}, \boldsymbol{b})
\] を最小にすることは同値になります.

そこで CG 法では, \(f(\boldsymbol{x})\) を目的関数とする非線形最適化問題を解き最小点を求めることにより連立一次方程式の解を求めます. そのためには降下法 (10 章を参照) を適用します.

ここで, \(\boldsymbol{p_k}\) を方向ベクトル, \(\alpha_k\) をステップ長として降下法の反復を行うことにします.

ステップ長は \(f(\boldsymbol{x})\) が最小になるための条件
\[
{\partial}f(\boldsymbol{x_k})/{\partial}\alpha_k = 0
\] から次のように求められます.
\[
\alpha_k = (\boldsymbol{p_k}, \boldsymbol{r_k})/(\boldsymbol{p_k}, \boldsymbol{Ap_k})
\] これを使って \(\boldsymbol{x_{k+1}}\) を更新することできます. また, 残差 \(\boldsymbol{r_k} = \boldsymbol{b} – \boldsymbol{Ax}\) も更新できます.
\[
\begin{align}
& \boldsymbol{x_{k+1}} = \boldsymbol{x_k} + \alpha_k \boldsymbol{p_k} \\
& \boldsymbol{r_{k+1}} = \boldsymbol{r_k} – \alpha_k \boldsymbol{Ap_k} \\
\end{align}
\]

方向ベクトルとして次のように最大勾配方向をとるのがよいように思われます (最急降下法).
\[
-\nabla f(\boldsymbol{x_k}) = \boldsymbol{b} – \boldsymbol{Ax} (= \boldsymbol{r_k})
\] しかし, 最急降下法では実際の計算においては収束性が悪くなることがあるため, これに修正を加えた方法が使われます. なお, 最大勾配方向は残差 \(\boldsymbol{r_k}\) に一致しています.

初期値は最大勾配方向 \(\boldsymbol{p_0} = \boldsymbol{r_0}\) から始めて, \(\boldsymbol{p_k}\) を次のように更新していきます.
\[
\begin{align}
& \beta_k = -(\boldsymbol{r_{k+1}}, \boldsymbol{Ap_k})/(\boldsymbol{p_k}, \boldsymbol{Ap_k}) \\
& \boldsymbol{p_{k+1}} = \boldsymbol{r_{k+1}} + \beta_k \boldsymbol{p_k} \\
\end{align}
\] このようにすると \(\boldsymbol{p_{k+1}}\) はすべての \(\boldsymbol{p_i} (i = 1 \sim k)\) と A-共役の関係になります. すなわち, 次式が成り立ちます.
\[
(\boldsymbol{p_i}, \boldsymbol{Ap_j}) = 0 \space (i \neq j)
\] また, 残差 \(\boldsymbol{r_{k+1}}\) はすべての \(\boldsymbol{r_i} (i = 1 \sim k)\) と直交します.
\[
(\boldsymbol{r_i}, \boldsymbol{r_j}) = 0 \space (i \neq j)
\] 従って, 計算誤差がなければ \(\boldsymbol{r_{n+1}}\) 以降は \(0\) になるはずなので, このアルゴリズムは \(n\) 回の反復で収束することになります. 実際にはそうはならず, むしろ \(n\) 回以下で収束することが多いようです.

上の計算手順をよく見ると係数行列 \(\boldsymbol{A}\) は専ら \(\boldsymbol{Ax}\) のように行列×ベクトルの演算に使われているだけです. すなわち, \(\boldsymbol{A}\) 自体は変更することがないので疎行列格納形式の行列と相性がよいことがわかります.

CG 法には別の解釈もあります. 詳細は述べませんが, クリロフ部分空間と呼ばれる空間を生成して連立一次方程式の近似解を求める解法をクリロフ部分空間法といいます. CG 法はこの方法の 1 つのバリエーションとみなすこともできます.

16.6.1.2 双共役勾配 (BICG) 法

係数行列 \(\boldsymbol{A}\) が一般行列 (非対称行列) の場合には CG 法が使えません. そのため, 種々の解法が提案されていますが, 有力な解法のひとつが双共役勾配 (BICG) 法です.

方程式 \(\boldsymbol{Ax} = \boldsymbol{b}\) の他に補助的な方程式 \(\boldsymbol{A}^T\boldsymbol{x}^* = \boldsymbol{b}^*\) を考えます.

これを使って次のような行列とベクトルを定義します.
\[
\boldsymbol{\hat{A}} =
\begin{pmatrix}
\boldsymbol{A} & \boldsymbol{0} \\
\boldsymbol{0} & \boldsymbol{A}^T \\
\end{pmatrix}
\] \[
\boldsymbol{\hat{x}} =
\begin{pmatrix}
\boldsymbol{x} \\
\boldsymbol{x}^* \\
\end{pmatrix}
\] \[
\boldsymbol{\hat{b}} =
\begin{pmatrix}
\boldsymbol{b} \\
\boldsymbol{b}^* \\
\end{pmatrix}
\] \[
\boldsymbol{\hat{H}} =
\begin{pmatrix}
\boldsymbol{0} & \boldsymbol{I} \\
\boldsymbol{I} & \boldsymbol{0} \\
\end{pmatrix}
\] 次のように \(\boldsymbol{\hat{H}}\) 内積を定義します.
\[
<\boldsymbol{\hat{x}}, \boldsymbol{\hat{y}}> \equiv (\boldsymbol{\hat{x}}, \boldsymbol{\hat{H}}\boldsymbol{\hat{y}}) = (\boldsymbol{\hat{H}}\boldsymbol{\hat{x}}, \boldsymbol{\hat{y}})
\] そうすると, \(\boldsymbol{\hat{x}}\) が方程式
\[
\boldsymbol{\hat{A}}\boldsymbol{\hat{x}} = \boldsymbol{\hat{b}}
\] の解であることと, 関数
\[
f(\boldsymbol{\hat{x}}) = (1/2)<\boldsymbol{\hat{x}}, \boldsymbol{\hat{A}}\boldsymbol{\hat{x}}> – <\boldsymbol{\hat{x}}, \boldsymbol{\hat{b}}>
\] を最小にすることは同値になります.

これに CG 法と同様の計算手順を適用して得られるのが BICG 法になります.

BICG 法もクリロフ部分空間法の 1 つのバリエーションとみなすことができます.

16.6.2 XLPack を使った連立一次方程式の解き方

XLPack (基本機能) では, 対称行列用の共役勾配 (CG) 法 (サブルーチン Cg1) と一般 (非対称) 行列用の双共役勾配 (BICG) 法 (サブルーチン Bicg1)を使用して連立一次方程式を解くことができます.

上で説明したラプラス方程式を 5 点差分近似により解く例を説明します. 係数行列が対称なので CG 法を使用します.

16.6.2.1 係数行列および右辺行列の生成

例題の説明に従って係数行列 \(\boldsymbol{A}\) および右辺ベクトル \(\boldsymbol{b}\) を用意します.

Sub Mat5DF(M As Long, Val() As Double, Ptr() As Long, Ind() As Long)
    Dim MM As Long, II As Long, I As Long, J As Long, K As Long
    MM = M * M
    K = 0
    For II = 0 To MM - 1
        I = Int(II / M)
        J = II - I * M
        Ptr(II) = K
        If I > 0 Then
            Ind(K) = II - M
            Val(K) = -1
            K = K + 1
        End If
        If J > 0 Then
            Ind(K) = II - 1
            Val(K) = -1
            K = K + 1
        End If
        Ind(K) = II
        Val(K) = 4
        K = K + 1
        If J < M - 1 Then
            Ind(K) = II + 1
            Val(K) = -1
            K = K + 1
        End If
        If I < M - 1 Then
            Ind(K) = II + M
            Val(K) = -1
            K = K + 1
        End If
    Next
    Ptr(MM) = K
End Sub

Sub Rhs5DF(M As Long, B() As Double)
    Dim H As Double, MM As Long, I As Long
    H = 1# / (M + 1)
    MM = M * M
    For I = 0 To MM - 1
        B(I) = 0
    Next
    For I = M - 1 To MM - 1 Step M
        B(I) = H * (Int(I / M) + 1)
    Next
    For I = M * (M - 1) To MM - 1
        B(I) = B(I) + H * (I - M * (M - 1) + 1)
    Next
End Sub

16.6.2.2 CG 法による解の計算

\(\boldsymbol{A}\) と \(\boldsymbol{b}\) が用意できたら Cg1 を呼び出して解を求めます. 初期値はすべて 0 としてみました. 解が求められたら適当なフォーマットで出力します.

Sub Start()
    Const N = 7, M = N - 1, MM = M * M
    Dim Val(5 * MM) As Double, Ptr(MM) As Long, Ind(5 * MM) As Long
    Dim B(MM - 1) As Double, X(MM - 1) As Double
    Dim Resid As Double, Iter As Long, Info As Long
    Dim I As Long, J As Long
    '-- Set coefficients
    Call Mat5DF(M, Val(), Ptr(), Ind())
    Call Rhs5DF(M, B())
    '-- Set initial approximation
    For I = 0 To MM - 1
        X(I) = 0
    Next
    '-- Solve equations
    Call Cg1(MM, Val(), Ptr(), Ind(), B(), X(), Info, Iter, Resid)
    '-- Output result
    If Info = 0 Then
        Call Output(N, X())
    Else
        MsgBox ("Error in Cg1: Info =" + Str(Info))
    End If
End Sub

このプログラムを実行すると次のような結果が得られました. ここでは, 境界点を含め全格子点の値を表示しています.

16.7 ファイル入出力

疎行列をテキストファイルとして格納する際によく使用される形式の一つに MM 形式 (Matrix Market Exchange Format) があります. テストマトリックスを流通させる際などに使用されています.

XLPack (基本機能) では, MM 形式ファイルに対して次のような入出力プログラムが用意されています.

 
  MMRead      Matrix Market形式ファイルの読み込み
  MMReadInfo  Matrix Market形式ファイルの行列情報の読み込み
  MMWrite     Matrix Market形式ファイルへの書き込み

付録. 有限要素法による計算例

XLPack では疎行列計算機能の一部として偏微分方程式関連のプログラムが含まれています (現時点では実験バージョンのため将来変更される可能性があります). その中の有限要素法プログラムを使用して上のラプラス方程式の例題を解いてみます.

三角要素による単純メッシュを使った例です. 格子点, 有限要素および境界要素の情報を Fem2p に入力すると, CSR 形式の有限要素行列と右辺ベクトルを出力します. これは対称行列なので CG 法で解けば偏微分方程式の数値解が求められます.

Sub Start()
    Const Nx = 7, Ny = 7, N = (Nx + 1) * (Ny + 1), Ne = 2 * Nx * Ny
    Const Nb = 2 * (Nx + Ny), Nb2 = 0, LdKnc = 4, LdKs = 3
    Dim X(N - 1) As Double, Y(N - 1) As Double
    Dim P(N - 1) As Double, Q(N - 1) As Double, F(N - 1) As Double
    Dim Knc(LdKnc - 1, Ne - 1) As Long, Ks(LdKs - 1, Nb - 1) As Long
    Dim Lb(Nb - 1) As Long, Ib(Nb - 1) As Long, Bv(Nb - 1) As Double
    Dim Ks2() As Long, Alpha() As Double, Beta() As Double
    Dim Val(20 * N) As Double, Ptr(N) As Long, Ind(20 * N) As Long
    Dim B(N - 1) As Double, U(N - 1) As Double
    Dim Ux As Double, Err As Double
    Dim Info As Long
    Dim I As Long, J As Long
    '-- Set mesh data
    Call SetData(Nx, Ny, N, Ne, Nb, X(), Y(), Knc(), Ks(), Lb(), P(), Q(), F(), Ib(), Bv())
    '-- Assemble FEM matrix
    Call Fem2p(N, Ne, X(), Y(), Knc(), P(), Q(), F(), Nb, Ib(), Bv(), Nb2, Ks2(), Alpha(), Beta(), Val(), Ptr(), Ind(), B(), Info)
    '-- Solve equation
    Dim Iter As Long, Res As Double
    Call Cg1(N, Val(), Ptr(), Ind(), B(), U(), Info, Iter, Res)
    Debug.Print "Info =" + Str(Info) + ", Iter =" + Str(Iter) + ", Res =" + Str(Res)
    '-- Display solution
    For I = 0 To Nx
        For J = 0 To Ny
            Cells(4 + Nx - J, 1 + I) = U((Ny + 1) * J + I)
        Next
    Next
End Sub

Sub SetData(Nx As Long, Ny As Long, N As Long, Ne As Long, Nb As Long, X() As Double, Y() As Double, Knc() As Long, Ks() As Long, Lb() As Long, P() As Double, Q() As Double, F() As Double, Ib() As Long, Bv() As Double)
    Dim Id() As Long
    Dim I As Long, J As Long, K As Long
    Call Mesh23(Nx, Ny, X(), Y(), Knc(), Ks(), Lb())
    For I = 0 To N - 1
        P(I) = 1
        Q(I) = 0
        F(I) = 0
    Next
    '-- Set boundary condition 1
    ReDim Id(N - 1)
    For I = 0 To Nb - 1
        For J = 1 To 2
            Id(Ks(J, I) - 1) = Lb(I)
        Next
    Next
    K = 0
    For I = 0 To N - 1
        If Id(I) = 1 Or Id(I) = 4 Then
            Bv(K) = 0
        ElseIf Id(I) = 2 Then
            Bv(K) = Y(I)
        ElseIf Id(I) = 3 Then
            Bv(K) = X(I)
        Else
            GoTo Continue
        End If
        Ib(K) = I + 1
        K = K + 1
Continue:
    Next
End Sub

このプログラムを実行すると差分法による解と (誤差の範囲で) 同じ結果が得られます.

偏微分方程式関連機能として VTK ファイルへの出力プログラム WriteVtkug が提供されています. これを使用して上のプログラムの最後の部分を変更します.

Sub Start2()
    Const Nx = 7, Ny = 7, N = (Nx + 1) * (Ny + 1), Ne = 2 * Nx * Ny
    Const Nb = 2 * (Nx + Ny), Nb2 = 0, LdKnc = 4, LdKs = 3
    Dim X(N - 1) As Double, Y(N - 1) As Double
    Dim P(N - 1) As Double, Q(N - 1) As Double, F(N - 1) As Double
    Dim Knc(LdKnc - 1, Ne - 1) As Long, Ks(LdKs - 1, Nb - 1) As Long
    Dim Lb(Nb - 1) As Long, Ib(Nb - 1) As Long, Bv(Nb - 1) As Double
    Dim Ks2() As Long, Alpha() As Double, Beta() As Double
    Dim Val(20 * N) As Double, Ptr(N) As Long, Ind(20 * N) As Long
    Dim B(N - 1) As Double, U(N - 1) As Double
    Dim Ux As Double, Err As Double
    Dim Info As Long
    Dim I As Long, J As Long
    '-- Set mesh data
    Call SetData(Nx, Ny, N, Ne, Nb, X(), Y(), Knc(), Ks(), Lb(), P(), Q(), F(), Ib(), Bv())
    '-- Assemble FEM matrix
    Call Fem2p(N, Ne, X(), Y(), Knc(), P(), Q(), F(), Nb, Ib(), Bv(), Nb2, Ks2(), Alpha(), Beta(), Val(), Ptr(), Ind(), B(), Info)
    '-- Solve equation
    Dim Iter As Long, Res As Double
    Call Cg1(N, Val(), Ptr(), Ind(), B(), U(), Info, Iter, Res)
    Debug.Print "Info =" + Str(Info) + ", Iter =" + Str(Iter) + ", Res =" + Str(Res)
    '-- Output solution
    Dim Z(N - 1) As Double
    Call WriteVtkug("Test.vtk", N, X(), Y(), Z(), Ne, Knc(), U(), Info)
    Debug.Print "Info =" + Str(Info)
End Sub

このプログラムを実行すると Test.vtk ファイルが得られます. vtk ファイルを外部プログラムに入力して結果を表示することができます. 例えば, ParaView プログラムで表示すると次のようになります.

XLPack(基本機能)では他にもメッシュ生成プログラム gmsh のファイルの入出力プログラム (ReadGmsh22, WriteGmsh22) が提供されます. これを使用すると gmsh で作成したいろいろな形の領域での計算をすることができます. 簡単な計算結果例を次に示します.