前回は cholmod_triplet について書きましたが、今回は疎行列 cholmod_sparse について。
SuiteSparseでは、疎行列を表現するために cholmod_sparse を使います。
その定義は以下のようになります。
--- CHOLMOD/Include/cholmod_core.h 1116行 ---
typedef struct cholmod_sparse_struct
{
size_t nrow; // 行数
size_t ncol; // 列数
size_t nzmax ; // 行列要素の最大数
void *p; // 各列までの非ゼロ成分の総数。 long * (ncol + 1)
void *i; // 非ゼロ成分の行(row)番号。 long * nzmax
int packed; // 行列を圧縮して格納するか否か
void *nz; // 各列内の非ゼロ成分数(packed == TRUE の場合のみ)。long * ncol
void *x; // データ
void *z; // 複素数の場合用
int stype; // 前回の記事参照
int itype;
int xtype;
int dtype;
int sorted; // 列についてソートするか否か
} cholmod_sparse;
cholmod_sparse オブジェクトのメモリ確保をするのが以下の関数
cholmod_sparse *
cholmod_l_allocate_sparse (size_t nrow, size_t ncol, size_t nzmax,
int sorted, int packed, int stype, int xtype,
cholmod_common *cc);
です。
この関数で、
A->i, A->x (及び複素数の場合 A->z ) に nzmax個の long*、double*の配列、
A->p に ncol + 1個のlong*の配列が確保されます。
引数 packed == FALSE なら A->nz に ncol 個のlong*が確保されます（packed == TRUE なら A->nz = NULL）。
また A->p、A->nz の各要素には 0 が代入され初期化されます。
A->i、A->x のそれぞれには非ゼロ成分の行番号、値が格納されます。
例えば
A = 0 2.2 3.3
1. 1 0 4.4
の場合、
A->i[0] = 1, A->x[0] = 1.1
A->i[1] = 0, A->x[1] = 2.2
A->i[2] = 0, A->x[2] = 3.3
A->i[3] = 1, A->x[3] = 4.4
となります。
A->p には、p[0] = 0として、各列までに何個の非ゼロ要素があるか、その総数が格納されます。
上の例の場合
一列目の最後までの非ゼロ成分の総数は1個(非ゼロ成分数は1個)、
二列目は2個(非ゼロ成分数は1個)、
三列目は4個(非ゼロ成分数は2個)
となっているので、
A->p[0] = 0
A->p[1] = 1
A->p[2] = 2
A->p[3] = 4
となります。
さらに5番目の引数 packed = FALSE を指定すると、A->nz に各列内の非ゼロ成分数が格納されます。
上の例なら
A->nz[0] = 1 // 1列目
A->nz[1] = 1 // 2列目
A->nz[2] = 2 // 3列目
となります。
以下の例では、先にあげた行列
A =
0 2.2 3.3
1.1 0 4.4
を cholmod_sparse* に格納し出力します。
#include <cholmod.h>
int
main (void)
{
cholmod_sparse *A;
long *Ai;
long *Ap;
long *Anz;
double *Ax;
int sorted;
int packed;
int stype;
cholmod_common Common, *cc;
cc = &Common;
cholmod_l_start (cc);
sorted = TRUE;
packed = FALSE;
stype = 0;
A = cholmod_l_allocate_sparse (2, 3, 4, sorted, packed,
stype, CHOLMOD_REAL, cc);
Ai = (long *) A->i;
Ap = (long *) A->p;
if (!packed) Anz = (long *) A->nz;
Ax = (double *) A->x;
Ai[0] = 1; Ax[0] = 1.1; // 1列目: 2行目の非ゼロの値が 1.1
Ap[1] = 1;
if (!packed) Anz[0] = 1; // 1列目には非ゼロは1個
Ai[1] = 0; Ax[1] = 2.2; // 2列目: 1列目の非ゼロの値が 2.2
Ap[2] = 2;
if (!packed) Anz[1] = 1; // 2列目には非ゼロは1個
Ai[2] = 0; Ax[2] = 3.3; // 3列目: 1列目の非ゼロの値が 3.3
Ai[3] = 1; Ax[3] = 4.4; // 2列目の非ゼロの値が 4.4
Ap[3] = 4;
if (!packed) Anz[2] = 2; // 3列目には非ゼロは2個
cholmod_l_write_sparse (stdout, A, NULL, NULL, cc);
cholmod_l_free_sparse (&A, cc);
cholmod_l_finish (cc);
return 0;
}
$ gcc -g -o hoge hoge.c -I/usr/include/suitesparse -lcholmod
でコンパイルして実行した結果は
%%MatrixMarket matrix coordinate integer general
2 3 4
2 1 1.1
1 2 2.2
1 3 3.3
2 3 4.4
と、前回と同じ結果になりました。