dposv dposvx dppsv dppsvx

dppsvx() void dppsvx ( char  fact, char  uplo, int  n, int  nrhs, double  ap[], double  afp[], char *  equed, double  s[], int  ldb, double  b[], int  ldx, double  x[], double *  rcond, double  ferr[], double  berr[], double  work[], int  iwork[], int *  info  ) (エキスパートドライバ) 連立一次方程式 AX = B の解 (正定値対称行列) (圧縮形式) 目的 本ルーチンはコレスキー分解 A = U^T*U または A = L*L^T を用いて次の連立一次方程式を解く. A * X = B ここで, Aは圧縮形式のn×n正定値対称行列, また, XおよびBはn×nrhs行列である. 解の誤差限界および条件数の推定値も求められる. 解説 以下の手順で計算が行われる. fact = 'E'の場合, 連立方程式を均衡化するために次のような実数のスケーリング係数を求める. diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B 均衡化を行うかどうかは行列Aのスケールに依存するが, 均衡化が行われた場合にはAはdiag(S)*A*diag(S)により, Bはdiag(S)*Bにより上書きされる. fact = 'N'または'E'であれば, コレスキー分解によりAを(fact = 'E'であれば, 均衡化の後で)次のように分解する. A = U^T * U (uplo = 'U'の場合) A = L * L^T (uplo = 'L'の場合) ここで, Uは上三角行列, Lは下三角行列である. 係数行列のi×i首座小行列が正定値でない場合, info = iを返す. そうでない場合, 分解されたAを用いて行列Aの条件数を推定する. 条件数の逆数がマシンイプシロンより小さければ警告としてinfo = n+1を返すが, 引き続き下記のように解Xを求め誤差限界を計算する. Aの分解形を用いて連立方程式を解きXを求める. 計算された解行列に反復改良を適用して精度向上を図り, その誤差限界および後退誤差推定値を計算する. 均衡化が行われた場合, 均衡化前の元の連立方程式の解を求めるために, 行列Xに左からdiag(S)を乗ずる. 引数 [in] fact 行列Aの分解形を入力するかどうかを指定. 入力しない場合, 行列Aを分解前に均衡化するかどうかを指定する. = 'F': afp[][]にAの分解形を入力する. equed = 'Y'の場合, s[]に与えられたスケーリング係数で行列Aが均衡化済であることを示す. ap[][]およびafp[][]は変更されない. = 'N': 行列Aをafp[][]にコピーしてから分解する. = 'E': 行列Aを必要に応じて均衡化し, 次にafp[][]にコピーしてから分解する. [in] uplo = 'U': Aの上三角部分を格納. = 'L': Aの下三角部分を格納. [in] n 連立方程式の数, すなわち, 行列Aの行および列数. (n >= 0) (n = 0 の場合, 処理を行わずに戻る) [in] nrhs 右辺の数, すなわち, 行列Bの列数. (nrhs >= 0) (nrhs = 0 の場合, 処理を行わずに戻る) [in,out] ap[] 配列 ap[lap] (lap >= n(n + 1)/2) [in] 圧縮形式のn×n正定値対称行列 A. ただし, fact = 'F'かつequed = 'Y'の場合, ap[]は均衡化された行列diag(S)*A*diag(S)であること. uploに従い上三角部分あるいは下三角部分を格納する. [out] fact = 'E'かつequed = 'Y'の場合, ap[]はdiag(S)*A*diag(S)で上書きされる. その他の場合, ap[]は変更されない. [in,out] afp[] 配列 afp[lafp] (lafp >= n(n + 1)/2) [in] fact = 'F'の場合, コレスキー分解 A = U^T*U または A = L*L^T の三角行列UまたはLをap[]と同様にuploに従い入力する. ただし, equed = 'Y'の場合, afp[]は均衡化された行列diag(S)*A*diag(S)の分解形であること. [out] fact = 'N'の場合, 行列Aのコレスキー分解 A = U^T*U または A = L*L^T の三角行列UまたはLをafp[]に返す.   fact = 'E'の場合, 均衡化された行列Aのコレスキー分解 A = U^T*U または A = L*L^T の三角行列UまたはLをafp[]に返す(均衡化された行列の形についてはap[]の説明を参照せよ). [in,out] equed 均衡化の形式の指定 = 'N': 均衡化なし. = 'Y': 均衡化あり. すなわち, Aはdiag(S)*A*diag(S)で置き換えられる. [in] fact = 'F'の場合, 入力する行列Aの均衡化の有無. [out] fact = 'E'の場合, 均衡化の結果を返す. fact = 'N'の場合, 常に'N'を返す. [in,out] s[] 配列 s[ls] (ls >= n) Aのスケール係数. equed = 'N'の場合, 参照されない. [in] fact = 'F'の場合, 入力された行列Aのスケール係数. (各要素 > 0) [out] fact = 'E'の場合, 均衡化結果のスケール係数. [in] ldb 二次元配列b[][]の整合寸法. (ldb >= max(1, n)) [in,out] b[][] 配列 b[lb][ldb] (lb >= nrhs) [in] n×nrhs 右辺行列 B. [out] equed = 'N'の場合, b[][]は変更されない. equed = 'Y'の場合, b[][]はdiag(S)*Bで上書きされる. [in] ldx 二次元配列x[][]の整合寸法. (ldx >= max(1, n)) [out] x[][] 配列 x[lx][ldx] (lx >= nrhs) info = 0またはinfo = n+1の場合, 元の連立方程式のn×nrhs解行列Xを返す. equed = 'Y'の場合, a[][]およびb[][]は変更され, 均衡化後の連立方程式の解はinv(diag(S))*Xとなる. [out] rcond (均衡化した場合には均衡化後の)行列Aの(1/条件数)の推定値. rcondがマシンイプシロンより小さければ(特に rcond = 0 であれば), 実用精度において行列は特異である. これは info > 0 を返すことにより通知される. [out] ferr[] 配列 ferr[lferr] (lferr >= nrhs) 各解ベクトルX(j)(解行列Xの第j列)の前進誤差限界. X(j)に対応する真の解をXtrueとするとき, ferr[j-1]は, (X(j) - Xtrue)の要素の最大絶対値をX(j)の要素の最大絶対値で割った値の上限の推定値である. この推定値はrcondの推定値と同程度の信頼性があり, ほぼ常に真の誤差よりも大きめに推定される. [out] berr[] 配列 berr[lberr] (lberr >= nrhs) 各解ベクトルX(j)の要素に関する後退相対誤差. (すなわち, X(j)を真の解にするためのAまたはBの任意の要素の相対変化の最小値) [out] work[] 配列 work[lwork] (lwork >= 3*n) 作業領域. [out] iwork[] 配列 iwork[liwork] (liwork >= n) 作業領域. [out] info = 0: 正常終了 = -1: 入力パラメータ fact の誤り (fact != 'F', 'N'および'E') = -2: 入力パラメータ uplo の誤り (uplo != 'U'および'L') = -3: 入力パラメータ n の誤り (n < 0) = -4: 入力パラメータ nrhs の誤り (nrhs < 0) = -7: 入力パラメータ equed の誤り (fact = 'F' かつ equed != 'N'および'Y') = -8: 入力パラメータ s の誤り (s[i] <= 0 (fact = 'F' かつ equed = 'Y'のときに)) = -9: 入力パラメータ ldb の誤り (ldb < max(1, n)) = -11: 入力パラメータ ldx の誤り (ldx < max(1, n)) = i (0 < i <= n): Aのi×i首座小行列が正定値でないため分解を完了できなかった. 解は計算されなかった. rcond = 0 を返す. = n+1: Uは非特異であるが, rcondがマシンイプシロンより小さく, これは実用精度において行列が特異であることを意味する. しかしながら, rcondの値が示すよりも計算値の精度が良いことがよくあるため, 解と誤差限界は計算される. 出典 LAPACK