ZGELSD(複素線形最小二乗) LAPACKサンプルソースコード

概要

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

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

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

最小ノルム解を

とします。そして

$\begin{displaymath} A = \left( \begin{array}{rrrrrr} 0.47 - 0.34 i & -0.32 - ... ...- 0.09 i & 1.14 - 0.78 i & 0.07 + 1.14 i \end{array} \right) \end{displaymath}$

及び

$\begin{displaymath} b = \left( \begin{array}{r} 2.15 - 0.20 i \\ -2.24 + 1.82 i \\ 4.45 - 4.28 i \\ 5.70 - 6.25 i \end{array} \right). \end{displaymath}$

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

入力データ

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

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

ZGELSD Example Program Data

   4             5                                            :Values of M and N

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

 ( 2.15,-0.20)
 (-2.24, 1.82)
 ( 4.45,-4.28)
 ( 5.70,-6.25)                                                  :End of vector b

出力結果

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

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

Warning: Floating underflow occurred
 ZGELSD Example Program Results

 Least squares solution
 ( 3.9747,-1.8377) (-0.9186, 0.8253) (-0.3105, 0.1477) ( 1.0050, 0.8626)
 (-0.2256,-1.9425)

 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

ソースコード

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

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

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

    Program zgelsd_example

!     ZGELSD Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgelsd
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: rcond
      Integer :: i, info, lda, liwork, lrwork, lwork, m, n, rank
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:), work(:)
      Complex (Kind=dp) :: lw(1)
      Real (Kind=dp) :: lrw(1)
      Real (Kind=dp), Allocatable :: rwork(:), s(:)
      Integer, Allocatable :: iwork(:)
      Integer :: liw(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: nint, real
!     .. Executable Statements ..
      Write (nout, *) 'ZGELSD Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      Allocate (a(lda,n), b(n), s(m))

!     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

!     Call zgelsd in workspace query mode.
      lwork = -1
      Call zgelsd(m, n, 1, a, lda, b, n, s, rcond, rank, lw, lwork, lrw, liw, &
        info)
      lwork = nint(real(lw(1)))
      lrwork = nint(lrw(1))
      liwork = liw(1)
      Allocate (work(lwork), rwork(lrwork), iwork(liwork))

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

      Call zgelsd(m, n, 1, a, lda, b, n, s, rcond, rank, work, lwork, rwork, &
        iwork, 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:m)
      Else If (info>0) Then
        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

ご案内

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

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

概要

入力データ

出力結果

ソースコード