◆ _dggsvd3()

void _dggsvd3	(	char	jobu,
		char	jobv,
		char	jobq,
		int	m,
		int	n,
		int	p,
		int *	k,
		int *	l,
		int	lda,
		double	a[],
		int	ldb,
		double	b[],
		double	alpha[],
		double	beta[],
		int	ldu,
		double	u[],
		int	ldv,
		double	v[],
		int	ldq,
		double	q[],
		double	work[],
		int	lwork,
		int	iwork[],
		int *	info
	)

一般化特異値分解 (GSVD)

目的: 本ルーチンはm×n実行列 A および p×n実行列 B の一般化特異値分解(GSVD)を求める.
U^T*A*Q = D1*(0 R), V^T*B*Q = D2*(0 R)

ここで, U, VおよびQは直交行列である.

k+l = (A^T, B^T)^T の有効ランク数とすると, R は (k+l)×(k+l)上三角行列, D1およびD2はそれぞれ m×(k+l) および p×(k+l)「対角」行列となり, 次のような構造をもつ.

m-k-l >= 0 の場合:
k l

D1 = k (I 0)

l (0 C)

m-k-l (0 0)

k l

D2 = l (0 S)

p-l (0 0)

n-k-l k l

(0 R) = k ( 0 R11 R12)

l ( 0 0 R22)

ただし,

C = diag(alpha[k], ... , alpha[k+l-1]),

S = diag(beta[k]), ... , beta[k+l-1]),

C^2 + S^2 = I.

終了時, R は a[n-k-l〜n-1][0〜k+l-1] に格納される.

m-k-l < 0 の場合:
k m-k k+l-m

D1 = k (I 0 0)

m-k (0 C 0)

k m-k k+l-m

D2 = m-k (0 S 0)

k+l-m (0 0 I)

p-l (0 0 0)

n-k-l k m-k k+l-m

(0 R) = k ( 0 R11 R12 R13)

m-k ( 0 0 R22 R23)

k+l-m ( 0 0 0 R33)

ただし,

C = diag(alpha[k], ... , alpha[m-1]),

S = diag(beta[k], ... , beta[m-1]),

C^2 + S^2 = I.

終了時, (R11 R12 R13) は a[n-k-l〜n-1][0〜m-1] に格納される. また, R33 は

( 0 R22 R23)

b[n+m-k-l〜n-1][m-k〜l-1] に格納される.

C, S, R および必要ならば直交変換行列 U, V, Q が求められる.

行列 B がn×n行列の場合, AとBのGSVDは暗黙のうちに A*inv(B) のSVDを与える.
A*inv(B) = U*(D1*inv(D2))*V^T

(A^T,B^T)^T の列が正規直交であれば, AとBのGSVDはAとBのCS分解に等しい. さらに, GSVDを固有値問題を解くために使うことができる.
A^T*A x = lambda B^T*B x

文献によってはAとBのGSVDを次のように表していることがある.
U^T*A*X = (0 D1), V^T*B*X = (0 D2)

ここで, UとVは直交行列で, Xは特異でないものとする. D1とD2は「対角」行列である. 非特異行列Xを次のようにとることにより, GSVDの前者の形式を後者の形式に変換することができる.
X = Q*(I 0 )

(0 inv(R))

引数

[in]	jobu	= 'U': 直交行列 U を求める. = 'N': U を求めない.
[in]	jobv	='V': 直交行列 V を求める. = 'N': V を求めない.
[in]	jobq	= "Q": 直交行列 Q を求める. = 'N': Q を求めない.
[in]	m	行列 A の行数. (m >= 0)
[in]	n	行列 A および B の列数. (n >= 0)
[in]	p	行列 B の行数. (p >= 0)
[out]	k
[out]	l	kおよびlは, 「目的」に説明されているサブブロックの次数を表す. (k + l = (A^T, B^T)^T の有効ランク数)
[in]	lda	二次元配列a[][]の整合寸法. (lda >= max(1, m))
[in,out]	a[][]	配列 a[la][lda] (la >= n) [in] m×n行列 A. [out] a[][]に三角行列R, あるいは, Rの一部が入る. 詳細は「目的」を参照せよ.
[in]	ldb	二次元配列b[][]の整合寸法. (ldb >= max(1, p))
[in,out]	b[][]	配列 b[lb][ldb] (lb >= n) [in] p×n行列 B. [out] m-k-l < 0 の場合, b[][]に三角行列Rが入る. 詳細は「目的」を参照せよ.
[out]	alpha[]	配列 alpha[lalpha] (lalpha >= n)
[out]	beta[]	配列 beta[lbeta] (lbeta >= n) alpha[]およびbeta[]にAおよびBの一般化特異値が入る. alpha[0〜k-1] = 1, beta[0〜k-1] = 0, m-k-l >= 0の場合, alpha[k〜k+l-1] = C, beta[k〜k+l-1] = S, m-k-l < 0の場合, alpha[k〜m-1] = C, alpha[m〜k+l-1] = 0, beta[k〜m-1] = S, beta[m〜k+l-1] = 1, また alpha[k+l〜n-1] = 0, beta[k+l〜n-1] = 0.
[in]	ldu	二次元配列u[][]の整合寸法. (ldu >= 1 (jobu = 'N'の場合), ldu >= max(1, m) (jobu = 'U'の場合))
[out]	u[][]	配列 u[lu][ldu] (lu >= m) jobu = 'U': u[][]にm×m直交行列Uが入る. jobu = 'N': u[][]は参照されない.
[in]	ldv	二次元配列v[][]の整合寸法. (ldv >= 1 (jobv = 'N'の場合), ldv >= max(1, p) (jobv = 'V'の場合))
[out]	v[][]	配列 v[lv][ldv] (lv >= p) jobv = 'V': v[][]にp×p直交行列Vが入る. jobv = 'N': v[][]は参照されない.
[in]	ldq	二次元配列q[][]の整合寸法. (ldq >= 1 (jobq = 'N'の場合), ldq >= max(1, n) (jobq = 'Q'の場合))
[out]	q[][]	配列 q[lq][ldq] (lq >= n) jobq = 'Q': q[][]にn×n直交行列Qが入る. jobq = 'N': q[][]は参照されない.
[out]	work[]	配列 work[max(1, lwork)] 作業領域. info = 0の場合, work[0]にlworkの最適値を返す.
[in]	lwork	配列 work[]のサイズ. lwork = -1 の場合, 作業領域サイズの問い合わせとみなし, 最適サイズを求める計算だけを行い, work[0]にその値を返す.
[out]	iwork[]	配列 iwork[liwork] (liwork >= n) 作業領域. 終了時, iwork[]は並べ替え情報を格納する. より正確には, 次のループによりalpha[]を alpha[0] >= alpha[1] >= ... >= alpha[n-1] となるように並べ替える. for i = k, min(m, k+l)-1 alpha[i] と alpha[iwork[i]-1] を交換. end for
[out]	info	= 0: 正常終了 = -1: 入力パラメータ jobu の誤り (jobu != 'U'および'N') = -2: 入力パラメータ jobv の誤り (jobv != 'V'および'N') = -3: 入力パラメータ jobq の誤り (jobq != 'Q'および'N') = -4: 入力パラメータ m の誤り (m < 0) = -5: 入力パラメータ n の誤り (n < 0) = -6: 入力パラメータ p の誤り (p < 0) = -9: 入力パラメータ lda の誤り (lda < max(1, m)) = -11: 入力パラメータ ldb の誤り (ldb < max(1, p)) = -15: 入力パラメータ ldu の誤り (lduが小さすぎる) = -17: 入力パラメータ ldv の誤り (ldvが小さすぎる) = -19: 入力パラメータ ldq の誤り (ldqが小さすぎる) = -24: 入力パラメータ lwork の誤り (lworkが小さすぎる) = 1: ヤコビ型手続きが収束しなかった.

出典: LAPACK