|
◆ dgelsd()
void dgelsd |
( |
int |
m, |
|
|
int |
n, |
|
|
int |
nrhs, |
|
|
int |
lda, |
|
|
double |
a[], |
|
|
int |
ldb, |
|
|
double |
b[], |
|
|
double |
s[], |
|
|
double |
rcond, |
|
|
int * |
rank, |
|
|
double |
work[], |
|
|
int |
lwork, |
|
|
int |
iwork[], |
|
|
int * |
info |
|
) |
| |
優決定または劣決定系連立一次方程式 Ax = b の解 (特異値分解 (SVD)) (分割統治法)
- 目的
- 本ルーチンは線形最小二乗問題の最小ノルム解を A の特異値分解(SVD)を使って求める.
2-norm(| b - A*x |) を最小化する.
A は m x n 行列で, ランク落ちしていてもよい.
いくつかの右辺ベクトル b および解ベクトル x を 1 回の呼び出しで扱うことができる. これらのベクトルは, m x nrhs 右辺行列 B および n x nrhs 解行列 X の列として格納される.
問題を以下の 3 ステップにより解く.
(1) ハウスホルダー変換により係数行列 A を二重対角形にする. これにより, 元の問題は「二重対角最小二乗問題」(BLS) に変換される.
(2) 分割統治法により BLS を解く.
(3) 元の最小二乗問題を解くためにハウスホルダー変換を逆に適用する.
A の有効ランク数は, 最大の特異値の rcond 倍よりも小さな特異値を 0 として扱うことにより決められる.
- 引数
-
[in] | m | 行列 A の行数. (m >= 0) (m = 0 の場合, rank = 0 として戻る) |
[in] | n | 行列 A の列数. (n >= 0) (n = 0 の場合, rank = 0 として戻る) |
[in] | nrhs | 右辺の数, すなわち, 行列 B および X の列数. (nrhs >= 0) |
[in] | lda | 二次元配列 a[][] の整合寸法. (lda >= max(1, m)) |
[in,out] | a[][] | 配列 a[la][lda] (la >= n)
[in] m x n 行列 A.
[out] a[][] は壊される. |
[in] | ldb | 二次元配列 b[][] の整合寸法. (ldb >= max(1, max(m, n))) |
[in,out] | b[][] | 配列 b[lb][ldb] (lb >= nrhs)
[in] m x nrhs 右辺行列 B.
[out] b[][] は n x nrhs 解行列 X により上書きされる. m >= n かつ rank = n の場合, i 列の解の残差二乗和は同じ列の n〜m-1 番要素の二乗和で与えられる. |
[out] | s[] | 配列 s[ls] (ls >= min(m, n))
行列 A の特異値 (降順).
A の2-ノルムによる条件数は s[0]/s[min(m, n)-1] である. |
[in] | rcond | rcond は A の有効ランク数を決めるために使われる.
s[i] <= rcond*s[0] となる特異値は 0 として扱われる. rcond < 0 の場合, 代わりにマシンイプシロンが使われる. |
[out] | rank | A の有効ランク数, すなわち, rcond*s[0] より大きい特異値の数. |
[out] | work[] | 配列 work[lwork]
作業領域.
終了時, info = 0 であれば, work[0] に lwork の最適値を返す. |
[in] | lwork | 配列 work[] のサイズ. (lwork は 1 以上でなければならない: lwork = max(1, lwork))
必要な作業領域の最小サイズは m, n および nrhs に依存する.
lwork が少なくとも
12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz + 1)^2 (m >= n の場合),
12*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz + 1)^2 (m < n の場合)
であればプログラムは正しく実行される.
smlsiz は ilaenv により返される値で, 計算ツリーの底の下位問題の最大サイズに等しい(通常, 約25). nlvl は次のとおりである.
nlvl = max(0, int(log2(min(m, n)/(smlsiz + 1))) + 1)
一般にパフォーマンスをよくするためには lwork をより大きくとらなければならない.
lwork = -1 の場合, 作業領域サイズの問い合わせとみなし, 配列 work[] の最適サイズ, および, 配列 iwork[] の最小サイズを求める計算だけを行い, それらの値を work[0] および iwork[0] に返す. |
[out] | iwork[] | 配列 iwork[liwork] (liwork >= max(1, 3*minmn*nlvl + 11*minmn). ただし, minmn = min(m, n).)
整数作業領域.
終了時, info = 0 であれば, iwork[0] に liwork の最小値を返す. |
[out] | info | = 0: 正常終了
= -1: 入力パラメータ m の誤り (m < 0)
= -2: 入力パラメータ n の誤り (n < 0)
= -3: 入力パラメータ nrhs の誤り (nrhs < 0)
= -4: 入力パラメータ lda の誤り (lda < max(1, m))
= -6: 入力パラメータ ldb の誤り (ldb < max(1, m, n))
= -12: 入力パラメータ lwork の誤り (lwork が小さすぎる)
= i > 0: SVD の計算アルゴリズムが収束しなかった. 中間結果の二重対角形の副対角要素のうち i 個が 0 に収束しなかった. |
- 出典
- LAPACK
|