CHOLMOD その１

MUMPS と並んでもう一つ考慮に入れている Sparse Cholesky 分解用のソフトウェアが CHOLMOD になる。CHOLMOD を使うために一緒に以下のソフトウェアも必要になるので、オプションも含めて同ホームページから全部（以下の５種類）を持ってきた。

# Requires the UFconfig package.
# Requires the AMD package.
# Requires the COLAMD package.
# Optionally requires the CAMD package.
# Optionally requires the CCOLAMD package.

SDPA のためならば fill-in の最小化のための Ordering から線形方程式を解くところまで一気に進めても良いのだが、Conversion のためには一旦 Cholesky 分解時に止めて下三角行列を取り出す必要があるので、CHOLMOD を試してみようということになった。

CHOLMOD には Demo ディレクトリ以下に cholmod_simple.c というソースファイルがあるので以下のように変更してみた。

int* index;
int* index_p;
int* perm;
double* element;
int j;

cholmod_start (&c) ; /* start CHOLMOD */
A = cholmod_read_sparse (stdin, &c) ; /* read in a matrix */
cholmod_print_sparse (A, "A", &c) ; /* print the matrix */
if (A == NULL || A->stype == 0) /* A must be symmetric */
{
cholmod_free_sparse (&A, &c) ;
cholmod_finish (&c) ;
return (0) ;
}
b = cholmod_ones (A->nrow, 1, A->xtype, &c) ; /* b = ones(n,1) */
L = cholmod_analyze (A, &c) ; /* analyze */
cholmod_factorize (A, L, &c) ; /* factorize */
printf("n = %dn", L->n);
index = (int *)L->i;
index_p = (int *)L->p;
perm = (int *)L->Perm;
element = (double *)L->x;
printf("nzmax = %dn", L->nzmax);
for (j = 0; j <L->n; j++)
printf("Perm[%d] = %dn", j, perm[j]);
for (j = 0; j <L->n; j++)
printf("p[%d] = %dn", j, index_p[j]);

for (j = 0; j <L->nzmax; j++) {
printf("i[%d] = %dn", j, index[j]);
printf("x[%d] = %lfn", j, element[j]);
}