Previous: zggrqf Up: ../lapack-z.html Next: zggsvp


zggsvd


 NAME
      ZGGSVD - compute the generalized singular value decomposi-
      tion (GSVD) of the M-by-N complex matrix A and P-by-N com-
      plex matrix B

 SYNOPSIS
      SUBROUTINE ZGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA,
                         B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q,
                         LDQ, WORK, RWORK, IWORK, INFO )

          CHARACTER      JOBQ, JOBU, JOBV

          INTEGER        INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M,
                         N, P

          INTEGER        IWORK( * )

          DOUBLE         PRECISION ALPHA( * ), BETA( * ), RWORK( *
                         )

          COMPLEX*16     A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U(
                         LDU, * ), V( LDV, * ), WORK( * )

 PURPOSE
      ZGGSVD computes the generalized singular value decomposition
      (GSVD) of the M-by-N complex matrix A and P-by-N complex
      matrix B:

            U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )
      (1)

      where U, V and Q are unitary matrices, R is an upper tri-
      angular matrix, and Z' means the conjugate transpose of Z.
      Let K+L = the numerical effective rank of the matrix
      (A',B')', then D1 and D2 are M-by-(K+L) and P-by-(K+L)
      "diagonal" matrices and of the following structures, respec-
      tively:

      If M-K-L >= 0,

         U'*A*Q = D1*( 0 R )

                = K     ( I  0 ) * (  0   R11  R12 ) K
                  L     ( 0  C )   (  0    0   R22 ) L
                  M-K-L ( 0  0 )    N-K-L  K    L
                          K  L

         V'*B*Q = D2*( 0 R )

                = L     ( 0  S ) * (  0   R11  R12 ) K
                  P-L   ( 0  0 )   (  0    0   R22 ) L
                          K  L      N-K-L  K    L

      where

        C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
        S = diag( BETA(K+1),  ... , BETA(K+L) ), C**2 + S**2 = I.
        The nonsingular triangular matrix R = ( R11 R12 ) is
      stored
                                              (  0  R22 )
        in A(1:K+L,N-K-L+1:N) on exit.

      If M-K-L < 0,

         U'*A*Q = D1*( 0 R )

                = K   ( I  0    0   ) * ( 0    R11  R12  R13  ) K
                  M-K ( 0  C    0   )   ( 0     0   R22  R23  )
      M-K
                        K M-K K+L-M     ( 0     0    0   R33  )
      K+L-M
                                         N-K-L  K   M-K  K+L-M

         V'*B*Q = D2*( 0 R )

                = M-K   ( 0  S    0   ) * ( 0    R11  R12  R13  )
      K
                  K+L-M ( 0  0    I   )   ( 0     0   R22  R23  )
      M-K
                  P-L   ( 0  0    0   )   ( 0     0    0   R33  )
      K+L-M
                          K M-K K+L-M      N-K-L  K   M-K  K+L-M
      where

        C = diag( ALPHA(K+1), ... , ALPHA(M) ),
        S = diag( BETA(K+1),  ... , BETA(M) ), C**2 + S**2 = I.
        R = ( R11 R12 R13 ) is a nonsingular upper triangular
      matrix,
            (  0  R22 R23 )
            (  0   0  R33 )
        (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is
      stored
        ( 0  R22 R23 )
        in B(M-K+1:L,N+M-K-L+1:N) on exit.

      The routine computes C, S, R, and optionally the unitary
      transformation matrices U, V and Q.

      In particular, if B is an N-by-N nonsingular matrix, then
      the GSVD of A and B implicitly gives the SVD of the matrix
      A*inv(B):
                           A*inv(B) = U*(D1*inv(D2))*V'.
      If ( A',B')' has orthonormal columns, then the GSVD of A and
      B is also equal to the CS decomposition of A and B.

      Furthermore, the GSVD can be used to derive the solution of
      the eigenvalue problem:
                           A'*A x = lambda* B'*B x.
      In some literature, the GSVD of A and B is presented in the
      form
                       U'*A*X = ( 0 D1 ),   V'*B*X = ( 0 D2 )
      (2) where U and V are orthogonal and X is nonsingular, and
      D1 and D2 are ``diagonal''.  It is easy to see that the GSVD
      form (1) can be converted to the form (2) by taking the non-
      singular matrix X as

                            X = Q*(  I   0    )
                                  (  0 inv(R) )

 ARGUMENTS
      JOBU    (input) CHARACTER*1
              = 'U':  Unitary matrix U is computed;
              = 'N':  U is not computed.

      JOBV    (input) CHARACTER*1
              = 'V':  Unitary matrix V is computed;
              = 'N':  V is not computed.

      JOBQ    (input) CHARACTER*1
              = 'Q':  Unitary matrix Q is computed;
              = 'N':  Q is not computed.

      M       (input) INTEGER
              The number of rows of the matrix A.  M >= 0.

      N       (input) INTEGER
              The number of columns of the matrices A and B.  N >=
              0.

      P       (input) INTEGER
              The number of rows of the matrix B.  P >= 0.

      K       (output) INTEGER
              L       (output) INTEGER On exit, K and L specify
              the dimension of the subblocks described in Purpose.
              K + L = effective numerical rank of (A',B')'.

      A       (input/output) COMPLEX*16 array, dimension (LDA,N)
              On entry, the M-by-N matrix A.  On exit, A contains
              the triangular matrix R, or part of R.  See Purpose
              for details.

      LDA     (input) INTEGER
              The leading dimension of the array A. LDA >=
              max(1,M).

      B       (input/output) COMPLEX*16 array, dimension (LDB,N)
              On entry, the P-by-N matrix B.  On exit, B contains
              part of the triangular matrix R if M-K-L < 0.  See
              Purpose for details.

      LDB     (input) INTEGER
              The leading dimension of the array B. LDB >=
              max(1,P).

      ALPHA   (output) DOUBLE PRECISION array, dimension (N)
              BETA    (output) DOUBLE PRECISION array, dimension
              (N) On exit, ALPHA and BETA contain the generalized
              singular value pairs of A and B; if M-K-L >= 0,
              ALPHA(1:K) = ONE,  ALPHA(K+1:K+L) = C,
              BETA(1:K)  = ZERO, BETA(K+1:K+L)  = S; or if M-K-L <
              0, ALPHA(1:K)= ONE,  ALPHA(K+1:M)= C,
              ALPHA(M+1:K+L)= ZERO,
              BETA(1:K) = ZERO, BETA(K+1:M) = S, BETA(M+1:K+L) =
              ONE.  and ALPHA(K+L+1:N) = ZERO
              BETA(K+L+1:N)  = ZERO

      U       (output) COMPLEX*16 array, dimension (LDU,M)
              If JOBU = 'U', U contains the M-by-M unitary matrix
              U.  If JOBU = 'N', U is not referenced.

      LDU     (input) INTEGER
              The leading dimension of the array U. LDU >=
              max(1,M).

      V       (output) COMPLEX*16 array, dimension (LDV,P)
              If JOBU = 'V', V contains the P-by-P unitary matrix
              V.  If JOBV = 'N', V is not referenced.

      LDV     (input) INTEGER
              The leading dimension of the array V. LDV >=
              max(1,P).

      Q       (output) COMPLEX*16 array, dimension (LDQ,N)
              If JOBU = 'Q', Q contains the N-by-N unitary matrix
              Q.  If JOBQ = 'N', Q is not referenced.

      LDQ     (input) INTEGER
              The leading dimension of the array Q. LDQ >=
              max(1,N).

      WORK    (workspace) COMPLEX*16 array, dimension (MAX(3*N,M,P)+N)

      RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)

      IWORK   (workspace) INTEGER array, dimension (N)

      INFO    (output)INTEGER

              = 0:  successful exit.
              < 0:  if INFO = -i, the i-th argument had an illegal
              value.
              > 0:  if INFO = 1, the Jacobi-type procedure failed
              to converge.  For further details, see subroutine
              ZTGSJA.

 PARAMETERS
      TOLA    DOUBLE PRECISION
              TOLB    DOUBLE PRECISION TOLA and TOLB are the
              thresholds to determine the effective rank of
              (A',B')'. Generally, they are set to TOLA =
              MAX(M,N)*norm(A)*MAZHEPS, TOLB =
              MAX(P,N)*norm(B)*MAZHEPS.  The size of TOLA and TOLB
              may affect the size of backward errors of the decom-
              position.