Zgssv Zgssvx

Zgssvx() Sub Zgssvx ( Trans As  String, N As  Long, Val() As  Complex, Ptr() As  Long, Ind() As  Long, Perm_c() As  Long, Perm_r() As  Long, Etree() As  Long, Equed As  String, R() As  Double, C() As  Double, B() As  Complex, X() As  Complex, Rpg As  Double, RCond As  Double, FErr() As  Double, BErr() As  Double, Optional Info As  Long, Optional Format As  Long = 0, Optional Base As  Long = -1, Optional Nrhs As  Long = 1, Optional ColPerm As  String = "A", Optional Refine As  String = "N", Optional Thresh As  Double = 1, Optional SymMode As  String = "N"  ) 連立一次方程式 A*X = B, A^T*X = B または A^H*X = B を解く (直接法) (複素行列) (SuperLU) (エキスパートドライバ) 目的 本ルーチンは LU 分解により連立一次方程式 A*X = B, A^T*X = B または A^H*X = B を解く. ただし, A は N x N 一般複素疎行列である. 方程式系の均衡化を行い, 解の誤差限界および条件数の推定値も求める. 処理ステップは以下のとおりである. A が圧縮列格納(CSC)形式の場合: 1.1. Equed = "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 を(Equed = "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. Equed = "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 を(Equed = "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] Trans 連立一次方程式の形(転置の有無)を指定する. = "N": A*X = B (転置なし). = "T": A^T*X = B (転置あり). = "C": A^H*X = B (共役転置あり). [in] N 行列 A の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in,out] Val() 配列 Val(LVal - 1) (LVal >= Nnz) (Nnz は行列 A の非ゼロ要素数) [in] 行列 A の非ゼロ要素の値. [out] Equed(入力) = "N", または, Equed(入力) = "E" であるが Equed(出力) = "N" であれば, 変更されない. Equed(入力) = "E" で Equed(出力) <> "N" であれば A は以下のようにスケーリングされる. CSC形式の場合: Equed = "R": A := diag(R)*A. Equed = "C": A := A*diag(C). Equed = "B": A := diag(R)*A*diag(C). CSR形式の場合: Equed = "R": A^T := diag(R)*A^T. Equed = "C": A^T := A^T*diag(C). Equed = "B": A^T := diag(R)*A^T*diag(C). [in] Ptr() 配列 Ptr(LPtr - 1) (LPtr >= N + 1) 行列 A の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind() 配列 Ind(LInd - 1) (LInd >= Nnz) 行列 A の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). [in,out] Perm_c() 配列 Perm_c(LPerm_c - 1) (LPerm_c >= N) [in] ColPerm = "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 - 1) (LPerm_r >= N) CSC形式の場合, ピボットの部分選択により求められる置換行列 Pr を表す行置換ベクトルである. Perm_r(i) = j は「A の行 i は Pr*A では行 j である」ことを表す. CSR形式の場合, 同様に A^T の行(A の列)の交換を表す行置換ベクトルである. [out] Etree() 配列 Etree(LEtree - 1) (LEtree >= N) Pc^T*A^T*A*Pc の消去ツリー. 注: Etree() は木構造における親ポインタのベクトルである. 木構造の頂点は整数値 0 〜 N - 1 で表される. Etree(根) = N である. [in,out] Equed 均衡化(スケーリング)の指定. [in] 均衡化するか否かを指定. = "N": 均衡化しない. = "E": 均衡化する. [out] Equed = "E" と入力した場合, 行われた均衡化を返す. = "N": 均衡化しなかった. = "R": 行均衡化を行った. すなわち, A = diag(R)*A とした. = "C": 列均衡化を行った. すなわち, A = A*diag(C) とした. = "B": 行および列均衡化を行った. すなわち, A = diag(R)*A*diag(C) とした. [out] R() 配列 R(LR - 1) (LR >= N) A または A^T の行スケーリング係数. Equed = "R" または "B" の場合, A (CSC形式の場合) または A^T (CSR形式の場合) に左から diag(R) を掛けた. Equed = "N" または "C" の場合, R() は使用されない. [out] C() 配列 C(LC - 1) (LC >= N) A または A^T の列スケーリング係数. Equed = "C" または "B" の場合, A (CSC形式の場合) または A^T (CSR形式の場合) に右から diag(C) を掛けた. Equed = "N" または "R" の場合, C() は使用されない. [in,out] B() 配列 B(LB1 - 1, LB2 - 1) (LB1 >= N, LB2 >= Nrhs) (2次元配列) または B(LB - 1) (LB >= N, Nrhs = 1) (1次元配列) [in] N x Nrhs 右辺行列 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 で上書きされる. [out] X() 配列 X(LX1 - 1, LX2 - 1) (LX1 >= N, LX2 >= Nrhs) (2次元配列) または X(LX - 1) (LX >= N, Nrhs = 1) (1次元配列) 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 - 1) (LBErr >= Nrhs) 各解ベクトル X(j) の要素に関する後退相対誤差 (すなわち, X(j) を真の解にするための A または B の任意の要素の相対変化の最小値). オプション Refine = "N" の場合, BErr(j) = 1 を返す. [out] Info (省略可) = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = i > 0 かつ <= N: U(i, i) が 0 である. 分解を完了したが, U が特異であるため解を計算できなかった. = N + 1: U は特異ではないが RCond がマシンイプシロンより小さい, これは実効精度において行列が特異であることを意味する. しかしながら, RCond の値が示すよりも計算値の精度のほうが良いことが多いため, 解と誤差限界は計算される. > N + 1: メモリ割り当てエラーが起きた. 割り当てバイト数 + N + 1 を返す. [in] Format (省略可) 行列の格納形式. (省略時 = 0) = 0: CSR 形式. = 1: CSC 形式. [in] Base (省略可) Rowptr() および Colind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Ptr(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Nrhs (省略可) 右辺の数, すなわち, 行列 B の列数. (Nrhs >= 0) (Nrhs = 0 の場合, 処理を行わずに戻る) (省略時 = 1) [in] ColPerm (省略可) フィルインを減らすための列の並べ替え方法を指定する. = "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": ユーザー指定. (省略時 = "A") [in] Refine (省略可) 反復改良を行うかどうか指定する. = "N": 反復改良を行わない. = "S": 反復改良を行う(単精度). = "D": 反復改良を行う(倍精度). = "X": 反復改良を行う(拡張精度). (省略時 = "N") [in] Thresh (省略可) ピボットとして採用する対角要素の基準値. (0 <= Thresh <= 1). (省略時 = 1) [in] SymMode (省略可) 対称モードを使用するかどうか指定する. 対称モードでは対角ピボットを優先し, A^T + A に基づく列交換アルゴリズムを使用する. = "N": 対称モードを使用しない. = "S": 対称モードを使用する. (省略時 = "N") 出典 SuperLU 使用例 連立一次方程式 Ax = B を解き, 同時にAの条件数の逆数の推定値(RCond)を求める. ただし, ( 0.78+0.16i -0.9-1.46i 0.48-1.08i ) A = ( 0.73+0.63i 1.58-1.24 -0.41-0.91i ) ( 0.23-1.37i 0.79+0.64i -0.73-1.5i ) ( 0.2126-0.2904i ) B = ( -0.3028+0.3346i ) ( -1.2905-1.0346i ) とする. Sub Ex_Zgssvx() Const N = 3, Nnz = N * N, Nrhs = 1 Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(N - 1) As Complex, X(N - 1) As Complex Dim R(N - 1) As Double, C(N - 1) As Double Dim Perm_c(N - 1) As Long, Perm_r(N - 1) As Long, Etree(N - 1) As Long Dim FErr(Nrhs - 1) As Double, BErr(Nrhs - 1) As Double Dim Equed As String, Rpg As Double, RCond As Double, Info As Long A(0) = Cmplx(0.78, 0.16): A(1) = Cmplx(-0.9, -1.46): A(2) = Cmplx(0.48, -1.08): A(3) = Cmplx(0.73, 0.63): A(4) = Cmplx(1.58, -1.24): A(5) = Cmplx(-0.41, -0.91): A(6) = Cmplx(0.23, -1.37): A(7) = Cmplx(0.79, 0.64): A(8) = Cmplx(-0.73, -1.5) Ia(0) = 0: Ia(1) = 3: Ia(2) = 6: Ia(3) = 9 Ja(0) = 0: Ja(1) = 1: Ja(2) = 2: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2: Ja(6) = 0: Ja(7) = 1: Ja(8) = 2 B(0) = Cmplx(0.2126, -0.2904): B(1) = Cmplx(-0.3028, 0.3346): B(2) = Cmplx(-1.2905, -1.0346) Equed = "E" Call Zgssvx("N", N, A(), Ia(), Ja(), Perm_c(), Perm_r(), Etree(), Equed, R(), C(), B(), X(), Rpg, RCond, FErr(), BErr(), Info) Debug.Print "X =" Debug.Print "(" + CStr(Creal(X(0))) + ", " + CStr(Cimag(X(0))) + ")" Debug.Print "(" + CStr(Creal(X(1))) + ", " + CStr(Cimag(X(1))) + ")" Debug.Print "(" + CStr(Creal(X(2))) + ", " + CStr(Cimag(X(2))) + ")" Debug.Print "RCond = " + CStr(RCond) + ", Info = " + CStr(Info) End Sub 実行結果 X = (0.59, -0.28) (-0.2, -3.99999999999999E-02) (0.24, -0.49) RCond = 0.170713882241957, Info = 0