dgssv dgssvx

dgssvx() void dgssvx ( char  equil, char  col_perm, char  trans, char  refine, double  thresh, char  sym_mode, int  n, double  val[], const int  ptr[], const int  ind[], int  format, int  base, int  perm_c[], int  perm_r[], int  etree[], char *  equed, double  r[], double  c[], int  nrhs, int  ldb, double  b[], int  ldx, double  x[], double *  rpg, double *  rcond, double  ferr[], double  berr[], int *  info  ) 連立一次方程式 A*X = B または A^T*X = B を解く (疎行列, 直接法(SuperLU)) (エキスパートドライバ) 目的 本ルーチンは LU 分解により連立一次方程式 A*X = B または A^T*X = B を解く. ただし, A は M x N 一般疎行列である. 解の誤差限界および条件数の推定値も求める. 処理ステップは以下のとおりである. A が圧縮列格納(CSC)形式の場合: 1.1. equil = 'E' の場合, 方程式系を均衡化するために次のようにスケーリング係数が求められる. trans = 'N' の場合: diag(R)*A*diag(C) * inv(diag(C))*X = diag(R)*B trans = 'T' または 'C' の場合: (diag(R)*A*diag(C))^T * inv(diag(R))*X = diag(C)*B 方程式系が均衡化されるかどうかは行列 A のスケーリングに依存する. 均衡化する場合, A は diag(R)*A*diag(C) により, B は diag(R)*B (trans = 'N' の場合) または diag(C)*B (trans = 'T' または 'C' の場合) により上書きされる. 1.2. A の列の交換を行い A*Pc を作る. ただし, Pc は置換行列である. 通常スパース性は損なわれない. 1.3. 行列 A を(equil = 'E' であれば, 均衡化の後で) Pr*A*Pc = L*U とLU分解を行う. Pr はピボットの部分選択により定められる. 1.4. Pivot growth factor の逆数を求める. 1.5. U(i, i) = 0 となるものがあれば U は完全に特異であり info = i を返す. そうでなければ, 行列 A の条件数の推定値を A の分解形を使って求める. 条件数の逆数がマシンイプシロンより小さければ警告として info = n + 1 を返すが, 引き続き下記のように解 X を求め誤差限界を計算する.. 1.6. A の分解形を用いて連立方程式を解き X を求める. 1.7. refine != 'N' であれば, 計算された解行列の精度向上を図るため反復改良を適用し, その誤差限界および後退誤差推定値を計算する. 1.8. 均衡化が行われた場合, 均衡化前の元の連立方程式の解を求めるために, 行列 X に左から diag(C) (trans = 'N' の場合) または diag(R) (trans = 'T' または 'C' の場合)を掛ける. Aが圧縮行格納(CSR)形式の場合, 上記のアルゴリズムを A の転置に適用する: 2.1. equil = 'E' の場合, 方程式系を均衡化するために次のようにスケーリング係数が求められる. trans = 'N' の場合: diag(R)*A*diag(C) * inv(diag(C))*X = diag(R)*B trans = 'T' または 'C' の場合: (diag(R)*A*diag(C))^T * inv(diag(R))*X = diag(C)*B 方程式系が均衡化されるかどうかは行列 A のスケーリングに依存する. 均衡化する場合, A^T は diag(R)*A^T*diag(C) により, B は diag(R)*B (trans = 'N' の場合) または diag(C)*B (trans = 'T' または 'C' の場合) により上書きされる. 2.2. A^T の列(A の行)の交換を行い A^T*Pc を作る. ただし, Pc は置換行列である. 通常スパース性は損なわれない. 2.3. 行列 A^T を(equil = 'E' であれば, 均衡化の後で) Pr*A^T*Pc = L*U とLU分解を行う. Pr はピボットの部分選択により定められる. 2.4. Pivot growth factor の逆数を求める. 2.5. U(i, i) = 0 となるものがあれば U は完全に特異であり info = i を返す. そうでなければ, 行列 A の条件数の推定値を A^T の分解形を使って求める. 条件数の逆数がマシンイプシロンより小さければ警告として info = n + 1 を返すが, 引き続き下記のように解 X を求め誤差限界を計算する.. 2.6. A^T の分解形を用いて連立方程式を解き X を求める. 2.7. refine != 'N' であれば, 計算された解行列の精度向上を図るため反復改良を適用し, その誤差限界および後退誤差推定値を計算する. 2.8. 均衡化が行われた場合, 均衡化前の元の連立方程式の解を求めるために, 行列 X に左から diag(C) (trans = 'N' の場合) または diag(R) (trans = 'T' または 'C' の場合)を掛ける. 引数 [in] equil = 'E': 均衡化を行う (A の行および列のノルムが 1 になるようにスケーリングする). = 'N': 均衡化を行わない. [in] col_perm フィルインを減らすための列の並べ替え方法を指定する. = 'A': Approximate minimum degree (AMD) column ordering. = 'N': 並べ替えをしない (Pc = I). = 'M': A^T*A に対する Multiple minimum degree (MMD) ordering. = 'P': A^T + A に対する Multiple minimum degree (MMD) ordering. = 'U': ユーザー指定. [in] trans = 'N': A*X = B (転置なし). = 'T' または 'C': A^T*X = B (転置あり). [in] refine 反復改良を行うかどうか指定する. = 'N': 反復改良を行わない. = 'S': 反復改良を行う(単精度). = 'D': 反復改良を行う(倍精度). = 'X': 反復改良を行う(拡張精度). [in] thresh ピボットとして採用する対角要素の基準値. (0 <= thresh <= 1) [in] sym_mode 対称モードを使用するかどうか指定する. 対称モードでは対角ピボットを優先し, A^T + A に基づく列交換アルゴリズムを使用する. = 'N': 対称モードを使用しない. = 'S': 対称モードを使用する. [in] n 行列 A の次数. (n >= 0) (n = 0 の場合, 処理を行わずに戻る) [in] val[] 配列 val[lval] (lval >= nnz) CSC/CSR形式の行列 A の非ゼロ要素の値. (nnz は非ゼロ要素数) [in] ptr[] 配列 ptr[lptr] (lptr >= n + 1) CSC/CSR形式の行列 A の列／行ポインタ. [in] ind[] 配列 ind[lind] (lind >= nnz) CSC/CSR形式の行列 A の行／列インデクス. (nnz は非ゼロ要素数) [in] format 行列の格納形式. = 0: CSR 形式. = 1: CSC 形式. [in] base ptr[] および ind[] のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. [in,out] perm_c[] 配列 perm_c[lperm_c] (lperm_c >= n) [in] オプション col_perm = 'U' の場合, 置換ベクトルを入力する. その他の場合, 本引数は出力用である. [out] CSC形式の場合, 置換行列 Pc を表す列置換ベクトルである. perm_c[i] = j は「A の列 i は A*Pc では列 j である」ことを表す. CSR の場合, 同様に A^T の列(A の行)の交換を表す列置換ベクトルである. 入力された置換ベクトルは, perm_c[] の入力と Pc^T*A^T*A*Pc の消去ツリーを事後置換する置換行列との積により上書きされることがある. 消去ツリーが事後置換済であれば perm_c[] は変更されない. [out] perm_r[] 配列 perm_r[lperm_r] (lperm_r >= n) CSC形式の場合, ピボットの部分選択により求められる置換行列 Pr を表す行置換ベクトルである. perm_r[i] = j は「A の行 i は Pr*A では行 j である」ことを表す. CSR形式の場合, 同様に A^T の行(A の列)の交換を表す行置換ベクトルである. [out] etree[] 配列 etree[letree] (letree >= n) Pc^T*A^T*A*Pc の消去ツリー. 注: etree[] は木構造における親ポインタのベクトルである. 木構造の頂点は整数値 0 〜 n - 1 で表される. Etree(根) = n である. [out] equed 行われた均衡化の形式を示す. = 'N': 均衡化なし. = 'R': 行均衡化. すなわち, A に diag(R) を左から掛けた. = 'C': 列均衡化. すなわち, A に diag(C) を右から掛けた. = 'B': 行および列均衡化. すなわち, A は diag(R)*A*diag(C) である. [out] r[] 配列 r[lr] (lr >= n) A または A^T の行スケーリング係数. equed = 'R' または 'B' の場合, A (CSC の場合) または A^T (CSR の場合) に左から diag(R) を掛けた. equed = 'N' または 'C' の場合, r[] は使用されない. [out] c[] 配列 c[lc] (lc >= n) A または A^T の列スケーリング係数. equed = 'C' または 'B' の場合, A (CSC の場合) または A^T (CSR の場合) に右から diag(C) を掛けた. equed = 'N' または 'R' の場合, c[] は使用されない. [in] nrhs 右辺の数, すなわち, 行列Bの列数. (nrhs >= 0) (nrhs = 0 の場合, 処理を行わずに戻る) [in] ldb 二次元配列b[][]の整合寸法. (ldb >= max(0, n)) [in,out] b[][] 配列 b[lb][ldb] (lb >= nrhs) [in] 右辺行列 B. [out] equed = 'N' (出力された値) の場合, b[][]は変更されない. そうでない場合, b[][]は以下のように変更される. A が CSC形式の場合: trans = 'N' かつ equed = 'R' または 'B' であれば, b[][]は diag(R)*B で上書きされる. trans = 'T' または 'C' かつ equed = 'C' または 'B' であれば, b[][]は diag(C)*B で上書きされる. A が CSR形式の場合: trans = 'N' かつ equed = 'C' または 'B' であれば, b[][]は diag(C)*B で上書きされる. trans = 'T' または 'C' かつ equed = 'R' または 'B' であれば, b[][]は diag(R)*B で上書きされる. [in] ldx 二次元配列x[][]の整合寸法. (ldx >= max(0, n)) [out] x[][] 配列 x[lx][ldx] (lx >= nrhs) 解行列 X. info = 0 または info = n + 1 の場合, X は元の(均衡化前の)方程式の解行列である. equed != 'N' であれば A と B は変更されていることに注意せよ. また, 均衡化された方程式の解は inv(diag(C))*X (trans = 'N' かつ equed = 'C' または 'B' の場合) または inv(diag(R))*X (trans = 'T' または 'C' かつ equed = 'R' または 'B' の場合) である. [out] rpg Pivot growth factor の逆数 max_j(norm(A_j)/norm(U_j)). 無限ノルムが使われる. rpg が 1 に比べて非常に小さければ, LU分解は安定性に乏しいと考えられる. [out] rcond (もし均衡化されたならば, 均衡化後の)行列 A の条件数の逆数の推定値. rcond がマシンイプシロンより小さければ(特に rcond = 0 ならば), 行列は実効精度において特異である. その場合 info > 0 を返す. [out] ferr[] 配列 ferr(lferr - 1) (lferr >= nrhs) 各解ベクトル X(j)(解行列 X の第 j 列) の前進誤差限界. X(j) に対応する真の解を Xtrue とするとき, ferr[j] は (X(j) - Xtrue) の要素の最大絶対値を X(j) の要素の最大絶対値で割った値の上限の推定値である. この推定値は rcond の推定値と同程度の信頼性があり, ほぼ常に真の誤差よりも大きめに推定される. refine = 'N' の場合, ferr[j] = 1 を返す. [out] berr[] 配列 berr[lberr] (lberr >= nrhs) 各解ベクトル X(j) の要素に関する後退相対誤差 (すなわち, X(j) を真の解にするための A または B の任意の要素の相対変化の最小値). refine = 'N' の場合, berr[j] = 1 を返す. [out] info = 0: 正常終了. < 0: (-info)番目の入力パラメータの誤り. = -10000: SuperLUプログラムで回復不能なエラーが発生した. = i > 0 かつ <= n: U(i, i) が 0 である. 分解を完了したが, U が特異であるため解を計算できなかった. = n + 1: U は特異ではないが rcond がマシンイプシロンより小さい, これは実効精度において行列が特異であることを意味する. しかしながら, rcond の値が示すよりも計算値の精度のほうが良いことが多いため, 解と誤差限界は計算される. = i > n + 1: メモリ割り当てエラーが起きた. i - n - 1 は割り当てバイト数を表す. 出典 SuperLU