F08YEF (DTGSJA) (PDF version)
F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

F08YEF (DTGSJA)

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

F08YEF (DTGSJA) computes the generalized singular value decomposition (GSVD) of two real upper trapezoidal matrices A and B, where A is an m by n matrix and B is a p by n matrix.
A and B are assumed to be in the form returned by F08VEF (DGGSVP).

2  Specification

SUBROUTINE F08YEF ( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, NCYCLE, INFO)
INTEGER  M, P, N, K, L, LDA, LDB, LDU, LDV, LDQ, NCYCLE, INFO
REAL (KIND=nag_wp)  A(LDA,*), B(LDB,*), TOLA, TOLB, ALPHA(N), BETA(N), U(LDU,*), V(LDV,*), Q(LDQ,*), WORK(2*N)
CHARACTER(1)  JOBU, JOBV, JOBQ
The routine may be called by its LAPACK name dtgsja.

3  Description

F08YEF (DTGSJA) computes the GSVD of the matrices A and B which are assumed to have the form as returned by F08VEF (DGGSVP)
A= n-k-lklk(0A12A13) l 0 0 A23 m-k-l 0 0 0 ,   if ​ m-k-l 0; n-k-lklk(0A12A13) m-k 0 0 A23 ,   if ​ m-k-l < 0 ; B= n-k-lkll(00B13) p-l 0 0 0 ,
where the k by k matrix A12 and the l by l matrix B13 are nonsingular upper triangular, A23 is l by l upper triangular if m-k-l0 and is m-k by l upper trapezoidal otherwise.
F08YEF (DTGSJA) computes orthogonal matrices Q, U and V, diagonal matrices D1 and D2, and an upper triangular matrix R such that
UTAQ = D1 0 R ,   VTBQ = D2 0 R .
Optionally Q, U and V may or may not be computed, or they may be premultiplied by matrices Q1, U1 and V1 respectively.
If m-k-l0 then D1, D2 and R have the form
D1= klk(I0) l 0 C m-k-l 0 0 ,
D2= kll(0S) p-l 0 0 ,
R = klk(R11R12) l 0 R22 ,
where C=diagαk+1,,,,,,αk+l,  S=diagβk+1,,,,,,βk+l.
If m-k-l<0 then D1, D2 and R have the form
D1= km-kk+l-mk(I00) m-k 0 C 0 ,
D2= km-kk+l-mm-k(0S0) k+l-m 0 0 I p-l 0 0 0 ,
R = km-kk+l-mk(R11R12R13) m-k 0 R22 R23 k+l-m 0 0 R33 ,
where C=diagαk+1,,,,,,αm,  S=diagβk+1,,,,,,βm.
In both cases the diagonal matrix C has non-negative diagonal elements, the diagonal matrix S has positive diagonal elements, so that S is nonsingular, and C2+S2=1. See Section 2.3.5.3 of Anderson et al. (1999) for further information.

4  References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

5  Parameters

1:     JOBU – CHARACTER(1)Input
On entry: if JOBU='U', U must contain an orthogonal matrix U1 on entry, and the product U1U is returned.
If JOBU='I', U is initialized to the unit matrix, and the orthogonal matrix U is returned.
If JOBU='N', U is not computed.
Constraint: JOBU='I', 'U' or 'N'.
2:     JOBV – CHARACTER(1)Input
On entry: if JOBV='V', V must contain an orthogonal matrix V1 on entry, and the product V1V is returned.
If JOBV='I', V is initialized to the unit matrix, and the orthogonal matrix V is returned.
If JOBV='N', V is not computed.
Constraint: JOBV='V', 'I' or 'N'.
3:     JOBQ – CHARACTER(1)Input
On entry: if JOBQ='Q', Q must contain an orthogonal matrix Q1 on entry, and the product Q1Q is returned.
If JOBQ='I', Q is initialized to the unit matrix, and the orthogonal matrix Q is returned.
If JOBQ='N', Q is not computed.
Constraint: JOBQ='Q', 'I' or 'N'.
4:     M – INTEGERInput
On entry: m, the number of rows of the matrix A.
Constraint: M0.
5:     P – INTEGERInput
On entry: p, the number of rows of the matrix B.
Constraint: P0.
6:     N – INTEGERInput
On entry: n, the number of columns of the matrices A and B.
Constraint: N0.
7:     K – INTEGERInput
8:     L – INTEGERInput
On entry: K and L specify the sizes, k and l, of the subblocks of A and B, whose GSVD is to be computed by F08YEF (DTGSJA).
9:     A(LDA,*) – REAL (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array A must be at least max1,N.
On entry: the m by n matrix A.
On exit: if m-k-l0, A1:k+ln-k-l+1:n  contains the k+l by k+l upper triangular matrix R.
If m-k-l<0, A1:mn-k-l+1:n  contains the first m rows of the k+l by k+l upper triangular matrix R, and the submatrix R33 is returned in Bm-k+1:ln+m-k-l+1:n .
10:   LDA – INTEGERInput
On entry: the first dimension of the array A as declared in the (sub)program from which F08YEF (DTGSJA) is called.
Constraint: LDAmax1,M.
11:   B(LDB,*) – REAL (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array B must be at least max1,N.
On entry: the p by n matrix B.
On exit: if m-k-l<0 , Bm-k+1:ln+m-k-l+1:n  contains the submatrix R33 of R.
12:   LDB – INTEGERInput
On entry: the first dimension of the array B as declared in the (sub)program from which F08YEF (DTGSJA) is called.
Constraint: LDBmax1,P.
13:   TOLA – REAL (KIND=nag_wp)Input
14:   TOLB – REAL (KIND=nag_wp)Input
On entry: TOLA and TOLB are the convergence criteria for the Jacobi–Kogbetliantz iteration procedure. Generally, they should be the same as used in the preprocessing step performed by F08VSF (ZGGSVP), say
TOLA=maxM,NAε, TOLB=maxP,NBε,
where ε  is the machine precision.
15:   ALPHA(N) – REAL (KIND=nag_wp) arrayOutput
On exit: see the description of BETA.
16:   BETA(N) – REAL (KIND=nag_wp) arrayOutput
On exit: ALPHA and BETA contain the generalized singular value pairs of A and B;
  • ALPHAi=1 , BETAi=0 , for i=1,2,,k, and
  • if m-k-l0 , ALPHAi=αi , BETAi=βi , for i=k+1,,k+l, or
  • if m-k-l<0 , ALPHAi=αi , BETAi=βi , for i=k+1,,m and ALPHAi=0 , BETAi=1 , for i=m+1,,k+l.
Furthermore, if k+l<n, ALPHAi= BETAi=0 , for i=k+l+1,,n.
17:   U(LDU,*) – REAL (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array U must be at least max1,M if JOBU='U' or 'I', and at least 1 otherwise.
On entry: if JOBU='U', U must contain an m by m matrix U1 (usually the orthogonal matrix returned by F08VEF (DGGSVP)).
On exit: if JOBU='I', U contains the orthogonal matrix U.
If JOBU='U', U contains the product U1U.
If JOBU='N', U is not referenced.
18:   LDU – INTEGERInput
On entry: the first dimension of the array U as declared in the (sub)program from which F08YEF (DTGSJA) is called.
Constraints:
  • if JOBU'N', LDU max1,M ;
  • otherwise LDU1.
19:   V(LDV,*) – REAL (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array V must be at least max1,P if JOBV='V' or 'I', and at least 1 otherwise.
On entry: if JOBV='V', V must contain an p by p matrix V1 (usually the orthogonal matrix returned by F08VEF (DGGSVP)).
On exit: if JOBV='I', V contains the orthogonal matrix V.
If JOBV='V', V contains the product V1V.
If JOBV='N', V is not referenced.
20:   LDV – INTEGERInput
On entry: the first dimension of the array V as declared in the (sub)program from which F08YEF (DTGSJA) is called.
Constraints:
  • if JOBV'N', LDV max1,P ;
  • otherwise LDV1.
21:   Q(LDQ,*) – REAL (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array Q must be at least max1,N if JOBQ='Q' or 'I', and at least 1 otherwise.
On entry: if JOBQ='Q', Q must contain an n by n matrix Q1 (usually the orthogonal matrix returned by F08VEF (DGGSVP)).
On exit: if JOBQ='I', Q contains the orthogonal matrix Q.
If JOBQ='Q', Q contains the product Q1Q.
If JOBQ='N', Q is not referenced.
22:   LDQ – INTEGERInput
On entry: the first dimension of the array Q as declared in the (sub)program from which F08YEF (DTGSJA) is called.
Constraints:
  • if JOBQ'N', LDQ max1,N ;
  • otherwise LDQ1.
23:   WORK(2×N) – REAL (KIND=nag_wp) arrayWorkspace
24:   NCYCLE – INTEGEROutput
On exit: the number of cycles required for convergence.
25:   INFO – INTEGEROutput
On exit: INFO=0 unless the routine detects an error (see Section 6).

6  Error Indicators and Warnings

Errors or warnings detected by the routine:
INFO<0
If INFO=-i, argument i had an illegal value. An explanatory message is output, and execution of the program is terminated.
INFO=1
The procedure does not converge after 40 cycles.

7  Accuracy

The computed generalized singular value decomposition is nearly the exact generalized singular value decomposition for nearby matrices A+E and B+F, where
E2 = Oε A2   and   F2= Oε B2 ,
and ε is the machine precision. See Section 4.12 of Anderson et al. (1999) for further details.

8  Further Comments

The complex analogue of this routine is F08YSF (ZTGSJA).

9  Example

This example finds the generalized singular value decomposition
A = UΣ1 0 R QT ,   B= VΣ2 0 R QT ,
of the matrix pair A,B, where
A = 1 2 3 3 2 1 4 5 6 7 8 8   and   B= -2 -3 3 4 6 5 .

9.1  Program Text

Program Text (f08yefe.f90)

9.2  Program Data

Program Data (f08yefe.d)

9.3  Program Results

Program Results (f08yefe.r)


F08YEF (DTGSJA) (PDF version)
F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2012