Zgssv Zgssvx

Zgssv() Sub Zgssv ( N As  Long, Val() As  Complex, Ptr() As  Long, Ind() As  Long, B() As  Complex, RCond 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 Thresh As  Double = 1, Optional SymMode As  String = "N"  ) 連立一次方程式 A*X = B を解く (直接法) (複素行列) (SuperLU) (シンプルドライバ) 目的 本ルーチンは LU 分解により連立一次方程式 A*X = B を解く. ただし, A は N x N 一般複素疎行列である. 条件数の推定値も求める. 処理ステップは以下のとおりである. A が圧縮列格納(CSC)形式の場合: 1.1. A の列の交換を行い A*Pc を作る. ただし, Pc は置換行列である. 1.2. ピボットの部分選択を行うガウス消去法により求められる置換行列 Prを用いて A を Pr*A*Pc = L*U の形に分解する. L は対角要素が 1 で非対角要素の大きさが 1 を上限とする下三角行列である. U は上三角行列である. 1.3. 連立一次方程式 A*X = B を A の分解形を用いて解く. A が圧縮行格納(CSR)形式の場合, 上記のアルゴリズムをAの転置に適用する: 2.1. A の転置の列(A の行)の交換を行い A^T*Pc を作る. ただし, Pc は置換行列である. 2.2. ピボットの部分選択を行うガウス消去法により求められる置換行列 Pr を用いて A を Pr*A^T*Pc = L*U の形に分解する. L は対角要素が 1 で非対角要素の大きさが 1 を上限とする下三角行列である. U は上三角行列である. 1.3. 連立一次方程式 A*X = B を A の分解形を用いて解く. 引数 [in] N 行列 A の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) (Nnz は行列の非ゼロ要素数) 行列 A の非ゼロ要素の値. [in] Ptr() 配列 Ptr(LPtr - 1) (LPtr >= N + 1) 行列 A の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind() 配列 Ind(LInd - 1) (LInd >= Nnz) 行列 A の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). [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] 解行列 X. [out] Rcond 行列 A の条件数の逆数の推定値. [out] Info (省略可) = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = -10000: SuperLUプログラムで回復不能なエラーが発生した. = i > 0 かつ <= N: U(i, i)が0である. 分解を完了したが, U が特異であるため解を計算できなかった. = i > N: メモリ割り当てエラーが起きた. i - N は割り当てバイト数を表す. [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. (省略時 = "A") [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_Zgssv() Const N = 3, Nnz = N * N Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(N - 1) As Complex Dim 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) Call Zgssv(N, A(), Ia(), Ja(), B(), RCond, Info) Debug.Print "X =" Debug.Print "(" + CStr(Creal(B(0))) + ", " + CStr(Cimag(B(0))) + ")" Debug.Print "(" + CStr(Creal(B(1))) + ", " + CStr(Cimag(B(1))) + ")" Debug.Print "(" + CStr(Creal(B(2))) + ", " + CStr(Cimag(B(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.194425064735609, Info = 0