計算ルーチン: 複素行列の特異値分解 : (オプションで左および/または右特異ベクトル(高速ヤコビ))

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

概要

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

入力データ

(本ルーチンの詳細はZGESVJ のマニュアルページを参照)

このデータをダウンロード
ZGESVJ 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

出力結果

(本ルーチンの詳細はZGESVJ のマニュアルページを参照)

この出力例をダウンロード
 ZGESVJ Example Program Results

 Singular values
     3.9994  3.0003  1.9944  0.9995

 Left singular vectors
          1       2       3       4
 1   0.5632 -0.1814 -0.4517 -0.0171
     0.0151 -0.3390  0.2700  0.4820
 
 2  -0.1386 -0.3097  0.1991  0.0245
     0.6070  0.0231  0.4227 -0.0405
 
 3   0.0864 -0.2677 -0.4069  0.1363
    -0.1588  0.3268 -0.4452 -0.1991
 
 4  -0.1486  0.2390 -0.1554  0.5883
     0.1489 -0.1115  0.0051 -0.1077
 
 5   0.2459  0.5084  0.0337  0.0479
     0.1000  0.4916  0.2629  0.0891
 
 6   0.3780 -0.0068 -0.1676  0.0806
    -0.0681 -0.0870  0.1359 -0.5766

 Right singular vectors
          1       2       3       4
 1   0.6968  0.2307 -0.0278 -0.2602
     0.0207  0.0670 -0.5116 -0.3552
 
 2   0.0972  0.0045 -0.1897  0.8260
    -0.3520  0.2446 -0.2932  0.0511
 
 3  -0.0399 -0.3854  0.5192 -0.2068
    -0.5414  0.4391  0.0397 -0.2237
 
 4   0.1944  0.7355  0.4098  0.1193
    -0.2196 -0.0131  0.4201  0.1192

 Error estimates (as multiples of machine precision):

   for the singular values
       4

   for left singular vectors
       4   4   4   4

   for right singular vectors
       4   4   4   4

ソースコード

(本ルーチンの詳細はZGESVJ のマニュアルページを参照)

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


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

!     ZGESVJ Example Program Text

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

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen
      Use lapack_interfaces, Only: ddisna, zgesvj
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: eps, serrbd
      Integer :: i, ifail, info, j, lda, ldv, lrwork, lwork, m, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), cwork(:), v(:, :)
      Real (Kind=dp), Allocatable :: rcondu(:), rcondv(:), rwork(:), sva(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, epsilon, max, nint
!     .. Executable Statements ..
      Write (nout, *) 'ZGESVJ Example Program Results'
      Write (nout, *)
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      ldv = n
      lwork = n + m
      lrwork = max(6, n)
      Allocate (a(lda,n), rcondu(m), rcondv(m), sva(n), v(ldv,n), &
        cwork(lwork), rwork(lrwork))

!     Read the m by n matrix A from data file

      Read (nin, *)((a(i,j),j=1,n), i=1, m)

!     Compute the singular values and left and right singular vectors
!     of A (A = U*S*V, m.ge.n)

      Call zgesvj('G', 'U', 'V', m, n, a, lda, sva, 0, v, ldv, cwork, lwork, &
        rwork, lrwork, info)
      If (info==0) Then

!       Compute the approximate error bound for the computed singular values
!       using the 2-norm, s(1) = norm(A), and machine precision, eps.
        eps = epsilon(1.0E0_dp)
        serrbd = eps*sva(1)

!       Print solution
        Write (nout, *) 'Singular values'
        Write (nout, 100)(sva(j), j=1, n)

        If (abs(rwork(1)-1.0_dp)>eps) Then
          Write (nout, 130) 'Values need scaling by factor = ', rwork(1)
        End If
        Write (nout, *)
        Flush (nout)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_file_print_matrix_complex_gen('General', ' ', m, n, a, lda, &
          'Left singular vectors', ifail)

        Write (nout, *)
        Flush (nout)

        ifail = 0
        Call nagf_file_print_matrix_complex_gen('General', ' ', n, n, v, ldv, &
          'Right singular vectors', ifail)

!       Call DDISNA to estimate reciprocal condition
!       numbers for the singular vectors
        Call ddisna('Left', m, n, sva, rcondu, info)
        Call ddisna('Right', m, n, sva, rcondv, info)

!       Print the approximate error bounds for the singular values
!       and vectors
        Write (nout, '(/1X,A)') &
          'Error estimates (as multiples of machine precision):'
        Write (nout, '(/1X,A)') '  for the singular values'
        Write (nout, 110) nint(serrbd/epsilon(1.0E0_dp))
        Write (nout, '(/1X,A)') '  for left singular vectors'
        Write (nout, 110)(nint(serrbd/rcondu(i)/epsilon(1.0E0_dp)), i=1, n)
        Write (nout, '(/1X,A)') '  for right singular vectors'
        Write (nout, 110)(nint(serrbd/rcondv(i)/epsilon(1.0E0_dp)), i=1, n)
      Else
        Write (nout, 120) 'Failure in ZGESVJ. INFO =', info
      End If

100   Format (3X, (8F8.4))
110   Format (4X, 6I4)
120   Format (1X, A, I4)
130   Format (/, 1X, A, 1P, E13.5)
    End Program


ご案内
関連情報
Privacy Policy  /  Trademarks