計算ルーチン: 複素三角-五角行列の QR 分解

ZTPQRT Example Program Data

  6      4      2                                        : m, n and nrhs

 ( 0.96,-0.81) (-0.03, 0.96) (-0.91, 2.06) (-0.05, 0.41)
 (-0.98, 1.98) (-1.20, 0.19) (-0.66, 0.42) (-0.81, 0.56)
 ( 0.62,-0.46) ( 1.01, 0.02) ( 0.63,-0.17) (-1.11, 0.60)
 (-0.37, 0.38) ( 0.19,-0.54) (-0.98,-0.36) ( 0.22,-0.20)
 ( 0.83, 0.51) ( 0.20, 0.01) (-0.17,-0.46) ( 1.47, 1.59)
 ( 1.08,-0.28) ( 0.20,-0.12) (-0.07, 1.23) ( 0.26, 0.26) : matrix A

 (-2.09, 1.93) ( 3.26,-2.70)
 ( 3.34,-3.53) (-6.22, 1.16)
 (-4.94,-2.04) ( 7.94,-3.13)
 ( 0.17, 4.23) ( 1.04,-4.26)
 (-5.19, 3.63) (-2.31,-2.12)
 ( 0.98, 2.53) (-1.39,-4.05)                             : matrix B

 ZTPQRT Example Program Results

 solution(s) for n rows
                    1                 2
 1  (-0.5091,-1.2428) ( 0.7569, 1.4384)
 2  (-2.3789, 2.8651) ( 5.1727,-3.6193)
 3  ( 1.4634,-2.2064) (-2.6613, 2.1339)
 4  ( 0.4701, 2.6964) (-2.6933, 0.2724)

 Least squares solution(s) for all rows
                    1                 2
 1  (-0.5044,-1.2179) ( 0.7629, 1.4529)
 2  (-2.4281, 2.8574) ( 5.1570,-3.6089)
 3  ( 1.4872,-2.1955) (-2.6518, 2.1203)
 4  ( 0.4537, 2.6904) (-2.7606, 0.3318)

 Square root(s) of the residual sum(s) of squares
        6.88E-02   1.87E-01

    Program ztpqrt_example

!     ZTPQRT Example Program Text

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

!     .. Use Statements ..
      Use blas_interfaces, Only: dznrm2
      Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen_comp
      Use lapack_interfaces, Only: zgemqrt, zgeqrt, ztpmqrt, ztpqrt, ztrtrs
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nbmax = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, ifail, info, j, lda, ldb, ldt, lwork, m, n, nb, nrhs
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:, :), c(:, :), t(:, :), &
        work(:)
      Real (Kind=dp), Allocatable :: rnorm(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: max, min
!     .. Executable Statements ..
      Write (nout, *) 'ZTPQRT Example Program Results'
      Write (nout, *)
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, nrhs
      lda = m
      ldb = m
      nb = min(m, n, nbmax)
      ldt = nb
      lwork = nb*max(n, m)
      Allocate (a(lda,n), b(ldb,nrhs), c(ldb,nrhs), rnorm(nrhs), t(ldt,min(m, &
        n)), work(lwork))

!     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)

      c(1:m, 1:nrhs) = b(1:m, 1:nrhs)
!     Compute the QR factorization of first n rows of A
      Call zgeqrt(n, n, nb, a, lda, t, ldt, work, info)

!     Compute C = (C1) = (Q**H)*B, storing the result in C
!                 (C2)
      Call zgemqrt('Left', 'Conjugate Transpose', n, nrhs, n, nb, a, lda, t, &
        ldt, c, ldb, work, info)

      b(1:n, 1:nrhs) = c(1:n, 1:nrhs)
!     Compute least squares solutions for first n rows by back-substitution in
!     R*X = C1
      Call ztrtrs('Upper', 'No transpose', 'Non-Unit', n, nrhs, a, lda, c, &
        ldb, 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

!       Print solution using first n rows

!       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, &
          c, ldb, 'Bracketed', 'F7.4', 'solution(s) for n rows', 'Integer', &
          rlabs, 'Integer', clabs, 80, 0, ifail)

      End If

!     Now add the remaining rows and perform QR update
      Call ztpqrt(m-n, n, 0, nb, a, lda, a(n+1,1), lda, t, ldt, work, info)

!     Apply orthogonal transformations to C
      Call ztpmqrt('Left', 'Conjugate Transpose', m-n, nrhs, n, 0, nb, &
        a(n+1,1), lda, t, ldt, b, ldb, b(5,1), ldb, work, info)
!     Compute least squares solutions for first n rows by back-substitution in
!     R*X = C1
      Call ztrtrs('Upper', 'No transpose', 'Non-Unit', n, nrhs, a, lda, b, &
        ldb, 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

!       Print least squares solutions
        Write (nout, *)
        Flush (nout)
        ifail = 0
        Call nagf_file_print_matrix_complex_gen_comp('G', ' ', n, nrhs, b, &
          ldb, 'Bracketed', 'F7.4', 'Least squares solution(s) for all rows', &
          'Integer', rlabs, 'Integer', clabs, 80, 0, ifail)

!       Compute and print estimates of the square roots of the residual
!       sums of squares

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

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

100   Format (5X, 1P, 7E11.2)
    End Program