関連情報
ホーム > 製品 > NAG数値計算ライブラリ > サンプルソースコード集 > 実スパース対称行列の不完全コレスキー分解 (C言語/C++)

実スパース対称行列の不完全コレスキー分解

C言語によるサンプルソースコード
使用関数名:nag_sparse_sym_chol_fac (f11jac)

Keyword: スパース, 不完全コレスキー分解

概要

本サンプルは実スパース対称行列の不完全コレスキー分解を行うC言語によるサンプルプログラムです。 本サンプルは実スパース対称行列Aを読み込み、不完全コレスキー分解を計算し、行列Aと以下に示される下三角行列Cの非ゼロ要素を出力します。

実スパース非対称線形連立方程式のデータ 

※本サンプルはNAG Cライブラリに含まれる関数 nag_sparse_sym_chol_fac() のExampleコードです。本サンプル及び関数の詳細情報は nag_sparse_sym_chol_fac のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はnag_sparse_sym_chol_fac のマニュアルページを参照)

このデータをダウンロード
nag_sparse_sym_chol_fac (f11jac) Example Program Data
 7                               n
 16                              nnz
 0 0.0                           lfill, dtol
 Nag_SparseSym_UnModFact 0.0     mic, dscale
 Nag_SparseSym_MarkPiv           pstrat
  4.   1    1
  1.   2    1
  5.   2    2
  2.   3    3
  2.   4    2
  3.   4    4
 -1.   5    1
  1.   5    4
  4.   5    5
  1.   6    2
 -2.   6    5
  3.   6    6
  2.   7    1
 -1.   7    2
 -2.   7    3
  5.   7    7        a[i-1], irow[i-1], icol[i-1], i=1,...,nnz

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に行列Aの次数(n)を指定しています。
  • 3行目に行列Aの下三角部分の非ゼロ要素の数(nnz)を指定しています。
  • 4行目の1番目のパラメータ(lfill)は値が0以上の場合、コレスキー分解で許されたフィル(fill)の最大レベルを指定しています。2番目のパラメータ(dtol)は、1番目のパラメータの値が負の値の場合、フィルイン(fill-in)を制御するためのドロップ許容値を指定しています。1番目のパラメータの値が0以上の場合は2番目のパラメータの値は参照されません。
  • 5行目には行の合計を保持するためにコレスキー分解が修正されるかどうかを示すフラグ(mic)と対角スケーリング引数(dscale)を指定しています。 "Nag_SparseSym_UnModFact"はコレスキー分解が修正されないことを意味します。
  • 6行目にピボット法(pstrat)を指定しています。"Nag_SparseSym_MarkPiv"はMarkowitz法を用いて対角ピボットが実行されることを意味します。
  • 7〜22行目に行列Aの非ゼロ要素(a)、行列Aの非ゼロ要素の行インデックス(irow)、列インデックス(icol)を指定しています。

出力結果

(本関数の詳細はnag_sparse_sym_chol_fac のマニュアルページを参照)

この出力例をダウンロード
nag_sparse_sym_chol_fac (f11jac) Example Program Results

 Original Matrix 
 n     =      7
 nnz   =     16

        1      4.0000e+00       1       1
        2      1.0000e+00       2       1
        3      5.0000e+00       2       2
        4      2.0000e+00       3       3
        5      2.0000e+00       4       2
        6      3.0000e+00       4       4
        7     -1.0000e+00       5       1
        8      1.0000e+00       5       4
        9      4.0000e+00       5       5
       10      1.0000e+00       6       2
       11     -2.0000e+00       6       5
       12      3.0000e+00       6       6
       13      2.0000e+00       7       1
       14     -1.0000e+00       7       2
       15     -2.0000e+00       7       3
       16      5.0000e+00       7       7

 Factorization
 n =      7 
 nnz =     16
 npivm =      0

       17      5.0000e-01       1       1
       18      3.3333e-01       2       2
       19      3.3333e-01       3       2
       20      2.7273e-01       3       3
       21     -5.4545e-01       4       3
       22      5.2381e-01       4       4
       23     -2.7273e-01       5       3
       24      2.6829e-01       5       5
       25      6.6667e-01       6       2
       26      5.2381e-01       6       4
       27      2.6829e-01       6       5
       28      3.4788e-01       6       6
       29     -1.0000e+00       7       1
       30      5.3659e-01       7       5
       31     -5.3455e-01       7       6
       32      9.0461e-01       7       7

        i     ipiv(i) 
        1       3
        2       4
        3       5
        4       6
        5       1
        6       2
        7       7

  • 3〜22行目には元の行列の情報が出力されています。
  • 4行目に行列Aの次数が出力されています。
  • 5行目に行列Aの非ゼロ要素の数が出力されています。
  • 7〜22行目に行列Aの要素の値、行インデックスと列インデックスが出力されています。
  • 24〜44行目にはコレスキー分解して得られた行列の情報が出力されています。
  • 25行目にコレスキー分解して得られた行列の次数が出力されています。
  • 26行目にコレスキー分解して得られた行列の非ゼロ要素の数が出力されています。
  • 27行目にコレスキー分解で修正されたピボットの数が出力されています。
  • 29〜44行目にコレスキー分解して得られた下三角行列Cの要素の値、行インデックスと列インデックスが出力されています。
  • 47〜53行目にピボットインデックスが出力されています。

ソースコード

(本関数の詳細はnag_sparse_sym_chol_fac のマニュアルページを参照)

※本サンプルソースコードはNAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法


このソースコードをダウンロード
/* nag_sparse_sym_chol_fac (f11jac) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nagf11.h>

int main(void)
{
  double dtol;
  double *a;
  double dscale;
  Integer *irow, *icol;
  Integer *ipiv, nnzc, *istr;
  Integer exit_status = 0, i, n, lfill, npivm;
  Integer nnz;
  Integer num;
  char nag_enum_arg[40];
  Nag_SparseSym_Piv pstrat;
  Nag_SparseSym_Fact mic;
  Nag_Sparse_Comm comm;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_sparse_sym_chol_fac (f11jac) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n]");

  /* Read algorithmic parameters */
  scanf("%ld", &n);
  scanf("%*[^\n]");
  scanf("%ld%*[^\n]", &nnz);
  scanf("%ld%lf%*[^\n]", &lfill, &dtol);
  scanf("%39s%lf%*[^\n]", nag_enum_arg, &dscale);
  /* nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  mic = (Nag_SparseSym_Fact) nag_enum_name_to_value(nag_enum_arg);
  scanf("%39s%*[^\n]", nag_enum_arg);
  pstrat = (Nag_SparseSym_Piv) nag_enum_name_to_value(nag_enum_arg);

  /* Allocate memory */
  num = 2 * nnz;
  ipiv = NAG_ALLOC(n, Integer);
  istr = NAG_ALLOC(n + 1, Integer);
  irow = NAG_ALLOC(num, Integer);
  icol = NAG_ALLOC(num, Integer);
  a = NAG_ALLOC(num, double);

  if (!ipiv || !istr || !irow || !icol || !a) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read the matrix a */

  for (i = 1; i <= nnz; ++i)
    scanf("%lf%ld%ld%*[^\n]", &a[i - 1], &irow[i - 1],
          &icol[i - 1]);

  /* Calculate incomplete Cholesky factorization */

  /* nag_sparse_sym_chol_fac (f11jac).
   * Incomplete Cholesky factorization (symmetric)
   */
  nag_sparse_sym_chol_fac(n, nnz, &a, &num, &irow, &icol, lfill, dtol, mic,
                          dscale, pstrat, ipiv, istr, &nnzc, &npivm, &comm,
                          &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparse_sym_chol_fac (f11jac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Output original matrix */

  printf(" Original Matrix \n");
  printf(" n     = %6ld\n", n);
  printf(" nnz   = %6ld\n\n", nnz);
  for (i = 1; i <= nnz; ++i)
    printf(" %8ld%16.4e%8ld%8ld\n", i, a[i - 1],
           irow[i - 1], icol[i - 1]);
  printf("\n");

  /* Output details of the factorization */
  printf(" Factorization\n n = %6ld \n nnz = %6ld\n", n,
         nnzc);
  printf(" npivm = %6ld\n\n", npivm);
  for (i = nnz + 1; i <= nnz + nnzc; ++i)
    printf(" %8ld%16.4e%8ld%8ld\n", i, a[i - 1],
           irow[i - 1], icol[i - 1]);

  printf("\n        i     ipiv(i) \n");
  for (i = 1; i <= n; ++i)
    printf(" %8ld%8ld\n", i, ipiv[i - 1]);

END:
  NAG_FREE(irow);
  NAG_FREE(icol);
  NAG_FREE(a);
  NAG_FREE(istr);
  NAG_FREE(ipiv);
  return exit_status;
}


Results matter. Trust NAG.

Privacy Policy | Trademarks