計算ルーチン: 複素一般矩形行列の RQ 分解

ZGERQF Example Program Data

   3             4                               :Values of M, N and NRHS

 ( 0.28,-0.36) ( 0.50,-0.86) (-0.77,-0.48) ( 1.58, 0.66)
 (-0.50,-1.10) (-1.21, 0.76) (-0.32,-0.24) (-0.27,-1.15)
 ( 0.36,-0.51) (-0.07, 1.33) (-0.75, 0.47) (-0.08, 1.01) :End of matrix A

 (-1.35, 0.19)
 ( 9.41,-3.56)
 (-7.57, 6.93)                                           :End of vector b

 ZGERQF Example Program Results

 Minimum-norm solution
 ( -2.8501,  6.4683) (  1.6264, -0.7799) (  6.9290,  4.6481) (  1.4048,  3.2400)

    Program zgerqf_example

!     ZGERQF Example Program Text

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

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgerqf, ztrtrs, zunmrq
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=dp), Parameter :: zero = (0.0_dp, 0.0_dp)
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, info, lda, lwork, m, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:), tau(:), work(:), x(:)
!     .. Executable Statements ..
      Write (nout, *) 'ZGERQF Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      lwork = nb*m
      Allocate (a(lda,n), b(m), tau(m), work(lwork), x(n))

!     Read the matrix A and the vector b from data file

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

!     Compute the RQ factorization of A
      Call zgerqf(m, n, a, lda, tau, work, lwork, info)

!     Copy the m-element vector b into elements x(n-m+1), ..., x(n) of x

      x(n-m+1:n) = b(1:m)

!     Solve R*y2 = b, storing the result in x2
      Call ztrtrs('Upper', 'No transpose', 'Non-Unit', m, 1, a(1,n-m+1), lda, &
        x(n-m+1), m, info)

      If (info>0) Then
        Write (nout, *) 'The upper triangular factor, R, of A is singular, '
        Write (nout, *) 'the least squares solution could not be computed'
      Else

!       Set y1 to zero (stored in x(1:n-m))

        x(1:n-m) = zero

!       Compute minimum-norm solution x = (Q**H)*y
        Call zunmrq('Left', 'Conjugate transpose', n, 1, m, a, lda, tau, x, n, &
          work, lwork, info)

!       Print minimum-norm solution

        Write (nout, *) 'Minimum-norm solution'
        Write (nout, 100) x(1:n)
      End If

100   Format (4(' (',F8.4,',',F8.4,')',:))
    End Program