IMSL C ライブラリでの行列格納モード
本節で、行列（matrix)という単語は数学的オブジェクトを参照するために使用され、
配列(array) という単語は C データ構造としての表現のために使用されます。次の
配列リストの中で、IMSL C/Math/ライブラリの関数は、行列ディメンジョン値と行
列入力値のための全ての値から成る入力を必要します。これらの値はその配列の中
で、行並び順序で格納されます。
各関数は入力配列を処理して一般的には、
「結果」のポインターを返します。例えば、
線形代数連立方程式を解く際に、そのポインターはその解に対するものです。一般
的に、実数固有値問題では、そのポインターはその固有値に対するものです。 通常、
入力配列値はその関数によって変更されません。
IMSL C/Math/ライブラリの中で、配列はデータの連続するブロックへのポインター
です。 それらは、その行列の行へのポインターのポインターではありません。 一
般的な宣言は、以下の通りです。
float *a = {1, 2, 3, 4};
float b[2][2] = {1, 2, 3, 4};
float c[] = {1, 2, 3, 4};
注意：非 ANSI C を使用していて、その変数がタイプ auto の場合、上述の宜言は、
static float タイプとして宣言する必要があることにご注意ください。
一般モード
一般的な行列は、正方 n ×n 行列です。 一般的な配列のデータタイプは、float、
double、f_complex 又は d_complex にすることができます。
矩形モード
矩形行列は m ×n 行列です。 長方形配列のデータタイプは、float、double、
f_complex 又は d_complex にすることができます。
対称モード
対称行列は、AT = A のような正方 n×n 行列です（行列 AT は A の転置行列です）。
対称配列のデータタイプは、float、double にすることができます。
エルミートモード
エルミート行列は、次のような正方 n × n 行列 A です
H
T
A =A =A
行列 A は A の共役複素数で、
1
H
A ºA
T
H
は A の共役転置行列です。 エルミート行列 A = A では、エルミート配列のデ
ータタイプは、f_complex 又は d_complex にすることができます。
疎行列対応座標格納フォーマット
疎行列の非ゼロ要素だけを関数に知らせる必要があります。 疎行列対応座標格納
フォーマットは、入力の行と列の指標と共に各行列入力値を格納します。 次の 4
つの非同質データ構造体が、この概念を支持するために定義されています。
typedef struct{
int row;
int col;
float val;
} Imsl_f_sparse_elem:
typedef struct {
int row;
int col;
double val;
} Imsl_d_sparse_elem;
typedef struct {
int row;
int col;
f_complex val;
} Imsl_c_sparse_elem;
typedef struct {
int row;
int col;
d_complex val;
} Imsl_z_sparse_elem,
複素数データタイプ f_complex と d_complex の説明のための参照資料を見て
下さい。これらの構造体の唯一の違いは内部に含まれるデータタイプに違いが含ま
れる点を注意して下さい。疎行列は、これらのデータタイプの 1 つの配列を形成す
る事によって、疎行列対応座標フォーマットを受け取る関数に渡されます。この配
列の要素の数は疎行列中の非ゼロの数に等しくなる。
一例として、6 x 6 の行列を考えてみます。
2
0
0
0
0
0ù
é 2
ê0
9 - 3 -1 0
0 úú
ê
ê0
0
5
0
0
0ú
A=ê
ú
0 - 7 -1 0 ú
ê- 2 0
ê-1 0
0 - 5 1 - 3ú
ú
ê
0
0
6 úû
êë - 1 - 2 0
この行列 A は 15 の非ゼロ要素を持ち、その疎行列対応座標表現は次の様になりま
す。
行
列
値
0
0
2
1 1 1
1 2 3
9 -3 -1
2 3
2 0
5 -2
3 3 4
3 4 0
-7 -1 -1
4
3
-5
4 4
4 5
1 -3
5 5
0 1
-1 -2
5
5
6
この表現は順序に依存しないので、以下の形式で表現できます。
行
列
値
5 4 3
0 0 0
-1 -1 -2
0 5
0 1
2 -2
1
1
9
2 1 4
2 2 3
5 -3 -5
3 1
3 3
-7 -1
4
4
1
3
4
-1
5
5
6
4
5
-3
例えば、Imsl_f_sparse_elem の配列を初期化するために使用することができる
方法がいくつか存在します。 以下のプログラム部分を考えてみます。
#include <imsl.h>
main()
{
Imsl_f_sparse_elem a[] = {
{0, 0, 2.0},
{1, 1, 9.0},
{1, 2, -3.0},
{1, 3, -1.0},
{2, 2, 5.0},
{3, 0, -2.0},
{3, 3, -7.0},
{3, 4, -1.0},
{4, 0, -1.0},
{4, 3, -5.0},
{4, 4, 1.0},
{4, 5, -3.0},
{5, 0, -1.0},
{5, 1, -2.0},
{5, 5, 6.0} };
Imsl_f_sparse_elem b[15];
b[0].row
b[1].row
b[2].row
b[3].row
b[4].row
b[5].row
b[6].row
=
=
=
=
=
=
=
b[0].col = 0;
b[1].col = 1;
1; b[2].col = 2;
1; b[3].col = 3;
b[4].col = 2;
3; b[5].col = 0;
b[6].col = 3;
3
b[0].val
b[l].val
b[2].val
b[3].val
b[4].val
b[5].val
b[6].val
=
=
=
=
=
=
=
2.0;
9.0;
-3.0;
-1.0;
5.0;
-2.0;
-7.0;
b[7].row = 3; b[7].col = 4;
b[8].row = 4; b[8].col = 0;
b[9].row = 4; b[9].col = 3;
b[l0].row = b[10].col = 4;
b[11].row = 4; b[11].col = 5;
b[l2].row = 5; b[12].col = 0;
b[l3].row = 5; b[13] = 1;
b[l4].row = b[14].col = 5;
b[7].val = -1;
b[8].val = -1.0;
b[9].val = -5.0;
b[l0].val = 1.0;
b[11].val = -3.0;
b[l2].val = -1.0;
b[l3].val = -2.0;
b[l4].val = 6.0;
}
a と b の両方は、疎行列 A を表現して、そして、このモジュールの中の関数は、識
別子が引数リストを通して送られたかには関係なく全く同じ結果を生成します。
疎の対称、或いは、エルミート行列は、特別なケースで、これは対角項と、上、又
は、下三角行列を格納するためだけに必要です。 一例として、5 × 5 の線形連立
方程式を考えてみます。
0
0 ù
é(4,0) (1,-1)
ê (1,1) (4,0) (1,-1)
0 úú
H=ê
ê 0
(1,1) (4,0) (1,-1)ú
ê
ú
0
(1,1) (4,0) û
ë 0
このライブラリの中のエルミートと対称正定値連立方程式解法ルーチンは、その対
角項と下 3 角行列を指定されることを期待しています。 この下 3 角行列のための
疎の対応座標形式は以下の様に与えられます。
行
列
値
0
0
(4,0)
1
1
(4,0)
2
2
(4,0)
3
3
(4,0)
1
0
(1,1)
2
1
(1,1)
3
2
(1,1)
さきほどのように、以下の形式で表現できます。
行
列
値
0
0
(4,0)
1
0
(1,1)
1
1
(4,0)
2
1
(1,1)
2
2
(4,0)
以下のプログラム部分は、a と b の両方を H に初期化します。
#include <imsl.h>
main()
{
Imsl_c_sparse_elem a[] = {
{0, 0, {4.0, 0.0}},
{1, 1, {4.0, 0.0}},
{2, 2, {4.0, 0.0}},
{3, 3, {4.0, 0.0}},
{1, 0, {1.0, 1.0}},
{2, 1, {1.0, 1.0}},
{3, 2, {1.0, 1.0}}
}
4
3
2
(1,1)
3
3
(4,0)
Imsl_c_sparse_elem b[7];
b[0].row = b[0].col = 0,
b[0].val = imsl_cf_convert
b[1].row = 1; b[l].col = 0;
b[l].val = imsl_cf_convert
b[2].row = b[2].col =1;
b[2].val = imsl_cf_convert
b[3].row = 2; b[3].col = 1;
b[3].val = imsl_cf_convert
b[4].row = b[4].col = 2;
b[4].val = imsl_cf_convert
b[5].row = 3; b[5].col = 2;
b[5].val = imsl_cf_convert
b[6].row = b[6].col = 3;
b[6].val = imsl_cf_convert
(4.0, 0.0);
(1.0, 1.0);
(4.0, 0.0);
(1.0, 1.0);
(4.0, 0.0);
(1.0, 1.0);
(4.0, 0.0),
}
ここで注意すべき幾つかの重要な点があります。 H は対称ではなくてむしろエル
ミート形式です。 エルミートデータを受け取るこの関数はこれを理解して、次を
仮定することで操作します。
hij = hij
IMSL C/Math/ライブラリは、正定値ではない行列の中に対称を使用することはでき
ません。 ここに意味することは、たまたま不定である対称行列はこの簡潔な対称
形式に格納することができないということです。 上と下の、両方の 3 角行列が指
定されなければならなくて、そして、疎の一般の解法ルーチンが呼ばれなければな
りません。
帯格納フォーマット
帯行列は、その非ゼロ要素の全てが主対角項に「近い」M x N 行列です。 特に、
i-j > nlca 又は j-i> nuca であれば値は Aij = 0 です。 整数 m = nlca+ nuca
+１は、全体のバンド幅です。 主対角項以外の対角項は共対角項と名づけらます。
M x N 行列は帯行列ですが、帯格納フォーマットは、非ゼロ共対角項の数が N より
それほど多くないときだけに役立ちます。
帯格納フォーマットでは、nlca 下共対角項と、nuca 上共対角項はサイズ m × N
の配列の行に格納されます。 その要素は、それらがその行列の中にあるように、
その配列の同じ列に格納されます。 バンド幅の内側の値 Aij は、位置[（i-j +
nuca +１）* n + j ]で線形の配列に格納されます。
これは行並びの、その
行列についての 2 次元概念から一次元写像の結果になります。
例えば、1 つの下と 2 つの上共対角項を持つ 5×5 の行列 A を考えます。
5
é A0, 0
êA
ê 1, 0
A=ê 0
ê
ê 0
ê 0
ë
A0,1
A0, 2
0
A1,1
A1, 2
A1,3
A2,1
A2, 2
A2, 3
0
A3, 2
A3, 3
0
0
A4, 3
0 ù
0 úú
A2, 4 ú
ú
A3, 4 ú
A4, 4 úû
帯格納フォーマットの中では、このデータが次のように配列されます。
é 0
ê 0
ê
ê A0, 0
ê
ë A1,0
0
A0, 2
A1,3
A0,1
A1, 2
A2,3
A1,1
A2, 2
A3,3
A2,1
A3, 2
A4.3
A2, 4 ù
A3, 4 úú
A4, 4 ú
ú
0 û
このデータはそれから、長さ 20 の配列の中に、隣接して行主要の順番で保存され
ます。
一例として、次の三重対角行列を考えて下さい。
é10 1 0 0 0 ù
ê 5 20 2 0 0 ú
ú
ê
A = ê 0 6 30 3 0 ú
ú
ê
ê 0 0 7 40 4 ú
êë 0 0 0 8 50úû
以下の宜言は、この行列を帯格納フォーマットで保存します。
float a[] = {
0.0, 1.0, 2.0, 3.0, 4.0,
10.0, 20.0, 30.0, 40.0, 50.0,
5.0, 6.0, 7.0, 8.0, 0.0}；
疎行列対応座標表現の中のように、帯格納のスペース節約対称バージョンがありま
す。
一例として以下の 5 × 5 の対称の問題を見て下さい。
é A0, 0
êA
ê 0,1
A = ê A0, 2
ê
ê 0
ê 0
ë
A0,1
A0, 2
0
A1,1
A1, 2
A1,3
A1, 2
A2, 2
A2,3
A1, 3
A2,3
A3,3
0
A2, 4
A3, 4
0 ù
0 úú
A2, 4 ú
ú
A3, 4 ú
A4, 4 úû
帯対称格納フォーマットでは、そのデータは次のように配列されます。
6
é 0
ê
ê 0
êë A0, 0
0
A0, 2
A1, 2
A0,1
A1,1
A1, 2
A2, 2
A2, 3
A3, 3
A2, 4 ù
ú
A3, 4 ú
A4, 4 úû
次のエルミート形式の例は、この手順を示します。
(1,1)
0
0 ù
é (8,0) (1,1)
ê (1,-1) (8,0) (1,1)
(1,1)
0 úú
ê
H = ê (1,-1) (1,-1) (8,0) (1,1) (1,1) ú
ú
ê
(1,-1) (1,-1) (8,0) (1,1) ú
ê 0
êë 0
0
(1,-1) (1,-1) (8,0)úû
以下のプログラム断片は、バンド対称の記憶フォーマットを使用して h に H を格納し
ます。
f_complex h[ ] = {
{0.0, 0.0}, {0.0, 0.0}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0},
{0.0, 0.0}, {1.0, 1.O}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0},
{8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}
あるいは、同等の
f_complex h[15]；
h[0] = h[l] = h[5] = imsl_cf_convert（0.0, 0.0）；
h[2] = h[3] = h[4] = h[6] = h[7] = h[8] = h[9] =
imsl_cf_convert（1.0, 1.0）；
h[10] = h[11] = h[12] = h[13] = h[14] =
imsl_cf_convert（0.0, 8.0）；
帯形式と対応座標形式の間の選択方法
どのような行列も疎の対応座標形式か帯形式で格納することができます。
どちらかの形式を使うかの選択は行列の疎の分布様式に依存します。主対角近くの
帯の中に全ての非ゼロのデータと持つ行列は、恐らく帯形式の方が良いでしょう。
非ゼロ情報が多かれ少なかれ一様に行列を通して散在するならば、疎の対応座標形
式が最良の選択です。極端な例として、次の 2 つのケースが考えられます。(1)主対
角に全ての要素を持ち、そして、（0、n -１）と（n - 1, 0）入力が非ゼロの n ×
n 行列。疎の対応座標ベクトルは、n + 2 単位長さです。長さ n (2n − １)の配
列が帯表現を格納するために必要とされ、密な解法ルーチンのほぼ 2 倍の記憶域が
必要になるかも知れません。 2 番目に、全ての対角項と上対角項と下対角項の入
力が非ゼロである三重対角行列。 帯形式で長さ 3n の配列が必要です。疎の対応
座標形式の中で長さ 3n− 2 のベクトルが必要になります。しかし問題は、例えば、
32-ビットマシン上の浮動小数点精度で、対応座標形式におけるそれらの 3n − 2 単
位のそれぞれに、帯表現に必要な 3n 単位と同様な 3 倍もの多くの記憶域が必要と
されることです。これは対応座標形式の中に、その行と列指標を伴うことによりま
す。 帯記憶は、本質的に順序付けられたリストやリストの中の位置によって元の
行列の中の場所を定義することによってこの要求を避けます。
7
圧縮した疎の列（CSC）フォーマット
対応座標形式の中のデータを受け取る関数も、Harwell-Boeing 疎行列収集のための
ユーザ・ガイドに記述されたフォーマットに格納されたデータを受けることができ
ます。その枠組は列指向で、各列は疎ベクトルとして保持され、整数配列の入力値
の行指標のリストと別個になった float(double,f_complex,d_complex) 配列
の対応する値のリストによって表現されます。各列のデータは、連続的に格納され
ます。 別々の整数配列は、各列の最初の入力位置と最初の自由位置を保持する。
下三角形と対角項の中の入力値だけは、対称とエルミート行列のために格納されま
す。全ての配列は、Harwell-Boeing テストセットの、1− 基準配列に対比して、ゼロ
を基準にする。
Harnell-Boeing ユーザ・ガイドの中のように、記憶域枠組は、以下の例によって説
明されます。 5 × 5 行列
é1 - 3 0 - 1
ê0 0 - 2 0
ê
ê2 0
0
0
ê
0 -4
ê0 4
êë5 0 - 5 0
0ù
3úú
0ú
ú
0ú
6úû
は次のように、配列 colptr（最初の入力値の位置）
、rowind（行指標）と values
（非ゼロ入力値）の中に格納されます。
添字
0
1
2
3
4
5
6
7
8
9
10
colptr
0
3
5
7
9
11
rowind
0
4
2
3
0
1
4
0
3
4
1
values
1
5
2
4
-3
-2
-5
-1
-4
6
3
以下のプログラム部分は、CSC 記憶域フォーマットと対応座標表現の関係を示しま
す。
k = 0;
for (i=0; i<n; i++) {
start = colptr[i];
stop = colptr[i+1];
for (j=start; j<stop; j++) {
a[k].row = rowind[j];
a[k].col = i;
a[k++].val = values [j];
}
}
nz = k;
8
© Copyright 2024 ExpyDoc
