ZGESVD(複素特異値分解) LAPACKサンプルソースコード

ホーム > LAPACKサンプルプログラム目次 > 複素特異値分解 > 特異値分解

概要

本サンプルはFortran言語によりLAPACKルーチンZGESVDを利用するサンプルプログラムです。

6x4の行列の特異値と左および右特異ベクトルを求めます。

$\begin{displaymath} A = \left( \begin{array}{rrrrrr} 0.96 - 0.81i & -0.03 + 0... ...0 - 0.12i & -0.07 + 1.23i & 0.26 + 0.26i \end{array} \right), \end{displaymath}$

計算された特異値と特異ベクトルの誤差限界近似値も合わせて求めます。

ZGESDDの例題プログラムは $m \le n$ の場合の特異値分解を示します。

入力データ

（本ルーチンの詳細はZGESVD のマニュアルページを参照）

このデータをダウンロード

ZGESVD Example Program Data

   6             4                                       : m and n

 ( 0.96,-0.81) (-0.03, 0.96) (-0.91, 2.06) (-0.05, 0.41)
 (-0.98, 1.98) (-1.20, 0.19) (-0.66, 0.42) (-0.81, 0.56)
 ( 0.62,-0.46) ( 1.01, 0.02) ( 0.63,-0.17) (-1.11, 0.60)
 (-0.37, 0.38) ( 0.19,-0.54) (-0.98,-0.36) ( 0.22,-0.20)
 ( 0.83, 0.51) ( 0.20, 0.01) (-0.17,-0.46) ( 1.47, 1.59)
 ( 1.08,-0.28) ( 0.20,-0.12) (-0.07, 1.23) ( 0.26, 0.26) : Matrix A(1:m,1:n)

 ( 1.00, 0.00) ( 1.00, 0.00) ( 1.00, 0.00)
 ( 1.00, 0.00) ( 1.00, 0.00) ( 1.00, 0.00)               : RHS b(1:n)

出力結果

（本ルーチンの詳細はZGESVD のマニュアルページを参照）

この出力例をダウンロード

 ZGESVD Example Program Results

 Singular values of A:
         3.9994        3.0003        1.9944        0.9995

 Least squares solution:
    ( -0.0719, -0.2761)
    (  0.6796,  0.4326)
    ( -0.3832, -0.5447)
    (  0.0096, -0.3489)

 Norm of Residual:
         1.5266

ソースコード

（本ルーチンの詳細はZGESVD のマニュアルページを参照）

※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。

このソースコードをダウンロード

    Program zgesvd_example

!     ZGESVD Example Program Text

!     Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgesvd
      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 ..
      Complex (Kind=dp), Allocatable :: a(:, :), a_copy(:, :), b(:), u(:, :), &
        vt(:, :), work(:)
      Complex (Kind=dp) :: dummy(1, 1)
      Real (Kind=dp), Allocatable :: rwork(:), s(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: max, min, nint, real
!     .. Executable Statements ..
      Write (nout, *) 'ZGESVD 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), u(ldu,m), vt(ldvt,n), b(m), &
        rwork(5*n))

!     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 zgesvd('A', 'S', m, n, a, lda, s, u, ldu, vt, ldvt, dummy, lwork, &
        rwork, info)

!     Make sure that there is enough workspace for block size nb.
      lwork = max(m+3*n+nb*(m+n), nint(real(dummy(1,1))))
      Allocate (work(lwork))

!     Compute the singular values and left and right singular vectors
!     of A.

      Call zgesvd('A', 'S', m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, &
        rwork, info)

      If (info/=0) Then
        Write (nout, 100) 'Failure in ZGESVD. 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: dznrm2, zgemv
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In) :: lda, ldu, ldvt, m, n
!       .. Array Arguments ..
        Complex (Kind=dp), Intent (In) :: a(lda, n), u(ldu, m), vt(ldvt, n)
        Complex (Kind=dp), Intent (Inout) :: b(m)
        Real (Kind=dp), Intent (In) :: s(n)
!       .. Local Scalars ..
        Complex (Kind=dp) :: alpha, beta
        Real (Kind=dp) :: norm
!       .. Local Arrays ..
        Complex (Kind=dp), Allocatable :: x(:), y(:)
!       .. Intrinsic Procedures ..
        Intrinsic :: allocated, cmplx
!       .. Executable Statements ..
        Allocate (x(n), y(n))

!       Compute V*Inv(S)*U^H * b to get least squares solution.

!       y = U^T b
        alpha = cmplx(1.0_dp, 0.0_dp, kind=dp)
        beta = cmplx(0.0_dp, 0.0_dp, kind=dp)
        Call zgemv('C', m, n, alpha, u, ldu, b, 1, beta, y, 1)

        y(1:n) = y(1:n)/s(1:n)

!       x = V y
        Call zgemv('C', 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 = cmplx(-1.0_dp, 0.0_dp, kind=dp)
        beta = cmplx(1._dp, 0.0_dp, kind=dp)
        Call zgemv('N', m, n, alpha, a, lda, x, 1, beta, b, 1)

        norm = dznrm2(m, b, 1)

        Write (nout, *)
        Write (nout, *) 'Norm of Residual:'
        Write (nout, 110) norm

        If (allocated(x)) Then
          Deallocate (x)
        End If
        If (allocated(y)) Then
          Deallocate (y)
        End If

100     Format (4X, '(', F8.4, ',', F8.4, ')')
110     Format (4X, 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

ご案内

複素特異値分解: 特異値分解 : (オプショナルで左および右特異ベクトル)

LAPACKサンプルソースコード : 使用ルーチン名：ZGESVD

概要

入力データ

出力結果

ソースコード