計算ルーチン: ZGEBRD により決まる準対角形への縮約からのユニタリ変換の適用

ZUNMBR Example Program Data
  6  4                                        :Values of M and N, Example 1
 ( 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)   :End of matrix A
  3  4                                        :Values of M and N, Example 2
 ( 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

 ZUNMBR Example Program Results

 Example 1: bidiagonal matrix B
 Diagonal
    -3.0870 -2.0660 -1.8731 -2.0022
 Superdiagonal
     2.1126 -1.2628  1.6126

 Example 1: matrix Q
                    1                 2                 3                 4
 1  (-0.3110, 0.2624) ( 0.6521, 0.5532) ( 0.0427, 0.0361) (-0.2634,-0.0741)
 2  ( 0.3175,-0.6414) ( 0.3488, 0.0721) ( 0.2287, 0.0069) ( 0.1101,-0.0326)
 3  (-0.2008, 0.1490) (-0.3103, 0.0230) ( 0.1855,-0.1817) (-0.2956, 0.5648)
 4  ( 0.1199,-0.1231) (-0.0046,-0.0005) (-0.3305, 0.4821) (-0.0675, 0.3464)
 5  (-0.2689,-0.1652) ( 0.1794,-0.0586) (-0.5235,-0.2580) ( 0.3927, 0.1450)
 6  (-0.3499, 0.0907) ( 0.0829,-0.0506) ( 0.3202, 0.3038) ( 0.3174, 0.3241)

 Example 2: bidiagonal matrix B
 Diagonal
     2.7615  1.6298 -1.3275
 Superdiagonal
    -0.9500 -1.0183

 Example 2: matrix P**H
                    1                 2                 3                 4
 1  (-0.1258, 0.1618) (-0.2247, 0.3864) ( 0.3460, 0.2157) (-0.7099,-0.2966)
 2  ( 0.4148, 0.1795) ( 0.1368,-0.3976) ( 0.6885, 0.3386) ( 0.1667,-0.0494)
 3  ( 0.4575,-0.4807) (-0.2733, 0.4981) (-0.0230, 0.3861) ( 0.1730, 0.2395)

    Program zunmbr_example

!     ZUNMBR Example Program Text

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

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_blas_zmload, &
        nagf_file_print_matrix_complex_gen_comp
      Use lapack_interfaces, Only: zgebrd, zgelqf, zgeqrf, zlacpy, zunglq, &
        zungqr, zunmbr
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=dp), Parameter :: zero = (0.0E0_dp, 0.0E0_dp)
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, ic, ifail, info, lda, ldph, ldu, lwork, m, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), ph(:, :), tau(:), taup(:), &
        tauq(:), u(:, :), work(:)
      Real (Kind=dp), Allocatable :: d(:), e(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Executable Statements ..
      Write (nout, *) 'ZUNMBR Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Do ic = 1, 2
        Read (nin, *) m, n
        lda = m
        ldph = n
        ldu = m
        lwork = 64*(m+n)
        Allocate (a(lda,n), ph(ldph,n), tau(n), taup(n), tauq(n), u(ldu,n), &
          work(lwork), d(n), e(n-1))

!       Read A from data file

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

        If (m>=n) Then

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

!         Copy A to U
          Call zlacpy('Lower', m, n, a, lda, u, ldu)

!         Form Q explicitly, storing the result in U
          Call zungqr(m, n, n, u, ldu, tau, work, lwork, info)

!         Copy R to PH (used as workspace)
          Call zlacpy('Upper', n, n, a, lda, ph, ldph)

!         Set the strictly lower triangular part of R to zero
          Call nagf_blas_zmload('Lower', n-1, n-1, zero, zero, ph(2,1), ldph)

!         Bidiagonalize R
          Call zgebrd(n, n, ph, ldph, d, e, tauq, taup, work, lwork, info)

!         Update Q, storing the result in U
          Call zunmbr('Q', 'Right', 'No transpose', m, n, n, ph, ldph, tauq, &
            u, ldu, work, lwork, info)

!         Print bidiagonal form and matrix Q

          Write (nout, *)
          Write (nout, *) 'Example 1: bidiagonal matrix B'
          Write (nout, *) 'Diagonal'
          Write (nout, 100) d(1:n)
          Write (nout, *) 'Superdiagonal'
          Write (nout, 100) e(1:n-1)
          Write (nout, *)
          Flush (nout)

!         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', ' ', m, n, &
            u, ldu, 'Bracketed', 'F7.4', 'Example 1: matrix Q', 'Integer', &
            rlabs, 'Integer', clabs, 80, 0, ifail)

        Else

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

!         Copy A to PH
          Call zlacpy('Upper', m, n, a, lda, ph, ldph)

!         Form Q explicitly, storing the result in PH
          Call zunglq(n, n, m, ph, ldph, tau, work, lwork, info)

!         Copy L to U (used as workspace)
          Call zlacpy('Lower', m, m, a, lda, u, ldu)

!         Set the strictly upper triangular part of L to zero
          Call nagf_blas_zmload('Upper', m-1, m-1, zero, zero, u(1,2), ldu)

!         Bidiagonalize L
          Call zgebrd(m, m, u, ldu, d, e, tauq, taup, work, lwork, info)

!         Update P**H, storing the result in PH
          Call zunmbr('P', 'Left', 'Conjugate transpose', m, n, m, u, ldu, &
            taup, ph, ldph, work, lwork, info)

!         Print bidiagonal form and matrix P**H

          Write (nout, *)
          Write (nout, *) 'Example 2: bidiagonal matrix B'
          Write (nout, *) 'Diagonal'
          Write (nout, 100) d(1:m)
          Write (nout, *) 'Superdiagonal'
          Write (nout, 100) e(1:m-1)
          Write (nout, *)
          Flush (nout)

          ifail = 0
          Call nagf_file_print_matrix_complex_gen_comp('General', ' ', m, n, &
            ph, ldph, 'Bracketed', 'F7.4', 'Example 2: matrix P**H', &
            'Integer', rlabs, 'Integer', clabs, 80, 0, ifail)

        End If
        Deallocate (a, ph, tau, taup, tauq, u, work, d, e)
      End Do

100   Format (3X, (8F8.4))
    End Program