概要
本サンプルはFortran言語によりLAPACKルーチンDGESVDを利用するサンプルプログラムです。
6x4の行列の特異値と左および右特異ベクトルを求めます。
計算された特異値と特異ベクトルの誤差限界近似値も合わせて求めます。
DGESDDの例題プログラムは の場合の特異値分解を示します。
入力データ
(本ルーチンの詳細はDGESVD のマニュアルページを参照)このデータをダウンロード |
DGESVD Example Program Data 6 4 :Values of M and N 2.27 -1.54 1.15 -1.94 0.28 -1.67 0.94 -0.78 -0.48 -3.09 0.99 -0.21 1.07 1.22 0.79 0.63 -2.35 2.93 -1.45 2.30 0.62 -7.39 1.03 -2.57 :End of matrix A 1.0 1.0 1.0 1.0 1.0 1.0 :RHS b(1:m)
出力結果
(本ルーチンの詳細はDGESVD のマニュアルページを参照)この出力例をダウンロード |
DGESVD Example Program Results Singular values of A: 9.9966 3.6831 1.3569 0.5000 Least squares solution: -0.0563 -0.1700 0.8202 0.5545 Norm of Residual: 1.7472
ソースコード
(本ルーチンの詳細はDGESVD のマニュアルページを参照)※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。
このソースコードをダウンロード |
Program dgesvd_example ! DGESVD Example Program Text ! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com ! .. Use Statements .. Use lapack_interfaces, Only: dgesvd Use lapack_precision, Only: dp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nb = 64, nin = 5, nout = 6, prerr = 0 ! .. Local Scalars .. Integer :: i, info, lda, ldu, ldvt, lwork, m, n ! .. Local Arrays .. Real (Kind=dp), Allocatable :: a(:, :), a_copy(:, :), b(:), s(:), & u(:, :), vt(:, :), work(:) Real (Kind=dp) :: dummy(1, 1) ! .. Intrinsic Procedures .. Intrinsic :: max, min, nint ! .. Executable Statements .. Write (nout, *) 'DGESVD Example Program Results' Write (nout, *) ! Skip heading in data file Read (nin, *) Read (nin, *) m, n lda = m ldu = m ldvt = n Allocate (a(lda,n), a_copy(m,n), s(n), vt(ldvt,n), u(ldu,m), b(m)) ! Read the m by n matrix A from data file Read (nin, *)(a(i,1:n), i=1, m) ! Read the right hand side of the linear system Read (nin, *) b(1:m) a_copy(1:m, 1:n) = a(1:m, 1:n) ! Use routine workspace query to get optimal workspace. lwork = -1 Call dgesvd('A', 'S', m, n, a, lda, s, u, ldu, vt, ldvt, dummy, lwork, & info) ! Make sure that there is enough workspace for block size nb. lwork = max(m+4*n+nb*(m+n), nint(dummy(1,1))) Allocate (work(lwork)) ! Compute the singular values and left and right singular vectors ! of A. Call dgesvd('A', 'S', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, & info) If (info/=0) Then Write (nout, 100) 'Failure in DGESVD. INFO =', info 100 Format (1X, A, I4) Go To 120 End If ! Print the significant singular values of A Write (nout, *) 'Singular values of A:' Write (nout, 110) s(1:min(m,n)) 110 Format (1X, 4(3X,F11.4)) If (prerr>0) Then Call compute_error_bounds(m, n, s) End If If (m>n) Then ! Compute V*Inv(S)*U^T * b to get least squares solution. Call compute_least_squares(m, n, a_copy, m, u, ldu, vt, ldvt, s, b) End If 120 Continue Contains Subroutine compute_least_squares(m, n, a, lda, u, ldu, vt, ldvt, s, b) ! .. Use Statements .. Use blas_interfaces, Only: dgemv, dnrm2 ! .. Implicit None Statement .. Implicit None ! .. Scalar Arguments .. Integer, Intent (In) :: lda, ldu, ldvt, m, n ! .. Array Arguments .. Real (Kind=dp), Intent (In) :: a(lda, n), s(n), u(ldu, m), vt(ldvt, n) Real (Kind=dp), Intent (Inout) :: b(m) ! .. Local Scalars .. Real (Kind=dp) :: alpha, beta, norm ! .. Local Arrays .. Real (Kind=dp), Allocatable :: x(:), y(:) ! .. Intrinsic Procedures .. Intrinsic :: allocated ! .. Executable Statements .. Allocate (x(n), y(n)) ! Compute V*Inv(S)*U^T * b to get least squares solution. ! y = U^T b alpha = 1._dp beta = 0._dp Call dgemv('T', m, n, alpha, u, ldu, b, 1, beta, y, 1) y(1:n) = y(1:n)/s(1:n) ! x = V y Call dgemv('T', n, n, alpha, vt, ldvt, y, 1, beta, x, 1) Write (nout, *) Write (nout, *) 'Least squares solution:' Write (nout, 100) x(1:n) ! Find norm of residual ||b-Ax||. alpha = -1._dp beta = 1._dp Call dgemv('N', m, n, alpha, a, lda, x, 1, beta, b, 1) norm = dnrm2(m, b, 1) Write (nout, *) Write (nout, *) 'Norm of Residual:' Write (nout, 100) norm If (allocated(x)) Then Deallocate (x) End If If (allocated(y)) Then Deallocate (y) End If 100 Format (1X, 4(3X,F11.4)) End Subroutine Subroutine compute_error_bounds(m, n, s) ! Error estimates for singular values and vectors is computed ! and printed here. ! .. Use Statements .. Use lapack_interfaces, Only: ddisna Use lapack_precision, Only: dp ! .. Implicit None Statement .. Implicit None ! .. Scalar Arguments .. Integer, Intent (In) :: m, n ! .. Array Arguments .. Real (Kind=dp), Intent (In) :: s(n) ! .. Local Scalars .. Real (Kind=dp) :: eps, serrbd Integer :: i, info ! .. Local Arrays .. Real (Kind=dp), Allocatable :: rcondu(:), rcondv(:), uerrbd(:), & verrbd(:) ! .. Intrinsic Procedures .. Intrinsic :: epsilon ! .. Executable Statements .. Allocate (rcondu(n), rcondv(n), uerrbd(n), verrbd(n)) ! 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, n 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, 100) serrbd Write (nout, *) Write (nout, *) 'Error estimates for the left singular vectors' Write (nout, 100) uerrbd(1:n) Write (nout, *) Write (nout, *) 'Error estimates for the right singular vectors' Write (nout, 100) verrbd(1:n) 100 Format (4X, 1P, 6E11.1) End Subroutine End Program