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.