プ ン引数
IMSL_EQUILIBRATE, int equilibrate (入力)
入力行列 A 因子分解 前に平衡させ に指定す equilibrate 説明
0 平衡化し い
1 平衡化す
フ : equilibrate = 0
IMSL_COLUMN_ORDERING_METHOD, Imsl_col_ordering method (入力)
分解 前に 疎性 保 ようにす 列 並び え方法 次 一 method に設定し 並 び え法 選択す
method 説明
IMSL_NATURAL 自然順序 す わ 入力行列 列順
IMSL_MMD_ATA ATA 構造体 最小次 順 IMSL_MMD_AT_PLUS_A AT + A 構造体 最小次 順
IMSL_COLAMD 列近似最小次順
IMSL_PERMC プ ン引数 IMSL_COLPERM_VECTOR ー が入力した置換ベク permc 与え
順 ベク permc 数 0,1,…,n-1 置換 す
フ : method = IMSL_COLAMD
IMSL_COLPERM_VECTOR, int permc[] (入力)
後並び え(postordering) 前 置換行列 P 定義す 長さ n 配列 こ 引数
IMSL_COLUMN_ORDERING_METHOD が method = IMSL_PERMC 用い 必 須 あ そ 以外 無視さ
IMSL_TRANSPOSE, int transpose (入力)
転置 問題 ATx = b 解く場合に指定す こ プ ン 因子分解 指定す プ
ン 一緒に使うこ が transpose 説明
0 Ax = b 解く
1 ATx = b 解く.
フ : transpose = 0
IMSL_ITERATIVE_REFINEMENT, int refine (入力) 反復改良 す う 指定す
refine 説明
0 反復改良 し い
1 反復改良 す
フ : refine = 1
IMSL_FACTOR_SOLVE, int factsol (入力)
LU 分解 け 線形方程式 解く け 両方行う 指定す
factsol 説明
0 A LU 分解 方程式 Ax = b 解も求
1 LU 分解 け 行 う LU 分解 プ ン引数
IMSL_RETURN_SPARSE_LU_FACTOR 引 渡さ
factsol 説明
入力引数 b 無視さ
2 与え たLU分解 Ax = b 解く LU分解 プ
ン引数IMSL_SUPPLY_SPARSE_LU_FACTOR 渡 さ い 入力引数 A 反復改良 条件数
計算 あ い ッ 増大因子 逆数 計算以外 無視さ
フ : factsol = 0
IMSL_DIAG_PIVOT_THRESH, double diag_pivot_thresh (入力) 対角項が ッ し 扱わ う 閾値 指定す 0.0
diag_pivot_thresh
1.0.フ : diag_pivot_thresh = 1.0
IMSL_SYMMETRIC_MODE, int symm_mode (入力)
対称 ー プ ン 使う う 指定す こ ー 入力行列 A が対角項が 主要 場合に適用さ こ 場合 ー 小さ 対角 ッ 閾値 (e.g. 0.0 あ い
0.01) プ ン IMSL_DIAG_PIVOT_THRESH 指定し (AT + A)ベー 列置換ア
(例え 列置換法 IMSL_MMD_AT_PLUS_A) 選定す symm_mode 説明
0 対称 ー プ ン 使わ い
1 対称 ー プ ン 使う
フ : symm_mode = 0
IMSL_PERFORMANCE_TUNING, int sp_ienv[] (入力)
長さ 6 配列 ー が因子分解 ア 性能 チ ーニングす 際 正 パ ータ す
i sp_ienv[i] 説明
0 パネ イ
フ : sp_ienv[0] = 10
1 ーパー ー 融合 制御す 緩和パ ータ
フ : sp_ienv[1] = 5
2 ーパー ー 最大許容 イ
フ : sp_ienv[2] = 100
3 2D ロッキングに使わ 最小 行次元
フ : sp_ienv[3] = 200
4 2D ロッキングに使わ 最小 列次元
フ : sp_ienv[4] = 40
5 A 比べた L U に対す 推定充填因子
フ : sp_ienv[5] = 20
IMSL_CSC_FORMAT, int HB_col_ptr[], int HB_row_ind[], float HB_values[] (入力)
縮さ た疎 列 ー (CSC) 係数行列 受け入 こ ー 説明 イン ロ
ダク ン 参照
IMSL_SUPPLY_SPARSE_LU_FACTOR, Imsl_f_super_lu_factor lu_factor_supplied (入力) プ ン IMSL_RETURN_SPARSE_LU_FACTOR 計算さ た入力行列 LU 分解 含
Imsl_f_super_lu_factor 型 構造体 こ 構造体に い 説明 参照 こ 構造体
内に割 当 た領域 解放す に 関数 imsl_f_superlu_factor_free 使う IMSL_RETURN_SPARSE_LU_FACTOR, Imsl_f_super_lu_factor *lu_factor_returned (出力)
入力行列 LU 分解 含 Imsl_f_super_lu_factor 型 構造体 ア こ 構造体に
い 説明 参照 こ 構造体内に割 当 た領域 解放す に 関数 imsl_f_superlu_factor_free 使う
IMSL_CONDITION, float *condition (出力)
平衡化がさ た場合そ 後 行列 条件数 逆数 推定
IMSL_PIVOT_GROWTH_FACTOR, float *recip_pivot_growth (出力) ッ 増大因子 逆数
min (j P D AD Pr r c c j) /Uj
recip_pivot_growth が 1 よ 非常に小さい場合 LU 分解 安定性 良く い IMSL_FORWARD_ERROR_BOUND, float *ferr (出力)
解ベク x に対す 推定さ 前方誤差範囲 こ プ ン 引数 IMSL_ITERATIVE_REFINEMENT 1 に設定す こ が必要 す IMSL_BACKWARD_ERROR, float *berr (出力)
解ベク x 成分毎 相対後方誤差 こ プ ン 引数 IMSL_ITERATIVE_REFINEMENT 1 に設定す こ が必要 す IMSL_RETURN_USER, float x[] (出力)
線形方程式 解 x 含 ー 割 当 長さ n 配列
説明
疎 線形方程式 Ax = b 考えます ここ A 一般 n x n 非特異矩形行列 x b 長さ n
ベク す Ａ, x, b 全 ン 実数型 す
こ に対す ウ 消去 簡単に書く 次 ように ます：
1. 角分解 計算します： Pγ Dγ A Dc Pc = LU
ここ Dγ Dc 平衡にす 正定値 対角行列 Pγ Pc 数値的安定性 確 実にし 疎 形 保持す 置換行列 す L 単位下 角行列 U 角行列 す 2. Ax = b 次 ようにし 解 ます：
1 c
( (
c 1( ( (
1 r r)))))
xA b D P U L P D bこ 最後 表式 右 左に乗 こ によ 効果的に計算さ ます：
始 に b 行 Dγ ー します 次に Pγ
(D
γ b) Dγ b 行 置換す こ 意味 します L-1 (Pγ Dγ b) 行列 L 置 換えによ 角方程式 解くこ 意味します 次に U 同様に 角方程式 解 ます関数 imsl_f_superlu フ ップ 1 行 う あ い プ ン引数
IMSL_FACTOR_SOLVE が使わ 1 に設定さ す よ 詳しく述べ Ax = b 解く前に 次 ップが行 わ ます：
1. 行列 A 平衡化します す わ 対角行列 Dγ Dc à = Dγ ADc が A よ “よ 良 い条件” あ ま A 摂動に対す A-1 比べ à 摂動に対し à -1 が敏感 い ように計算します
2. Ã 列 計算さ た L U 因子 疎 度合いが増加す ように並び えます す わ Ã
à Pc 置 換えます ここ Pc 列置換行列 す
3. Ã Pc LU 分解 計算します 数値的 安定性 ÃPc 行 因子分解 過程 ー さ
た部分 ッ によ 置換さ ます こ によ 分解 Ã := Pγ Ã Pc = LU ます LU 分
解 2D ロッキング 左方 ーパー ーダ パネ ア によ 行 わ ます こ 手 法 詳細 Demmel, Eisenstat, Gilbert et al.(1999) 参照し 下さい
4. 逆 ッ 成長因子：
min1 j n j j
A U
,
計算します ここ Ã j Uj そ 行列 Ã U j 番目 列 す
5. 行列 Ã 条件数 逆数 推定します
解 過程 こ 情報 次 ップ 使わ ます：
1. 計算さ た L U 角因子 用い Ax = b 解く
2. 解 反復改良す ここ もまた計算さ た 角因子 使います こ 反復 ニ ー ン法に同等 す
3. 解ベク x に対し 前方 後方 誤差範囲 計算す
述べたいく ップ プ す こ 設定 関数 imsl_f_superlu プ ン引数 制御 ます
関数 imsl_f_superlu 行列 A LU 分解に対し ーパー ーダ 格納方式 使います 因子分
解 構造体 Imsl_f_super_lu_factor 2 構造体に含ま ます こ 構造体 以下に示します：
typedef struct{
int nnz; /* Number of nonzeros in the matrix */
float *nzval; /* Array of nonzero values packed by column */
int *rowind; /* Array of row indices of the nonzeros */
int *colptr; /* colptr[j] stores the location in nzval[]
and rowind[] which starts column j. It has ncol+1 entries, and colptr[ncol]
points to the first free location in arrays nzval[] and rowind[]. */
} Imsl_f_hb_format;
typedef struct{
int nnz; /* Number of nonzeros in the supernodal matrix */
int nsuper; /* Index of the last supernode */
float *nzval; /* Array of nonzero values packed by column */
int *nzval_colptr; /* Array of length ncol+1; nzval_colptr[j]
stores the location in nzval which starts column j. nzval_colptr[ncol] points to the first free location in arrays nzval[] and nzval_colptr[]. */
int *rowind; /* Array of compressed row indices of rectangular supernodes */
int *rowind_colptr; /* Array of length ncol+1;
rowind_colptr[sup_to_col[s]] stores the location in rowind[] which starts all columns in supernode s, and
rowind_colptr[ncol] points to the first free location in rowind[]. */
int *col_to_sup; /* Array of length ncol+1; col_to_sup[j] is the supernode number to which column j belongs. Only the first ncol entries in col_to_sup[] are defined. */
int *sup_to_col; /* Array of length ncol+1; sup_to_col[s]
points to the starting column of the s-th supernode. Only the first nsuper+2
entries in sup_to_col[] are defined, and sup_to_col[nsuper+1] = ncol+1. */
} Imsl_f_sc_format;
typedef struct{
int nrow; /* number of rows of matrix A */
int ncol; /* number of columns of matrix A */
int equilibration_method; /* The method used to equilibrate A:
0 – No equilibration 1 – Row equilibration.
2 – Column equilibration
3 – Both row and column equilibration */
float *rowscale; /* Array of length nrow containing the row scale factors for A */
float *columnscale; /* Array of length ncol containing the column scale factors for A */
int *rowperm; /* Row permutation array of length nrow describing the row permutation matrix Pr */
int *colperm; /* Column permutation array of length ncol describing the column permutation matrix Pc */
Imsl_f_hb_format *U; /* The part of the U factor of A outside the supernodal blocks, stored in Harwell-
Boeing format */
Imsl_f_sc_format *L; /* The L factor of A, stored in supernodal format as block lower triangular matrix */
} Imsl_f_super_lu_factor;
構造体 Imsl_d_super_lu_factor そ 2 構造体も float double に Imsl_f_hb_format
Imsl_d_hb_format に Imsl_f_sc_format Imsl_d_sc_format 置 換え こ によ 同様に定義さ ます
ーパー ー 定義 そ 疎LU 分解 使用に い SuperLU Users’ guide (1999) よびJ.W.
Demmel, S. C. Eisenstat et al. (1999) 参照し 下さい
例 し SuperLU Users’ guide (1999) 次 行列 考えます：
19 0 21 21 0
12 21 0 0 0
0 12 16 0 0
0 0 0 5 21
12 12 0 0 18
A
こ 行列 imsl_f_superlu 用い 自然列順 平衡化 せ sp_ienv[1] そ フ 値
5 1 に設定し 次 LU 分解 得ます：
1.00 19.00 21.00 21.00
0.63 1.00 21.00 13.26 13.26
.
0.57 1.00 23.58 7.58
1.00 5.00 21.00
0.63 0.57 0.24 0.77 1.00 34.20
A LU
次 充填さ た行列 F (I 単位行列)：
19.00 21.00 21.00
0.63 21.00 13.26 13.26 0.57 23.58 7.58
5.00 21.00 0.63 0.57 0.24 0.77 34.20
F L U I
,
行列 A 因子 ーパー ーダ 構造
1 3 4
1 2 2 4
2 2 4
3 3
3
1 2 2 3
s u u
s s s u
s s u
s s
s
s s s s
,表わさ ます ここ si i 番目 ーパー ー 非 ロ ン ui ーパー ーダ ロック 外 U i 番目 列 非 ロ ン 示します
従 ーパー ー 格納 キー 行列 F ーパー ー 部分 次 下 ロック対角行 列 し
snode
19.00
0.63 21.00 13.26 0.57 23.58
5.00 21.00
0.63 0.57 0.24 077 34.20
L
また ーパー ー 外 部分に 角行列 し 格納さ ます：
snode
21.00 21.00 13.26 7.58
U
.
こ 構造体 Imsl_f_super_lu_factor 出力 一致します：
Equilibration method: 0 Scale vectors:
rowscale: 1.000000 1.000000 1.000000 1.000000 1.000000 columnscale: 1.000000 1.000000 1.000000 1.000000 1.000000 Permutation vectors:
colperm: 0 1 2 3 4 rowperm: 0 1 2 3 4 Harwell-Boeing matrix U:
nrow 5, ncol 5, nnz 11
nzval: 21.000000 -13.263157 7.578947 21.000000 rowind: 0 1 2 0
colptr: 0 0 0 1 4 4 Supernodal matrix L:
nrow 5, ncol 5, nnz 11, nsuper 2 nzval:
0 0 1.900000e+001
1 0 6.315789e-001
4 0 6.315789e-001
1 1 2.100000e+001
2 1 5.714286e-001
4 1 5.714286e-001
1 2 -1.326316e+001
2 2 2.357895e+001
4 2 -2.410714e-001
3 3 5.000000e+000
4 3 -7.714285e-001
3 4 2.100000e+001
4 4 3.420000e+001
nzval_colptr: 0 3 6 9 11 13 rowind: 0 1 4 1 2 4 3 4 rowind_colptr: 0 3 6 6 8 8 col_to_sup: 0 1 1 2 2 sup_to_col: 0 1 3 5
関数 imsl_f_superlu Demmel, Gilbert, Li et al SuperLU コー に基 ます 分解 解法
ップ 詳細 SuperLU User’s Guide (1999) 参照し 下さい
Copyright (c) 2003, The Regents of the University of California, through Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from U.S. Dept. of Energy)
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
(1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
(2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
(3) Neither the name of Lawrence Berkeley National Laboratory, U.S. Dept. of Energy nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
例題 1
次 6×6 疎行列 LU 分解 行 います
10 0 0 0 0 0 0 10 3 1 0 0 0 0 15 0 0 0 2 0 0 10 1 0
1 0 0 5 1 3
1 2 0 0 0 6 A
y = (1, 2, 3, 4, 5, 6)T す b1 := Ay = (10, 7, 45, 33, -34, 31)T よび b2 := ATy = (-9, 8, 39, 13, 1, 21)T す A LU 分解 疎 線形方程式 Ax = b1 よび ATx = b2 解く に使わ ます
#include <imsl.h>
int main(){
Imsl_f_sparse_elem a[] = { 0, 0, 10.0, 1, 1, 10.0, 1, 2, -3.0, 1, 3, -1.0, 2, 2, 15.0, 3, 0, -2.0, 3, 3, 10.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};
float b1[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
float b2[] = { -9.0, 8.0, 39.0, 13.0, 1.0, 21.0 };
int n = 6, nz = 15;
float *x = NULL;
x = imsl_f_superlu (n, nz, a, b1, 0);
imsl_f_write_matrix ("solution to A*x = b1", 1, n, x, 0);
imsl_free (x);
x = imsl_f_superlu (n, nz, a, b2, IMSL_TRANSPOSE, 1, 0);
imsl_f_write_matrix ("solution to A^T*x = b2", 1, n, x, 0);
imsl_free (x);
}
出力
solution to A*x = b
1 2 3 4 5 6 1 2 3 4 5 6
solution to A^T*x = b2
1 2 3 4 5 6 1 2 3 4 5 6
例題 2
こ 例題 行列 A = E(1000,10) 用い A LU 分解 異 右辺 線形方程式 解くこ
示します 最大絶対値誤差 プ ン します 計算 後 LU 分解 た に割 あ いた領域 関数 imsl_f_superlu_factor_free 解放します
#include <imsl.h>
int main(){
Imsl_f_sparse_elem *a;
Imsl_f_super_lu_factor lu_factor;
float *b, *x, *mod_five, *mod_ten;
float error_factor_solve, error_solve;
int n = 1000, c = 10;
int i, nz, index;
/* Get the coefficient matrix */
a = imsl_f_generate_test_coordinate (n, c, &nz, 0);
/* Set two different predetermined solutions */
mod_five = (float*) malloc (n*sizeof(*mod_five));
mod_ten = (float*) malloc (n*sizeof(*mod_ten));
for (i=0; i<n; i++) {
mod_five[i] = (float) (i % 5);
mod_ten[i] = (float) (i % 10);
}
/* Choose b so that x will approximate mod_five */
b = imsl_f_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_five, 0);
/* Solve Ax = b */
x = imsl_f_superlu (n, nz, a, b,
IMSL_RETURN_SPARSE_LU_FACTOR, &lu_factor, 0);
/* Compute max absolute error */
error_factor_solve = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index, 0);
imsl_free (mod_five);
imsl_free (b);
imsl_free (x);
/* Get new right hand side -- b = A * mod_ten */
b = imsl_f_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_ten, 0);
/* Use the previously computed factorization to solve Ax = b */
x = imsl_f_superlu (n, nz, a, b,
IMSL_SUPPLY_SPARSE_LU_FACTOR, lu_factor, IMSL_FACTOR_SOLVE, 2,
0);
error_solve = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_ten, IMSL_INF_NORM, &index, 0);
imsl_free (mod_ten);
imsl_free (b);
imsl_free (x);
imsl_free (a);
/* Free sparse LU structure */
imsl_f_superlu_factor_free (&lu_factor);
/* Print errors */
printf ("absolute error (factor/solve) = %e\n", error_factor_solve);
printf ("absolute error (solve) = %e\n", error_solve);
}
出力
absolute error (factor/solve) = 1.502037e-005 absolute error (solve) = 1.621246e-005
警告 ー
IMSL_ILL_CONDITIONED 入力行列 非常に悪条件 す そ L1 条件数 逆数 推定値
“rcond” = # す 解 正確 い もし ませ
重大 ー
IMSL_SINGULAR_MATRIX 入力行列 特異 す