計算ルーチン: BLAS-3 を用いた複素一般矩形行列の QR 分解 : (列ピボット)

ZGEQP3 Example Program Data

   5             4             2                 :Values of M, N and NRHS

 ( 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.22, 2.35)
 (-2.61,-1.49) ( 1.62,-1.48)
 ( 3.13,-3.61) ( 1.65, 3.43)
 ( 7.33,-8.01) (-0.98, 3.08)
 ( 9.12, 7.63) (-2.84, 2.78)                             :End of matrix B

 ZGEQP3 Example Program Results

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

 Least squares solution(s)
                    1                 2
 1  ( 0.0000, 0.0000) ( 0.0000, 0.0000)
 2  ( 2.7020, 8.0911) (-2.2682,-2.9884)
 3  ( 2.8888, 2.5012) ( 0.9779, 1.3565)
 4  ( 2.7100, 0.4791) (-1.3734, 0.2212)

 Square root(s) of the residual sum(s) of squares
      2.51E-01   8.10E-02

    Program zgeqp3_example

!     ZGEQP3 Example Program Text

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

!     .. Use Statements ..
      Use blas_interfaces, Only: dznrm2, ztrsm
      Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen_comp
      Use lapack_interfaces, Only: zgeqp3, zunmqr
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=dp), Parameter :: one = (1.0E0_dp, 0.0E0_dp)
      Complex (Kind=dp), Parameter :: zero = (0.0E0_dp, 0.0E0_dp)
      Integer, Parameter :: inc1 = 1, nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: tol
      Integer :: i, ifail, info, j, k, lda, ldb, lwork, m, n, nrhs
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:, :), tau(:), work(:)
      Real (Kind=dp), Allocatable :: rnorm(:), rwork(:)
      Integer, Allocatable :: jpvt(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs
!     .. Executable Statements ..
      Write (nout, *) 'ZGEQP3 Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, nrhs
      lda = m
      ldb = m
      lwork = (n+1)*nb
      Allocate (a(lda,n), b(ldb,nrhs), tau(n), work(lwork), rnorm(nrhs), &
        rwork(2*n), jpvt(n))

!     Read A and B from data file

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

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

      jpvt(1:n) = 0

!     Compute the QR factorization of A

      Call zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)

!     Compute C = (C1) = (Q**H)*B, storing the result in B
!                  (C2)

      Call zunmqr('Left', 'Conjugate Transpose', m, nrhs, n, a, lda, tau, b, &
        ldb, work, lwork, info)

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

      tol = 0.01_dp

!     Determine and print the rank, K, of R relative to TOL

loop: Do k = 1, n
        If (abs(a(k,k))<=tol*abs(a(1,1))) Then
          Exit loop
        End If
      End Do loop
      k = k - 1

      Write (nout, *) 'Tolerance used to estimate the rank of A'
      Write (nout, 100) tol
      Write (nout, *) 'Estimated rank of A'
      Write (nout, 110) k
      Write (nout, *)
      Flush (nout)

!     Compute least squares solutions by back-substitution in
!     R(1:K,1:K)*Y = C1, storing the result in B

      Call ztrsm('Left', 'Upper', 'No transpose', 'Non-Unit', k, nrhs, one, a, &
        lda, b, ldb)

!     Compute estimates of the square roots of the residual sums of
!     squares (2-norm of each of the columns of C2)

      Do j = 1, nrhs
        rnorm(j) = dznrm2(m-k, b(k+1,j), inc1)
      End Do

!     Set the remaining elements of the solutions to zero (to give
!     the basic solutions)

      b(k+1:n, 1:nrhs) = zero

!     Permute the least squares solutions stored in B to give X = P*Y

      Do j = 1, nrhs
        work(jpvt(1:n)) = b(1:n, j)
        b(1:n, j) = work(1:n)
      End Do

!     Print least squares solutions

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, nrhs, b, &
        ldb, 'Bracketed', 'F7.4', 'Least squares solution(s)', 'Integer', &
        rlabs, 'Integer', clabs, 80, 0, ifail)

!     Print the square roots of the residual sums of squares

      Write (nout, *)
      Write (nout, *) 'Square root(s) of the residual sum(s) of squares'
      Write (nout, 100) rnorm(1:nrhs)

100   Format (3X, 1P, 7E11.2)
110   Format (1X, I8)
    End Program