計算ルーチン: DGEBRD により決まる準対角形への縮約からの直交変換の適用

DORMBR Example Program Data
  6  4                                      :Values of M and N, Example 1
 -0.57  -1.28  -0.39   0.25
 -1.93   1.08  -0.31  -2.14
  2.30   0.24   0.40  -0.35
 -1.93   0.64  -0.66   0.08
  0.15   0.30   0.15  -2.13
 -0.02   1.03  -1.43   0.50                 :End of matrix A
  4  6                                      :Values of M and N, Example 2
 -5.42   3.28  -3.68   0.27   2.06   0.46
 -1.65  -3.40  -3.20  -1.03  -4.06  -0.01
 -0.37   2.35   1.90   4.31  -1.76   1.13
 -3.15  -0.11   1.99  -2.70   0.26   4.50   :End of matrix A

 DORMBR Example Program Results

 Example 1: bidiagonal matrix B
 Diagonal
     3.6177 -2.4161  1.9213 -1.4265
 Superdiagonal
     1.2587 -1.5262  1.1895

 Example 1: matrix Q
          1       2       3       4
 1  -0.1576 -0.2690  0.2612  0.8513
 2  -0.5335  0.5311 -0.2922  0.0184
 3   0.6358  0.3495 -0.0250 -0.0210
 4  -0.5335  0.0035  0.1537 -0.2592
 5   0.0415  0.5572 -0.2917  0.4523
 6  -0.0055  0.4614  0.8585 -0.0532

 Example 2: bidiagonal matrix B
 Diagonal
    -7.7724  6.1573 -6.0576  5.7933
 Superdiagonal
     1.1926  0.5734 -1.9143

 Example 2: matrix P**T
          1       2       3       4       5       6
 1  -0.7104  0.4299 -0.4824  0.0354  0.2700  0.0603
 2   0.3583  0.1382 -0.4110  0.4044  0.0951 -0.7148
 3  -0.0507  0.4244  0.3795  0.7402 -0.2773  0.2203
 4   0.2442  0.4016  0.4158 -0.1354  0.7666 -0.0137

    Program dormbr_example

!     DORMBR Example Program Text

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

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen
      Use lapack_interfaces, Only: dgebrd, dgelqf, dgeqrf, dlacpy, dlaset, &
        dorglq, dorgqr, dormbr
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=dp), Parameter :: zero = 0.0E0_dp
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, ic, ifail, info, lda, ldpt, ldu, lwork, m, n
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), d(:), e(:), pt(:, :), tau(:), &
        taup(:), tauq(:), u(:, :), work(:)
!     .. Executable Statements ..
      Write (nout, *) 'DORMBR Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Do ic = 1, 2
        Read (nin, *) m, n
        lda = m
        ldpt = n
        ldu = m
        lwork = 64*(m+n)
        Allocate (a(lda,n), d(n), e(n-1), pt(ldpt,n), tau(n), taup(n), &
          tauq(n), u(ldu,n), work(lwork))

!       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 dgeqrf(m, n, a, lda, tau, work, lwork, info)

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

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

!         Copy R to PT (used as workspace)
          Call dlacpy('Upper', n, n, a, lda, pt, ldpt)

!         Set the strictly lower triangular part of R to zero
          Call dlaset('Lower', n-1, n-1, zero, zero, pt(2,1), ldpt)

!         Bidiagonalize R
          Call dgebrd(n, n, pt, ldpt, d, e, tauq, taup, work, lwork, info)

!         Update Q, storing the result in U
          Call dormbr('Q', 'Right', 'No transpose', m, n, n, pt, ldpt, 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_real_gen('General', ' ', m, n, u, ldu, &
            'Example 1: matrix Q', ifail)

        Else

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

!         Copy A to PT
          Call dlacpy('Upper', m, n, a, lda, pt, ldpt)

!         Form Q explicitly, storing the result in PT
          Call dorglq(n, n, m, pt, ldpt, tau, work, lwork, info)

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

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

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

!         Update P**T, storing the result in PT
          Call dormbr('P', 'Left', 'Transpose', m, n, m, u, ldu, taup, pt, &
            ldpt, work, lwork, info)

!         Print bidiagonal form and matrix P**T

          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_real_gen('General', ' ', m, n, pt, ldpt, &
            'Example 2: matrix P**T', ifail)

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

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