計算ルーチン: 実行列ペアの一般化特異値分解の前処理として直交行列を計算する : (BLAS-3 を用いて)

DGGSVP3 Example Program Data

  4    3    2   :Values of M, N and P

  1.0  2.0  3.0
  3.0  2.0  1.0
  4.0  5.0  6.0
  7.0  8.0  8.0 :End of matrix A

 -2.0 -3.0  3.0
  4.0  6.0  5.0 :End of matrix B

 DGGSVP3 Example Program Results

 Number of infinite generalized singular values (k)
     1
 Number of finite generalized singular values (l)
     2
 Effective Numerical rank of (A; B) (k+l)
     3

 Finite generalized singular values
     1.3151E+00  8.0185E-02

 Orthogonal matrix U
              1           2           3           4
 1  -1.3484E-01  5.2524E-01 -2.0924E-01  8.1373E-01
 2   6.7420E-01 -5.2213E-01 -3.8886E-01  3.4874E-01
 3   2.6968E-01  5.2757E-01 -6.5782E-01 -4.6499E-01
 4   6.7420E-01  4.1615E-01  6.1014E-01  1.5127E-15

 Orthogonal matrix V
              1           2
 1   3.5539E-01 -9.3472E-01
 2   9.3472E-01  3.5539E-01

 Orthogonal matrix Q
              1           2           3
 1  -8.3205E-01 -9.4633E-02 -5.4657E-01
 2   5.5470E-01 -1.4195E-01 -8.1985E-01
 3   0.0000E+00 -9.8534E-01  1.7060E-01

 Nonsingular upper triangular matrix R
              1           2           3
 1  -2.0569E+00 -9.0121E+00 -9.3705E+00
 2              -1.0882E+01 -7.2688E+00
 3                          -6.0405E+00

 Number of cycles of the Kogbetliantz method
     2

    Program dggsvp3_example

!     DGGSVP3 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_comp
      Use lapack_interfaces, Only: dggsvp3, dlange, dtgsja
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: eps, tola, tolb
      Integer :: i, ifail, info, irank, j, k, l, lda, ldb, ldq, ldu, ldv, &
        lwork, m, n, ncycle, p
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), alpha(:), b(:, :), beta(:), &
        q(:, :), tau(:), u(:, :), v(:, :), work(:)
      Real (Kind=dp) :: wdum(1)
      Integer, Allocatable :: iwork(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: epsilon, max, nint, real
!     .. Executable Statements ..
      Write (nout, *) 'DGGSVP3 Example Program Results'
      Write (nout, *)
      Flush (nout)

!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, p
      lda = m
      ldb = p
      ldq = n
      ldu = m
      ldv = p
      Allocate (a(lda,n), alpha(n), b(ldb,n), beta(n), q(ldq,n), tau(n), &
        u(ldu,m), v(ldv,p), iwork(n))

!     Perform workspace query to get optimal size of work
      lwork = -1
      Call dggsvp3('U', 'V', 'Q', m, p, n, a, lda, b, ldb, tola, tolb, k, l, &
        u, ldu, v, ldv, q, ldq, iwork, tau, wdum, lwork, info)
      lwork = nint(wdum(1))
      Allocate (work(lwork))

!     Read the m by n matrix A and p by n matrix B from data file

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

!     Compute tola and tolb as
!         tola = max(m,n)*norm(A)*macheps
!         tolb = max(p,n)*norm(B)*macheps

      eps = epsilon(1.0E0_dp)
      tola = real(max(m,n), kind=dp)*dlange('One-norm', m, n, a, lda, work)* &
        eps
      tolb = real(max(p,n), kind=dp)*dlange('One-norm', p, n, b, ldb, work)* &
        eps

!     Compute the factorization of (A, B)
!         (A = U*S*(Q**T), B = V*T*(Q**T))

      Call dggsvp3('U', 'V', 'Q', m, p, n, a, lda, b, ldb, tola, tolb, k, l, &
        u, ldu, v, ldv, q, ldq, iwork, tau, work, lwork, info)

!     Given the factors above find the generalized SVD of (A, B)

      Call dtgsja('U', 'V', 'Q', m, p, n, k, l, a, lda, b, ldb, tola, tolb, &
        alpha, beta, u, ldu, v, ldv, q, ldq, work, ncycle, info)

!     Print solution

      irank = k + l
      Write (nout, *) 'Number of infinite generalized singular values (k)'
      Write (nout, 100) k
      Write (nout, *) 'Number of finite generalized singular values (l)'
      Write (nout, 100) l
      Write (nout, *) 'Effective Numerical rank of (A; B) (k+l)'
      Write (nout, 100) irank
      Write (nout, *)
      Write (nout, *) 'Finite generalized singular values'
      Write (nout, 110)(alpha(j)/beta(j), j=k+1, irank)

      Write (nout, *)
      Flush (nout)

      Call nagf_file_print_matrix_real_gen_comp('General', ' ', m, m, u, ldu, &
        '1P,E12.4', 'Orthogonal matrix U', 'Integer', rlabs, 'Integer', clabs, &
        80, 0, ifail)

      Write (nout, *)
      Flush (nout)

      Call nagf_file_print_matrix_real_gen_comp('General', ' ', p, p, v, ldv, &
        '1P,E12.4', 'Orthogonal matrix V', 'Integer', rlabs, 'Integer', clabs, &
        80, 0, ifail)

      Write (nout, *)
      Flush (nout)

      Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, q, ldq, &
        '1P,E12.4', 'Orthogonal matrix Q', 'Integer', rlabs, 'Integer', clabs, &
        80, 0, ifail)

      Write (nout, *)
      Flush (nout)

      Call nagf_file_print_matrix_real_gen_comp('Upper triangular', &
        'Non-unit', irank, irank, a(1,n-irank+1), lda, '1P,E12.4', &
        'Nonsingular upper triangular matrix R', 'Integer', rlabs, 'Integer', &
        clabs, 80, 0, ifail)

      Write (nout, *)
      Write (nout, *) 'Number of cycles of the Kogbetliantz method'
      Write (nout, 100) ncycle

100   Format (1X, I5)
110   Format (3X, 8(1P,E12.4))
    End Program