Purpose ======= LA_GGSVD computes the generalized singular values and, optionally, the transformation matrices from the generalized singular value decomposition (GSVD) of a real or complex matrix pair (A,B), where A is m by n and B is p by n. The GSVD of (A,B) is written A = U * SIGMA1(0, R)*Q^H , B = V * SIGMA2(0, R)*Q^H where U , V and Q are orthogonal (unitary) matrices of dimensions m by m, p by p and n by n, respectively. Let l be the rank of B and r the rank of the (m + p) * n matrix ( A ) ( B ) , and let k = r-l. Then SIGMA1 and SIGMA2 are m*(k + l) and p * (k + l) "diagonal" matrices, respectively, and R is a (k + l) * (k + l) nonsingular triangular matrix. The detailed structure of SIGMA1 ,SIGMA2 and R depends on the sign of (m - k - l) as follows: The case m-k-l>=0: k l k ( I 0 ) SIGMA1 = l ( 0 C ) m-k-l ( 0 0 ) k l SIGMA2 = l ( 0 S ) p - l ( 0 S ) n-k-l k l (0, R) = k ( 0 R11 R12 ) l ( 0 0 R22 ) where C^2 + S^2 = I . We define alpha(1)=alpha(2)=...=alpha(k) = 1, alpha(k+i)=c(i i), i=1,2,...,l beta(1) = beta(2) = ... = beta(k) = 0, beta(k+i) = s(i i), i=1,2,...,l The case m-k-l < 0: k m-k k+l-m SIGMA1 = k ( I 0 0 ) m-k ( 0 C 0 ) k m-k k+l-m m-k ( 0 S 0 ) SIGMA2 = k+l-m ( 0 0 I ) p-l ( 0 0 0 ) n-k-l k m-k k+l-m k ( 0 R11 R12 R13 ) (0,R) = m-k ( 0 0 R22 R23 ) k+l-m ( 0 0 0 R33 ) where C^2 + S^2 = I . We define alpha(1)=alpha(2)=...=alpha(k) = 1, alpha(k+i)=c(i i), i =1,2,...,m-k, alpha(m+1) = alpha(m+2)=...= alpha(k+l) = 0 beta(1)=beta(2)= ... =beta(k)=0, beta(k+i)=s(i i), i=1,2,...,m-k, beta(m+1) = beta(m+2) = ... = beta(k+l) = 1 In both cases the generalized singular values of the pair (A,B) are the ratios sigma(i) = alpha(i)/beta(i), i = 1,2, ... ,k+l The first k singular values are infinite. The finite singular values are real and nonnegative. LA_GGSVD computes the real (nonnegative) scalars alpha(i), beta(i), i=1,2,..., k+l , the matrix R, and, optionally, the transformation matrices U , V and Q. ========= SUBROUTINE LA_GGSVD( A, B, ALPHA, BETA, K=k, L=l, & U=u, V=v, Q=q, IWORK=iwork, INFO=info ) (), INTENT(INOUT) :: A(:,:), B(:,:) REAL(), INTENT(OUT) :: ALPHA(:), BETA(:) INTEGER, INTENT(OUT), OPTIONAL :: K, L (), INTENT(OUT), OPTIONAL :: U(:,:), V(:,:), Q(:,:) INTEGER, INTENT(IN), OPTIONAL :: IWORK(:) INTEGER, INTENT(OUT), OPTIONAL :: INFO where ::= REAL | COMPLEX ::= KIND(1.0) | KIND(1.0D0) Arguments ========= A (input/output) REAL or COMPLEX array, shape (:,:) with size(A,1) = m and size(A,2) = n. On entry, the matrix A. On exit, A contains the triangular matrix R, or part of R, as follows: If m-k-l >= 0, then R is stored in A(1:k+l,n-k-l+1:n). If m-k-l < 0, then the matrix ( R11 R12 R13 ) ( 0 R22 R23 ) is stored in A(1:m,n-k-l+1:n). B (input/output) REAL or COMPLEX array, shape (:,:) with size(B,1) = p and size(B,2) = n. On entry, the matrix B. On exit, if m-k-l < 0, then R33 is stored in B(m-k+1:l,n+m-k-l+1:n). ALPHA (output) REAL array, shape (:) with size(ALPHA) = n The real scalars alpha(i) , i = 1, 2,..., k+l. BETA (output) REAL array, shape (:) with size(BETA) = n. The real scalars beta(i) , i = 1, 2, ..., k+l. Note: The generalized singular values of the pair (A,B) are sigma(i) = ALPHA(i)/BETA(i), i = 1, 2, ..., k+l. If k + l < n, then ALPHA(k+l+1:n) = BETA(k+l+1:n) = 0. K, L Optional (output) INTEGER. The dimension parameters k and l. U Optional (output) REAL or COMPLEX square array, shape (:,:) with size(U,1) = m. The matrix U . V Optional (output) REAL or COMPLEX square array, shape (:,:) with size(V,1) = p. The matrix V . Q Optional (output) REAL or COMPLEX square array, shape (:,:) with size(Q,1) = n. The matrix Q. IWORK Optional (output) INTEGER array, shape(:) with size(IWORK) = n. IWORK contains sorting information. More precisely, the loop for i = k + 1, min(m, k + l) swap ALPHA(i) and ALPHA(IWORK(i)) end will sort ALPHA so that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(n). INFO Optional (output) INTEGER. = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, the algorithm failed to converge. If INFO is not present and an error occurs, then the program is terminated with an error message.