hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_lapack_ztgsja (f08ys)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_lapack_ztgsja (f08ys) computes the generalized singular value decomposition (GSVD) of two complex 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 nag_lapack_zggsvp (f08vs).

Syntax

[a, b, alpha, beta, u, v, q, ncycle, info] = f08ys(jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, 'm', m, 'p', p, 'n', n)
[a, b, alpha, beta, u, v, q, ncycle, info] = nag_lapack_ztgsja(jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, 'm', m, 'p', p, 'n', n)

Description

nag_lapack_ztgsja (f08ys) computes the GSVD of the matrices A and B which are assumed to have the form as returned by nag_lapack_zggsvp (f08vs)
A= n-k-lklk0A12A13l00A23m-k-l000() ,   if ​ m-k-l 0; n-k-lklk0A12A13m-k00A23() ,   if ​ m-k-l < 0 ; B= n-k-lkll00B13p-l000() ,  
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.
nag_lapack_ztgsja (f08ys) computes unitary matrices Q, U and V, diagonal matrices D1 and D2, and an upper triangular matrix R such that
UHAQ = D1 0 R ,   VHBQ = 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= klkI0l0Cm-k-l00() ,  
D2= kll0Sp-l00() ,  
R = klkR11R12l0R22() ,  
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-mkI00m-k0C0() ,  
D2= km-kk+l-mm-k0S0k+l-m00Ip-l000() ,  
R = km-kk+l-mkR11R12R13m-k0R22R23k+l-m00R33() ,  
where C=diagαk+1,,,,,,αm,  S=diagβk+1,,,,,,βm.
In both cases the diagonal matrix C has real non-negative diagonal elements, the diagonal matrix S has real 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.

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

Parameters

Compulsory Input Parameters

1:     jobu – string (length ≥ 1)
If jobu='U', u must contain a unitary matrix U1 on entry, and the product U1U is returned.
If jobu='I', u is initialized to the unit matrix, and the unitary matrix U is returned.
If jobu='N', U is not computed.
Constraint: jobu='U', 'I' or 'N'.
2:     jobv – string (length ≥ 1)
If jobv='V', v must contain a unitary matrix V1 on entry, and the product V1V is returned.
If jobv='I', v is initialized to the unit matrix, and the unitary matrix V is returned.
If jobv='N', V is not computed.
Constraint: jobv='V', 'I' or 'N'.
3:     jobq – string (length ≥ 1)
If jobq='Q', q must contain a unitary matrix Q1 on entry, and the product Q1Q is returned.
If jobq='I', q is initialized to the unit matrix, and the unitary matrix Q is returned.
If jobq='N', Q is not computed.
Constraint: jobq='Q', 'I' or 'N'.
4:     k int64int32nag_int scalar
5:     l int64int32nag_int scalar
k and l specify the sizes, k and l, of the subblocks of A and B, whose GSVD is to be computed by nag_lapack_ztgsja (f08ys).
6:     alda: – complex array
The first dimension of the array a must be at least max1,m.
The second dimension of the array a must be at least max1,n.
The m by n matrix A.
7:     bldb: – complex array
The first dimension of the array b must be at least max1,p.
The second dimension of the array b must be at least max1,n.
The p by n matrix B.
8:     tola – double scalar
9:     tolb – double scalar
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 nag_lapack_zggsvp (f08vs), say
tola=maxm,nAε, tolb=maxp,nBε,  
where ε  is the machine precision.
10:   uldu: – complex array
The first dimension, ldu, of the array u must satisfy
  • if jobu='U' or 'I', ldu max1,m ;
  • otherwise ldu1.
The second dimension of the array u must be at least max1,m if jobu='U' or 'I', and at least 1 otherwise.
If jobu='U', u must contain an m by m matrix U1 (usually the unitary matrix returned by nag_lapack_zggsvp (f08vs)).
11:   vldv: – complex array
The first dimension, ldv, of the array v must satisfy
  • if jobv='V' or 'I', ldv max1,p ;
  • otherwise ldv1.
The second dimension of the array v must be at least max1,p if jobv='V' or 'I', and at least 1 otherwise.
If jobv='V', v must contain an p by p matrix V1 (usually the unitary matrix returned by nag_lapack_zggsvp (f08vs)).
12:   qldq: – complex array
The first dimension, ldq, of the array q must satisfy
  • if jobq='Q' or 'I', ldq max1,n ;
  • otherwise ldq1.
The second dimension of the array q must be at least max1,n if jobq='Q' or 'I', and at least 1 otherwise.
If jobq='Q', q must contain an n by n matrix Q1 (usually the unitary matrix returned by nag_lapack_zggsvp (f08vs)).

Optional Input Parameters

1:     m int64int32nag_int scalar
Default: the first dimension of the array a.
m, the number of rows of the matrix A.
Constraint: m0.
2:     p int64int32nag_int scalar
Default: the first dimension of the array b.
p, the number of rows of the matrix B.
Constraint: p0.
3:     n int64int32nag_int scalar
Default: the second dimension of the array a.
n, the number of columns of the matrices A and B.
Constraint: n0.

Output Parameters

1:     alda: – complex array
The first dimension of the array a will be max1,m.
The second dimension of the array a will be max1,n.
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 .
2:     bldb: – complex array
The first dimension of the array b will be max1,p.
The second dimension of the array b will be max1,n.
If m-k-l<0 , bm-k+1:ln+m-k-l+1:n contains the submatrix R33 of R.
3:     alphan – double array
See the description of beta.
4:     betan – double array
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.
5:     uldu: – complex array
The first dimension, ldu, of the array u will be
  • if jobu='U' or 'I', ldu= max1,m ;
  • otherwise ldu=1.
The second dimension of the array u will be max1,m if jobu='U' or 'I' and 1 otherwise.
If jobu='U', u contains the product U1U.
If jobu='I', u contains the unitary matrix U.
If jobu='N', u is not referenced.
6:     vldv: – complex array
The first dimension, ldv, of the array v will be
  • if jobv='V' or 'I', ldv= max1,p ;
  • otherwise ldv=1.
The second dimension of the array v will be max1,p if jobv='V' or 'I' and 1 otherwise.
If jobv='I', v contains the unitary matrix V.
If jobv='V', v contains the product V1V.
If jobv='N', v is not referenced.
7:     qldq: – complex array
The first dimension, ldq, of the array q will be
  • if jobq='Q' or 'I', ldq= max1,n ;
  • otherwise ldq=1.
The second dimension of the array q will be max1,n if jobq='Q' or 'I' and 1 otherwise.
If jobq='I', q contains the unitary matrix Q.
If jobq='Q', q contains the product Q1Q.
If jobq='N', q is not referenced.
8:     ncycle int64int32nag_int scalar
The number of cycles required for convergence.
9:     info int64int32nag_int scalar
info=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   info=-i
If info=-i, parameter i had an illegal value on entry. The parameters are numbered as follows:
1: jobu, 2: jobv, 3: jobq, 4: m, 5: p, 6: n, 7: k, 8: l, 9: a, 10: lda, 11: b, 12: ldb, 13: tola, 14: tolb, 15: alpha, 16: beta, 17: u, 18: ldu, 19: v, 20: ldv, 21: q, 22: ldq, 23: work, 24: ncycle, 25: info.
It is possible that info refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.
   info=1
The procedure does not converge after 40 cycles.

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.

Further Comments

The real analogue of this function is nag_lapack_dtgsja (f08ye).

Example

This example finds the generalized singular value decomposition
A = UΣ1 0 R QH ,   B= VΣ2 0 R QH ,  
of the matrix pair A,B, where
A = 0.96-0.81i -0.03+0.96i -0.91+2.06i -0.05+0.41i -0.98+1.98i -1.20+0.19i -0.66+0.42i -0.81+0.56i 0.62-0.46i 1.01+0.02i 0.63-0.17i -1.11+0.60i 0.37+0.38i 0.19-0.54i -0.98-0.36i 0.22-0.20i 0.83+0.51i 0.20+0.01i -0.17-0.46i 1.47+1.59i 1.08-0.28i 0.20-0.12i -0.07+1.23i 0.26+0.26i  
and
B = 1 0 -1 0 0 1 0 -1 .  
function f08ys_example


fprintf('f08ys example results\n\n');

% Compute generalized SVD of (A,B)
A = [ 0.96-0.81i -0.03+0.96i -0.91+2.06i -0.05+0.41i;
     -0.98+1.98i -1.20+0.19i -0.66+0.42i -0.81+0.56i;
      0.62-0.46i  1.01+0.02i  0.63-0.17i -1.11+0.60i;
      0.37+0.38i  0.19-0.54i -0.98-0.36i  0.22-0.20i;
      0.83+0.51i  0.20+0.01i -0.17-0.46i  1.47+1.59i;
      1.08-0.28i  0.20-0.12i -0.07+1.23i  0.26+0.26i];
B = complex([ 1 0 -1  0;
              0 1  0 -1]); 
[m, n] = size(A);
p      = size(B,1);

% Reduce A and B to upper triangular form S = U^H A Q, T = V^H B Q
tola = max(m,n)*norm(A,1)*x02aj;
tolb = max(p,n)*norm(B,1)*x02aj;
[S, T, k, l, U, V, Q, info] = ...
  f08vs(...
        'U', 'V', 'Q', A, B, tola, tolb);

% Compute singular values
[S, T, alpha, beta, U, V, Q, ncycle, info] = ...
  f08ys('U', 'V', 'Q', k, l, S, T, tola, tolb, U, V, Q);

fprintf('Number of infinite generalized singular values = %3d\n',k);
fprintf('Number of   finite generalized singular values = %3d\n',l);
fprintf('Effective rank of the matrix pair  (A^T B^T)^T = %3d\n',k+l);
fprintf('Number of cycles of the Kogbetliantz method    = %3d\n\n',ncycle);
disp('Finite generalized singular values');
disp(alpha(k+1:k+l)./beta(k+1:k+l));


f08ys example results

Number of infinite generalized singular values =   2
Number of   finite generalized singular values =   2
Effective rank of the matrix pair  (A^T B^T)^T =   4
Number of cycles of the Kogbetliantz method    =   2

Finite generalized singular values
    2.0720
    1.1058


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015