!-----------------------------------------------------------------
!
! Purpose
! =======
!
!  LA_GGSVD computes the generalized singular value decomposition
!  (GSVD) of an M-by-N real matrix A and P-by-N real matrix B:
!
!      U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )
!
!  where U, V and Q are orthogonal matrices, and Z' is the transpose
!  of Z.  Let K+L = the effective numerical rank of the matrix (A',B')',
!  then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
!  D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
!  following structures, respectively:
!
!  If M-K-L >= 0,
!
!                      K  L
!         D1 =     K ( I  0 )
!                  L ( 0  C )
!              M-K-L ( 0  0 )
!
!                    K  L
!         D2 =   L ( 0  S )
!              P-L ( 0  0 )
!
!                  N-K-L  K    L
!    ( 0 R ) = K (  0   R11  R12 )
!              L (  0    0   R22 )
!
!  where
!
!    C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
!    S = diag( BETA(K+1),  ... , BETA(K+L) ),
!    C**2 + S**2 = I.
!
!    R is stored in A(1:K+L,N-K-L+1:N) on exit.
!
!  If M-K-L < 0,
!
!                    K M-K K+L-M
!         D1 =   K ( I  0    0   )
!              M-K ( 0  C    0   )
!
!                      K M-K K+L-M
!         D2 =   M-K ( 0  S    0  )
!              K+L-M ( 0  0    I  )
!                P-L ( 0  0    0  )
!
!                     N-K-L  K   M-K  K+L-M
!    ( 0 R ) =     K ( 0    R11  R12  R13  )
!                M-K ( 0     0   R22  R23  )
!              K+L-M ( 0     0    0   R33  )
!
!  where
!
!    C = diag( ALPHA(K+1), ... , ALPHA(M) ),
!    S = diag( BETA(K+1),  ... , BETA(M) ),
!    C**2 + S**2 = I.
!
!    (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 orthogonal
!  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 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 )
!  where U and V are orthogonal and X is nonsingular, D1 and D2 are
!  ``diagonal''.  The former GSVD form can be converted to the latter
!  form by taking the nonsingular matrix X as
!
!                       X = Q*( I   0    )
!                             ( 0 inv(R) ).
!
! =========
!
!   SUBROUTINE LA_GGSVD( A, B, ALPHA, BETA, K, L, U, V, Q, INFO )
!      <type>(<wp>), INTENT(INOUT) :: A(:,:), B(:,:)
!      REAL(<wp>), INTENT(OUT) :: ALPHA(:), BETA(:)
!      INTEGER, INTENT(OUT), OPTIONAL :: K, L
!      <type>(<wp>), INTENT(OUT), OPTIONAL :: U(:,:), V(:,:), Q(:,:)
!      INTEGER, INTENT(OUT), OPTIONAL :: INFO
!      <type>   ::= REAL | COMPLEX
!      <wp>     ::= KIND(1.0) | KIND(1.0D0)
!
! Arguments
! =========
!
! A    (input/output) REAL / COMPLEX array, shape (:,:),
!      SIZE(A,1) == m, SIZE(A,2) == 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.
!
! B    (input/output) REAL / COMPLEX array, shape (:,:),
!      SIZE(A,1) == p, SIZE(A,2) == n.
!      On entry, the p-by-n matrix B.
!      On exit, B contains the triangular matrix R if m-K-L < 0.
!      See Purpose for details.
!
! ALPHA(output) REAL array, shape (:), SIZE(ALPHA) == n.
! BETA (output) REAL array, SIZE(BETA) == n.
!      On exit, ALPHA and BETA contain the generalized singular
!      value pairs of A and B;
!        ALPHA(1:K) = 1,
!        BETA(1:K)  = 0,
!      and if m-K-L >= 0,
!        ALPHA(K+1:K+L) = C,
!        BETA(K+1:K+L)  = S,
!      or if m-K-L < 0,
!        ALPHA(K+1:m)=C, ALPHA(m+1:K+L)=0
!        BETA(K+1:m) =S, BETA(m+1:K+L) =1
!      and
!        ALPHA(K+L+1:n) = 0
!        BETA(K+L+1:n)  = 0
!
! K    Optional (output) INTEGER
! L    Optional (output) INTEGER
!      On exit, K and L specify the dimension of the subblocks
!      described in the Purpose section.
!      K + L = effective numerical rank of (A',B')'.
!
! U    Optional (output) REAL / COMPLEX array, shape (:,:),
!      SIZE(U,1) == SIZE(U,2) == m.
!      U contains the m-by-m orthogonal matrix U.
!
! V    Optional (output) REAL / COMPLEX array, shape (:,:),
!      SIZE(AV1) == SIZE(V,2) == p.
!      V contains the p-by-p orthogonal matrix V.
!
! Q    Optional (output) REAL / COMPLEX array, shape (:,:),
!      SIZE(Q,1) == SIZE(Q,2) == n.
!      Q contains the n-by-n orthogonal matrix Q.
!
! INFO Optional (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 STGSJA.
!      If INFO is not present and an error occurs, then the program is
!         terminated with an error message.
!
!  Internal Parameters
!  ===================
!
!  TOLA    REAL
!  TOLB    REAL
!          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)*MACHEPS,
!                   TOLB = MAX(P,N)*norm(B)*MACHEPS.
!          The size of TOLA and TOLB may affect the size of backward
!          errors of the decomposition.
!
!  Further Details
!  ===============
!
!  2-96 Based on modifications by
!     Ming Gu and Huan Ren, Computer Science Division, University of
!     California at Berkeley, USA
!-----------------------------------------------------------------

.