複素線形最小二乗: 直行分解

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

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

概要

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

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

\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が使われています。

入力データ

（本ルーチンの詳細はZGELSY のマニュアルページを参照）
このデータをダウンロード 
ZGELSY 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

出力結果

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

 Least squares solution
 ( 1.1669,-3.3224) ( 1.3486, 5.5027) ( 4.1764, 2.3435) ( 0.6467, 0.0107)

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

ソースコード

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

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

!     ZGELSY Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgelsy
      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
      Integer :: i, info, lda, lwork, m, n, rank
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:), work(:)
      Real (Kind=dp), Allocatable :: rwork(:)
      Integer, Allocatable :: jpvt(:)
!     .. Executable Statements ..
      Write (nout, *) 'ZGELSY Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      lwork = nb*(n+1)
      Allocate (a(lda,n), b(m), work(lwork), rwork(2*n), jpvt(n))

!     Read A and B from data file

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

!     Initialize JPVT to be zero so that all columns are free

      jpvt(1:n) = 0

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

      rcond = 0.01_dp

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

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

!     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

100   Format (4(' (',F7.4,',',F7.4,')',:))
110   Format (1X, 1P, E10.2)
120   Format (1X, I6)
    End Program


ご案内
関連情報
Privacy Policy  /  Trademarks