|
|
◆ dgesvdq()
| void dgesvdq |
( |
char |
joba, |
|
|
char |
jobp, |
|
|
char |
jobr, |
|
|
char |
jobu, |
|
|
char |
jobv, |
|
|
int |
m, |
|
|
int |
n, |
|
|
int |
lda, |
|
|
double |
a[], |
|
|
double |
s[], |
|
|
int |
ldu, |
|
|
double |
u[], |
|
|
int |
ldv, |
|
|
double |
v[], |
|
|
int * |
numrank, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
double |
rwork[], |
|
|
int |
lrwork, |
|
|
int |
iwork[], |
|
|
int |
liwork, |
|
|
int * |
info |
|
) |
| |
特異値分解 (SVD) (前処理付きQR法)
- 目的
- 本ルーチンは m x n 実行列 A の特異値分解(SVD)を求める. ただし, m >= n である. A の SVD は次のように表される.
[++] [xx] [x0] [xx]
A = U * SIGMA * V^T, [++] = [xx] * [ox] * [xx]
[++] [xx]
ただし, SIGMA は n x n 対角行列, U は m x n 正規直交行列, V は n x n 直交行列である. SIGMA の対角要素は A の特異値である. U および V の列は A のそれぞれ右および左特異ベクトルである.
- 引数
-
| [in] | joba | 計算された SVD の精度のレベルを指定する.
= 'A': 要求精度は || delta A ||_F <= f(m, n) * EPS * || A ||_F で上限が示される後退誤差を持つことに対応する. ただし, EPS はマシンイプシロン (= dlamch('Epsilon')) である. これにより, Rank revealing QR分解で計算された三角行列が EPS * ||A||_F のオーダーの基準値以下であれば, このルーチンがそれを切り捨てることができるようにする (積極的な切り捨てレベル).
= 'M': 'A' と同様であるが, 切り捨てはより穏やかで, QR分解の三角行列の対角要素にドロップがある場合にのみ許される (中間的な切り捨てレベル).
= 'H': 高精度を要求する. Rank revealing QR分解に基づくランクの決定を行わない.
= 'E': 'H' と同様であるが, 追加で列スケーリングした A の条件数の推定を行い rwork[0] に返す. n^(-1/4)*rwork[0] <= ||pinv(A_scaled)||_2 <= n^(1/4)*rwork[0]. |
| [in] | jobp | = 'P': A の行は ||A(i,:)||_infinity に関して降順に並べ替えられる. これにより, 余分なデータ移動分のコストをもって数値精度を改善することができる. 数値的安定性のために推奨される.
= 'N': 行の並べ替えを行わない. |
| [in] | jobr | = 'T': 初期並べ替え付きQR分解の後, 求められた三角行列 R の転置 R^T に dgesvd を適用する. 余分なデータ移動(行列の転置)を含む. 実験, 研究, 開発に有用である.
= 'N': 三角行列 R を dgesvd に入力する. 余分なデータ移動が少ない. |
| [in] | jobu | = 'A': m 本の左特異ベクトルすべてを計算し, 配列 u[][]に返す. u[][] の項を参照せよ.
= 'S' または 'U': min(m, n) (= n) 本の左特異ベクトルを計算し, 配列 u[][]に返す. u[][] の項を参照せよ.
= 'R': ランク数 numrank を決定し numrank 本の左特異ベクトルだけを計算し, 配列 u[][]に返す.
= 'F': n 本の左特異ベクトルを, 初期QR分解の行列 Q と (R^T, 0)^T の n 本の左特異ベクトルの積とする分解形で返す. 行の並べ替えが行われた場合, 並べ替えに関する必要な情報が iwork[n:n+m-2] に返される.
= 'N': 左特異ベクトルを計算しない. |
| [in] | jobv | = 'A' または 'V': n 本の右特異ベクトルすべてを計算し, 配列 v[][]に返す.
= 'R': ランク数 numrank を決定し numrank 本の右特異ベクトルだけを計算し, 配列 v[][]に返す. このオプションは jobu = 'R' または jobu = 'N' の場合にだけ許される. それ以外であれば不正である.
= 'N': 右特異ベクトルを計算しない. |
| [in] | m | 入力行列 A の行数. (m >= 0) (m = 0 の場合, 処理を行わずに戻る) |
| [in] | n | 入力行列 A の列数. (m >= n >= 0) (n = 0 の場合, 処理を行わずに戻る) |
| [in] | lda | 二次元配列 a[][] の整合寸法. (lda >= max(1, m)) |
| [in,out] | a[][] | 配列 a[la][lda] (la >= n)
[in] 入力行列 A.
[out] jobu != 'N' または jobv != 'N' の場合, a[][] の下三角部分にハウスホルダーベクトルが dgeqp3 によって格納される. jobu = 'F' の場合, これらのハウスホルダーベクトルと work[0:n-1] を初期の A の並べ替え付きQR分解の Q を復元するために使用することができる. u[][] の項を参照せよ. |
| [out] | s[] | 配列 s[ls] (ls >= n)
A の特異値 (s[i] >= s[i+1] となるように並んでいる). |
| [in] | ldu | 二次元配列 u[][] の整合寸法. (ldu >= max(1, m) (jobu = 'A', 'S', 'U' または 'R' の場合), ldu >= max(1, n) (jobu = 'F' の場合), ldu >= 1 (その他の場合)) |
| [out] | u[][] | 配列 u[lu][ldu] (lu >= m (jobu = 'A' の場合), lu >= n (jobu = 'S', 'U', 'R' または 'F' の場合))
jobu = 'A': u[][] に m 本の左特異ベクトルを返す.
jobu = 'S', 'U' または 'R': u[][] に先頭から n 本または numrank 本の左特異ベクトルを返す.
jobu = 'F': u[][] に左特異ベクトルを構成するのに使うことができる n x n 直交行列を返す.
jobu = 'N': u[][] は参照されない. |
| [in] | ldv | 配列 v[][] の整合寸法. (ldv >= max(1, n) (jobv = 'A', 'V' または 'R', または, joba = 'E' の場合), ldv >= 1 (その他の場合)) |
| [out] | v[][] | 配列 v[lv][ldv] (lv >= n (jobv = 'A', 'V' または 'R', または, joba = 'E' の場合)
jobv = 'A' または 'V': v[][] に n x n 直交行列 V^T を返す.
jobv = 'R': v[][] に V^T (大きい方 numrank 個の特異値の右特異ベクトルで, 行ごとに格納される) の最初の numrank 行を返す.
jobv = 'N' かつ joba = 'E': v[][] は作業領域として使用される.
jobv = 'N' かつ joba != 'E': v[][] は参照されない. |
| [out] | numrank | joba の指定に従って Rank revealing QR分解の後, 最初に決定されたランク数. jobv = 'R' かつ jobu = 'R' の場合, 先頭から numrank 個の特異値およびベクトルのみが dgesvd の呼び出しにおいて要求される. 計算された特異値に 0 のものがあれば numrank の最終的な値は減るかもしれない. |
| [out] | work[] | 配列 work[lwork]
作業領域.
lwork != -1 の場合, dgeqp3 により計算されたQR分解から Q を復元するのに必要なパラメータを work[0:n-1] に返す.
lwork, lrwork または liwork = -1 の場合, info = 0 であれば work[0] に lwork の最適値, work[1] に lwork の最小値を返す. |
| [in] | lwork | 配列 work[] のサイズ. (lwork は 2 以上でなければならない: lwork = max(2, lwork))
lwqp3 = 3*n + 1.
lwcon = 3*n.
lworq = max(1, n) (jobu = 'R', 'S' または 'U' の場合), lworq = max(1, m) (jobu = 'A' の場合).
lwsvd = max(1, 5*n).
lwlqf = max(1, n/2).
lwsvd2 = max(1, 5*(n/2)).
lworlq = max(1, n).
lwqrf = max(1, n/2). and
lworq2 = max(1, n).
とすると lwork の最小値は以下のようになる.
= max(n + lwqp3, lwsvd) (固有値のみ必要な場合).
= max(n + lwqp3, lwcon, lwsvd) (固有値およびスケーリングされた条件数の推定値が必要な場合).
= n + max(lwqp3, lwsvd, lworq) (固有値および左固有ベクトルが必要な場合).
= n + max(lwqp3, lwcon, lwsvd, lworq) (固有値, 左固有ベクトルおよびスケーリングされた条件数の推定値が必要な場合).
= n + max(lwqp3, lwsvd) (固有値および右固有ベクトルが必要な場合).
= n + max(lwqp3, lwcon, lwsvd) (固有値, 右固有ベクトルおよびスケーリングされた条件数の推定値が必要な場合).
= n + max(lwqp3, lwsvd, lworq) (jobr の値にかかわらず jobv = 'R' として特異値分解全体が必要な場合).
= n + max(lwqp3, lwcon, lwsvd, lworq) (jobr の値にかかわらず jobv = 'R' として特異値分解全体およびスケーリングされた条件数の推定値がが必要な場合).
= max(n + max(lwqp3, lwsvd, lworq), n + max(lwqp3, n/2+lwlqf, n/2+lwsvd2, n/2+lworlq, lworq)) (jobv = 'A' または 'V', jobr = 'N' として特異値分解全体が必要な場合).
= max(n + max(lwqp3, lwcon, lwsvd, lworq), n + max(lwqp3, lwcon, n/2+lwlqf, n/2+lwsvd2, n/2+lworlq, lworq)) (jobv = 'A' または 'V', jobr = 'N' として特異値分解全体およびスケーリングされた条件数の推定値が必要な場合).
= max(n + max(lwqp3, lwsvd, lworq), n + max(lwqp3, n/2+lwqrf, n/2+lwsvd2, n/2+lworq2, lworq)) (jobv = 'A' または 'V', jobr = 'T' として特異値分解全体が必要な場合).
= max(n + max(lwqp3, lwcon, lwsvd, lworq), n + max(lwqp3, lwcon, n/2+lwqrf, n/2+lwsvd2, n/2+lworq2, lworq)) (jobv = 'A' または 'V', jobr = 'T' として特異値分解全体およびスケーリングされた条件数の推定値が必要な場合).
lwork = -1 の場合, 作業領域サイズの問い合わせとみなす. 配列 work, rwork および iwork のサイズの最適値および最小値の計算のみを行う. |
| [out] | rwork[] | 配列 rwork[lrwork]
- joba = 'E' の場合, 列スケーリングをした A の条件数の推定値を rwork[0] に返す. A = C * D であれば (D は対角行列, C の各列はユークリッドノルムが 1 であるとする), 列がフルランクと仮定すると n^(-1/4) * rwork[0] <= ||pinv(c)||_2 <= n^(1/4) * rwork[0] となる. そうでなければ, rwork[0] = -1 を返す.
- 初期のQR分解で求めた上三角あるいは台形行列 R を 0 として dgesvd で計算した特異値の数を rwork[1] に返す. 途中終了の場合 (行列が 0 の場合など dgesvd を呼ばなかった), rwork[1] = -1 を返す.
lwork, lrwork または liwork = -1 の場合, info = 0 であれば rwork[0] に lrwork の最小値を返す.
|
| [in] | lrwork | 配列 rwork[] のサイズ. (lrwork >= max(2, m) (jobp = 'P' の場合). lrwork >= 2 (その他の場合).)
lwork = -1 の場合, 作業領域サイズの問い合わせとみなす. 配列 work, rwork および iwork のサイズの最適値および最小値の計算のみを行う. |
| [out] | iwork[] | 配列 iwork[liwork]
iwork[0:n-1] に Rank revealing QR分解の列ピボット選択の置換情報を返す.
jobp = 'P' の場合, iwork[n:n+m-2] に行ピボット選択で使われる行入れ替え順序のインデックスを返す. これらは, jobu = 'F' の場合に左特異ベクトルを復元するために使用することができる.
lwork, lrwork または liwork = -1 の場合, info = 0 であれば iwork[0] に liwork の最小値を返す. |
| [in] | liwork | 配列 iwork[] のサイズ. (liwork は 1 以上でなければならない: liwork = max(1, liwork))
liwork >= n + m - 1 (jobp = 'P' かつ joba != 'E' の場合).
liwork >= n (jobp = 'N' かつ joba != 'E' の場合).
liwork >= n + m - 1 + n (jobp = 'P' かつ joba = 'E' の場合).
liwork >= n + n (jobp = 'N' かつ joba = 'E' の場合).
lwork = -1 の場合, 作業領域サイズの問い合わせとみなす. 配列 work, rwork および iwork のサイズの最適値および最小値の計算のみを行う. |
| [out] | info | = 0: 正常終了
= -1: 入力パラメータ joba の誤り (joba != 'A', 'M', 'H' および 'E')
= -2: 入力パラメータ jobp の誤り (jobp != 'P' および 'N')
= -3: 入力パラメータ jobr の誤り (jobr != 'T' および 'N')
= -4: 入力パラメータ jobu の誤り (jobu != 'A', 'S', 'R', 'F' および 'N')
= -5: 入力パラメータ jobv の誤り (jobv != 'A', 'V', 'R' および 'N')
= -6: 入力パラメータ m の誤り (m < 0)
= -7: 入力パラメータ n の誤り (n < 0 or n > m)
= -8: 入力パラメータ lda の誤り (lda < max(1, m))
= -9: 入力パラメータ a の誤り (a[][] に NAN の要素があった)
= -11: 入力パラメータ ldu の誤り (ldu が小さすぎる)
= -13: 入力パラメータ ldvt の誤り (ldvt が小さすぎる)
= -17: 入力パラメータ lwork の誤り (lwork が小さすぎる)
= -19: 入力パラメータ lrwork の誤り (lrwork が小さすぎる)
= -21: 入力パラメータ liwork の誤り (liwork が小さすぎる)
> 0: dbdsqr が収束しなかった. 中間結果の二重対角形 B (dgesvd で計算される) の上副対角要素のうち info 個が 0 に収束しなかった. |
- 出典
- LAPACK
|