|
|
◆ Zgesvd()
| Sub Zgesvd |
( |
Jobu As |
String, |
|
|
Jobvt As |
String, |
|
|
M As |
Long, |
|
|
N As |
Long, |
|
|
A() As |
Complex, |
|
|
S() As |
Double, |
|
|
U() As |
Complex, |
|
|
Vt() As |
Complex, |
|
|
Info As |
Long |
|
) |
| |
(シンプルドライバ) 特異値分解 (SVD) (複素行列)
- 目的
- 本ルーチンはM×N複素行列Aの特異値分解(SVD), および, 必要により左および／または右特異ベクトルを求める. SVDは次のように表される. ここで, SIGMAはmin(M, N)個の対角要素を除き0のM×N行列, UはM×Mユニタリ行列, VはN×Nユニタリ行列である. SIGMAの対角要素がAの特異値である. 特異値は非負の実数で, 降順に返される. UおよびVのはじめのmin(M, N)列はAの左および右特異ベクトルである.
本ルーチンではVではなくV^Hを返すので注意せよ.
- 引数
-
| [in] | Jobu | 行列Uの全部あるいは一部を計算するか指定する.
= "A": UのM列全てを配列U()に返す.
= "S": Uのはじめのmin(M, N)列(左特異ベクトル)を配列U()に返す.
= "O": Uのはじめのmin(M, N)列(左特異ベクトル)を配列A()に上書きする.
= "N": Uの列(左特異ベクトル)を計算しない. |
| [in] | Jobvt | 行列V^Hの全部あるいは一部を計算するか指定する.
= "A": V^HのN行全てを配列Vt()に返す.
= "S": V^Hのはじめのmin(M, N)行(右特異ベクトル)を配列Vt()に返す.
= "O": V^Hのはじめのmin(M, N)行(右特異ベクトル)を配列A()に上書きする.
= "N": V^Hの行(右特異ベクトル)を計算しない.
JobvtとJobuを両方"O"にはできない. |
| [in] | M | 入力行列Aの行数. (M >= 0) (M = 0 の場合, 処理を行わずに戻る) |
| [in] | N | 入力行列Aの列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) |
| [in,out] | A() | 配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N)
[in] M×N 行列 A.
[out] Jobu = "O": A()はU(左特異ベクトル, 列ごとに格納)のはじめのmin(M, N)列で上書きされる.
Jobvt = "O": A()はV^H(右特異ベクトル, 行ごとに格納)のはじめのmin(M, N)行で上書きされる.
Jobu <> "O" かつ Jobvt <> "O": A()の内容は壊される. |
| [out] | S() | 配列 S(LS - 1) (LS >= min(M, N))
Aの特異値 (S(i) >= S(i+1) となるように並べ替えられる). |
| [out] | U() | 配列 U(LU1 - 1, LU2 - 1) (LU1 >= M, LU2 >= M (Jobu = "A") または LU2 >= min(M, N) (Jobu = "S"))
Jobu = "A": U()にM×Mユニタリ行列Uが入る.
Jobu = "S": U()にU(左特異ベクトル, 列ごとに格納)のはじめのmin(M, N)列が入る.
Jobu = "N"あるいは"O": U()は参照されない. |
| [out] | Vt() | 配列 Vt(LVt1 - 1, LVt2 - 1) (LVt1 >= N (JobVt = "A") または LVt1 >= min(M, N) (Jobvt = "S"), LVt2 >= N)
Jobvt = "A": Vt()にN×Nユニタリ行列V^Hが入る.
Jobvt = "S": Vt()にV^H(右特異ベクトル, 行ごとに格納)のはじめのmin(M, N)行が入る.
Jobvt = "N"あるいは"O": Vt()は参照されない. |
| [out] | Info | = 0: 正常終了.
= -1: パラメータ Jobu の誤り. (Jobu <> "A", "S", "O"および"N")
= -2: パラメータ Jobvt の誤り. (Jobvt <> "A", "S", "O"および"N")
= -3: パラメータ M の誤り. (M < 0)
= -4: パラメータ N の誤り. (N < 0)
= -5: パラメータ A() の誤り.
= -6: パラメータ S() の誤り.
= -7: パラメータ U() の誤り.
= -8: パラメータ Vt() の誤り.
> 0: Zbdsqrで収束しなかった場合, Infoは中間結果の二重対角形Bの上副対角要素のうち 何個が0に収束しなかったを示す. |
- 出典
- LAPACK
- 使用例
- 行列Aの特異値, 左および右特異ベクトルを求める. ただし,
( 0.20-0.11i -0.93-0.32i 0.81+0.37i )
A = ( -0.80-0.92i -0.29+0.86i 0.64+0.51i )
( 0.71+0.59i -0.15+0.19i 0.20+0.94i )
Sub Ex_Zgesvd()
Const M = 3, N = 3
Dim A(M - 1, N - 1) As Complex, U(M - 1, N - 1) As Complex, Vt(N - 1, N - 1) As Complex
Dim S(N - 1) As Double, Info As Long
A(0, 0) = Cmplx(0.2, -0.11): A(0, 1) = Cmplx(-0.93, -0.32): A(0, 2) = Cmplx(0.81, 0.37) A(1, 0) = Cmplx(-0.8, -0.92): A(1, 1) = Cmplx(-0.29, 0.86): A(1, 2) = Cmplx(0.64, 0.51) A(2, 0) = Cmplx(0.71, 0.59): A(2, 1) = Cmplx(-0.15, 0.19): A(2, 2) = Cmplx(0.2, 0.94) Call Zgesvd("S", "A", M, N, A(), S(), U(), Vt(), Info) Debug.Print "Singular values =", S(0), S(1), S(2)
Debug.Print "U ="
Debug.Print "V^H ="
Debug.Print "Info =", Info
End Sub
Function Cmplx(R As Double, Optional I As Double=0) As Complex 複素数の作成 Function Cimag(A As Complex) As Double 複素数の虚数部 Function Creal(A As Complex) As Double 複素数の実数部 Sub Zgesvd(Jobu As String, Jobvt As String, M As Long, N As Long, A() As Complex, S() As Double, U() As Complex, Vt() As Complex, Info As Long) (シンプルドライバ) 特異値分解 (SVD) (複素行列)
- 実行結果
Singular values = 2.07084030821889 1.23513084760163 0.901483337149816
U =
-0.322093837376809 0.402732043434636 9.95958553097355E-02 0.309610800223195
0.25572227156877 0.689932737910446 -0.631621088952191 3.70329357348017E-02
-0.438560869464762 1.80488745428174E-02 -7.00863751874233E-02 0.699280401316292
0.711819740685389 -0.348707521966052
-9.32902099996816E-02 0.222663951646943
-0.55115638604919 -9.82856061041459E-02
V^H =
-0.603023132575377 0 0.366551261475641 0.394522570687146
0.663818465313821 -0 0.134961467338043 -0.149629139353428
-0.44238913491109 0 -0.297134277147366 -0.762299060909463
0.16075433957414 -0.566138903287951
0.366911928609536 -0.61977189734833
0.331437452593181 -0.158279598849071
Info = 0
|
1.9.7