実線形最小二乗: 直行分解

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

ホーム > LAPACKサンプルプログラム目次 > 実線形最小二乗 > 直行分解

概要

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

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

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

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

\begin{displaymath} A = \left( \begin{array}{rrrrrr} -0.09 & -1.56 & -1.48 & ... ...7.4 \\ 4.3 \\ -8.1 \\ 1.8 \\ 8.7 \end{array} \right). \end{displaymath}

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

入力データ

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

  5      6                                :Values of M and N

 -0.09  -1.56  -1.48  -1.09   0.08  -1.59
  0.14   0.20  -0.43   0.84   0.55  -0.72
 -0.46   0.29   0.89   0.77  -1.13   1.06
  0.68   1.09  -0.71   2.11   0.14   1.24
  1.29   0.51  -0.96  -1.27   1.74   0.34 :End of matrix A

  7.4
  4.3
 -8.1
  1.8
  8.7                                     :End of vector b

出力結果

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

 Least squares solution
      1.5938    -0.1180    -3.1501     0.1554     2.5529    -1.6730

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

 Singular values of A
      3.9997     2.9962     2.0001     0.9988     0.0025

ソースコード

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

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

!     DGELSD Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: dgelsd
      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, lwork, m, n, rank
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), b(:), s(:), work(:)
      Real (Kind=dp) :: lw(1)
      Integer, Allocatable :: iwork(:)
      Integer :: liw(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: nint
!     .. Executable Statements ..
      Write (nout, *) 'DGELSD 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.01_dp

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

!     Now Solve the least squares problem min( norm2(b - Ax) ) for the
!     x of minimum norm.
      Call dgelsd(m, n, 1, a, lda, b, n, s, rcond, rank, work, lwork, 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, 100) s(1:m)
      Else
        Write (nout, *) 'The SVD algorithm failed to converge'
      End If

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


ご案内
関連情報
Privacy Policy  /  Trademarks