概要
本サンプルはFortran言語によりLAPACKルーチンDGESDDを利用するサンプルプログラムです。
4x6の行列の特異値と右特異ベクトルを求めます。
計算された特異値と特異ベクトルの誤差限界近似値も合わせて求めます。
DGESVDの例題プログラムは 
の場合の特異値分解を示します。
入力データ
（本ルーチンの詳細はDGESDD のマニュアルページを参照）| このデータをダウンロード |
DGESDD Example Program Data 4 6 :Values of M and N 2.27 0.28 -0.48 1.07 -2.35 0.62 -1.54 -1.67 -3.09 1.22 2.93 -7.39 1.15 0.94 0.99 0.79 -1.45 1.03 -1.94 -0.78 -0.21 0.63 2.30 -2.57 :End of matrix A
出力結果
（本ルーチンの詳細はDGESDD のマニュアルページを参照）| この出力例をダウンロード |
DGESDD Example Program Results
Singular values
9.9966 3.6831 1.3569 0.5000
Left singular vectors
1 2 3 4
1 0.1921 0.8030 0.0041 0.5642
2 -0.8794 0.3926 -0.0752 -0.2587
3 0.2140 0.2980 0.7827 -0.5027
4 -0.3795 -0.3351 0.6178 0.6017
Right singular vectors by row (first m rows of V**T)
1 2 3 4 5 6
1 0.2774 0.2020 0.2918 -0.0938 -0.4213 0.7816
2 0.6003 0.0301 -0.3348 0.3699 -0.5266 -0.3353
3 -0.1277 0.2805 0.6453 0.6781 0.0413 -0.1645
4 -0.1323 -0.7034 -0.1906 0.5399 0.0575 0.3957
Error estimate for the singular values
2.2E-15
Error estimates for the left singular vectors
3.5E-16 9.5E-16 2.6E-15 2.6E-15
Error estimates for the right singular vectors
3.5E-16 9.5E-16 2.6E-15 4.4E-15
ソースコード
（本ルーチンの詳細はDGESDD のマニュアルページを参照）※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。
| このソースコードをダウンロード |
Program dgesdd_example
! DGESDD Example Program Text
! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com
! .. Use Statements ..
Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen
Use lapack_interfaces, Only: ddisna, dgesdd
Use lapack_precision, Only: dp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: nb = 64, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=dp) :: eps, serrbd
Integer :: i, ifail, info, lda, ldu, lwork, m, n
! .. Local Arrays ..
Real (Kind=dp), Allocatable :: a(:, :), rcondu(:), rcondv(:), s(:), &
u(:, :), uerrbd(:), verrbd(:), work(:)
Real (Kind=dp) :: dummy(1, 1)
Integer, Allocatable :: iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: epsilon, max, min, nint
! .. Executable Statements ..
Write (nout, *) 'DGESDD Example Program Results'
Write (nout, *)
! Skip heading in data file
Read (nin, *)
Read (nin, *) m, n
lda = m
ldu = m
Allocate (a(lda,n), rcondu(m), rcondv(m), s(m), u(ldu,m), uerrbd(m), &
verrbd(m), iwork(8*min(m,n)))
! Use routine workspace query to get optimal workspace.
lwork = -1
Call dgesdd('Overwrite A by tranpose(V)', m, n, a, lda, s, u, ldu, &
dummy, 1, dummy, lwork, iwork, info)
! Make sure that there is enough workspace for block size nb.
lwork = max((5*m+9)*m+n+nb*(m+n), nint(dummy(1,1)))
Allocate (work(lwork))
! Read the m by n matrix A from data file
Read (nin, *)(a(i,1:n), i=1, m)
! Compute the singular values and left and right singular vectors
! of A (A = U*S*(V**T), m.le.n)
Call dgesdd('Overwrite A by tranpose(V)', m, n, a, lda, s, u, ldu, &
dummy, 1, work, lwork, iwork, info)
If (info==0) Then
! Print solution
Write (nout, *) 'Singular values'
Write (nout, 100) s(1:m)
Flush (nout)
! Normalize so that u(1,j)>=0
Do i = 1, m
If (u(1,i)<0.0_dp) Then
u(1:m, i) = -u(1:m, i)
a(i, 1:n) = -a(i, 1:n)
End If
End Do
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call nagf_file_print_matrix_real_gen('General', ' ', m, m, u, ldu, &
'Left singular vectors', ifail)
Write (nout, *)
Flush (nout)
Call nagf_file_print_matrix_real_gen('General', ' ', m, n, a, lda, &
'Right singular vectors by row '//'(first m rows of V**T)', ifail)
! Get the machine precision, EPS and compute the approximate
! error bound for the computed singular values. Note that for
! the 2-norm, S(1) = norm(A)
eps = epsilon(1.0E0_dp)
serrbd = eps*s(1)
! Call DDISNA to estimate reciprocal condition
! numbers for the singular vectors
Call ddisna('Left', m, n, s, rcondu, info)
Call ddisna('Right', m, n, s, rcondv, info)
! Compute the error estimates for the singular vectors
Do i = 1, m
uerrbd(i) = serrbd/rcondu(i)
verrbd(i) = serrbd/rcondv(i)
End Do
! Print the approximate error bounds for the singular values
! and vectors
Write (nout, *)
Write (nout, *) 'Error estimate for the singular values'
Write (nout, 110) serrbd
Write (nout, *)
Write (nout, *) 'Error estimates for the left singular vectors'
Write (nout, 110) uerrbd(1:m)
Write (nout, *)
Write (nout, *) 'Error estimates for the right singular vectors'
Write (nout, 110) verrbd(1:m)
Else
Write (nout, 120) 'Failure in DGESDD. INFO =', info
End If
100 Format (3X, (8F8.4))
110 Format (4X, 1P, 6E11.1)
120 Format (1X, A, I4)
End Program