複素線形最小二乗: 複素数 : (特異値分解)

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

ホーム > LAPACKサンプルプログラム目次 > 複素線形最小二乗 > 複素数

概要

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

以下の線形最小二乗問題を解きます。

\begin{displaymath} \min_{x} \left\Vert b - A x \right\Vert _{2} \end{displaymath}

最小ノルム解を $ x$とします。そして

\begin{displaymath} A = \left( \begin{array}{rrrr} 0.47 - 0.34 i & -0.40 + 0.... ... - 0.85 i & 1.44 + 0.80 i & 0.07 + 1.14 i \end{array} \right) \end{displaymath}

及び

\begin{displaymath} b = \left( \begin{array}{r} -1.08 - 2.59 i \\ -2.61 - 1... ... i \\ 7.33 - 8.01 i \\ 9.12 + 7.63 i \end{array} \right). \end{displaymath}

$ A$の有効階数を求めるために誤差0.01が使われています。

入力データ

（本ルーチンの詳細はZGELSS のマニュアルページを参照）
このデータをダウンロード 
ZGELSS Example Program Data

   5             4                                       :Values of M and N

 ( 0.47,-0.34) (-0.40, 0.54) ( 0.60, 0.01) ( 0.80,-1.02)
 (-0.32,-0.23) (-0.05, 0.20) (-0.26,-0.44) (-0.43, 0.17)
 ( 0.35,-0.60) (-0.52,-0.34) ( 0.87,-0.11) (-0.34,-0.09)
 ( 0.89, 0.71) (-0.45,-0.45) (-0.02,-0.57) ( 1.14,-0.78)
 (-0.19, 0.06) ( 0.11,-0.85) ( 1.44, 0.80) ( 0.07, 1.14) :End of matrix A

 (-1.08,-2.59)
 (-2.61,-1.49)
 ( 3.13,-3.61)
 ( 7.33,-8.01)
 ( 9.12, 7.63)                                           :End of vector b

出力結果

（本ルーチンの詳細はZGELSS のマニュアルページを参照）
この出力例をダウンロード 
 ZGELSS Example Program Results

 Least squares solution
 ( 1.1673,-3.3222) ( 1.3480, 5.5028) ( 4.1762, 2.3434) ( 0.6465, 0.0105)

 Tolerance used to estimate the rank of A
      1.00E-02
 Estimated rank of A
      3

 Singular values of A
      2.9979     1.9983     1.0044     0.0064

ソースコード

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

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

!     ZGELSS Example Program Text

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

!     .. Use Statements ..
      Use blas_interfaces, Only: dznrm2
      Use lapack_interfaces, Only: zgelss
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: rcond, rnorm
      Integer :: i, info, lda, lwork, m, n, rank
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:), work(:)
      Real (Kind=dp), Allocatable :: rwork(:), s(:)
!     .. Executable Statements ..
      Write (nout, *) 'ZGELSS Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      lwork = 2*n + nb*(m+n)
      Allocate (a(lda,n), b(m), work(lwork), rwork(5*n), s(n))

!     Read A and B from data file

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

!     Choose RCOND to reflect the relative accuracy of the input data

      rcond = 0.01E0_dp

!     Solve the least squares problem min( norm2(b - Ax) ) for the x
!     of minimum norm.

      Call zgelss(m, n, 1, a, lda, b, m, s, rcond, rank, work, lwork, rwork, &
        info)

      If (info==0) Then

!       Print solution

        Write (nout, *) 'Least squares solution'
        Write (nout, 100) b(1:n)

!       Print the effective rank of A

        Write (nout, *)
        Write (nout, *) 'Tolerance used to estimate the rank of A'
        Write (nout, 110) rcond
        Write (nout, *) 'Estimated rank of A'
        Write (nout, 120) rank

!       Print singular values of A

        Write (nout, *)
        Write (nout, *) 'Singular values of A'
        Write (nout, 130) s(1:n)

!       Compute and print estimate of the square root of the
!       residual sum of squares

        If (rank==n) Then
          rnorm = dznrm2(m-n, b(n+1), 1)
          Write (nout, *)
          Write (nout, *) 'Square root of the residual sum of squares'
          Write (nout, 110) rnorm
        End If
      Else
        Write (nout, *) 'The SVD algorithm failed to converge'
      End If

100   Format (4(' (',F7.4,',',F7.4,')',:))
110   Format (3X, 1P, E11.2)
120   Format (1X, I6)
130   Format (1X, 7F11.4)
    End Program


ご案内
関連情報
Privacy Policy  /  Trademarks