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_zggsvd (f08vn)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_lapack_zggsvd (f08vn) computes the generalized singular value decomposition (GSVD) of an m by n complex matrix A and a p by n complex matrix B.

Syntax

[k, l, a, b, alpha, beta, u, v, q, iwork, info] = f08vn(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = nag_lapack_zggsvd(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
Note: the interface to this routine has changed since earlier releases of the toolbox:
At Mark 23: iwork was made an output parameter

Description

The generalized singular value decomposition is given by
UH A Q = D1 0 R ,   VH B Q = D2 0 R ,  
where U, V and Q are unitary matrices. Let k+l be 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 structured as follows:
if m-k-l0,
D1= klkI0l0Cm-k-l00()  
D2= kll0Sp-l00()  
0R = n-k-lklk0R11R12l00R22()  
where
C = diagαk+1,,αk+l ,  
S = diagβk+1,,βk+l ,  
and
C2 + S2 = I .  
R is stored as a submatrix of A with elements Rij stored as Ai,n-k-l+j on exit.
If m-k-l<0 ,
D1= km-kk+l-mkI00m-k0C0()  
D2= km-kk+l-mm-k0S0k+l-m00Ip-l000()  
0R = n-k-lkm-kk+l-mk0R11R12R13m-k00R22R23k+l-m000R33()  
where
C = diagαk+1,,αm ,  
S = diagβk+1,,βm ,  
and
C2 + S2 = I .  
R11 R12 R13 0 R22 R23  is stored as a submatrix of A with Rij stored as Ai,n-k-l+j, and R33  is stored as a submatrix of B with R33ij stored as Bm-k+i,n+m-k-l+j.
The function computes C, S, R and, optionally, the unitary 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×B-1:
A B-1 = U D1 D2-1 VH .  
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:
AH Ax=λ BH Bx .  
In some literature, the GSVD of A and B is presented in the form
UH A X = 0D1 ,   VH B X = 0D2 ,  
where U and V are orthogonal and X is nonsingular, and 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 R-1 .  

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', the unitary matrix U is computed.
If jobu='N', U is not computed.
Constraint: jobu='U' or 'N'.
2:     jobv – string (length ≥ 1)
If jobv='V', the unitary matrix V is computed.
If jobv='N', V is not computed.
Constraint: jobv='V' or 'N'.
3:     jobq – string (length ≥ 1)
If jobq='Q', the unitary matrix Q is computed.
If jobq='N', Q is not computed.
Constraint: jobq='Q' or 'N'.
4:     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.
5:     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.

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:     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.
3:     p int64int32nag_int scalar
Default: the first dimension of the array b.
p, the number of rows of the matrix B.
Constraint: p0.

Output Parameters

1:     k int64int32nag_int scalar
2:     l int64int32nag_int scalar
k and l specify the dimension of the subblocks k and l as described in Description; k+l is the effective numerical rank of AB.
3:     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.
Contains the triangular matrix R, or part of R. See Description for details.
4:     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.
Contains the triangular matrix R if m-k-l<0. See Description for details.
5:     alphan – double array
See the description of beta.
6:     betan – double array
alpha and beta contain the generalized singular value pairs of A and B, αi and βi ;
  • alpha1:k = 1 ,
  • beta1:k = 0 ,
and if m-k-l0 ,
  • alphak+1:k+l = C ,
  • betak+1:k+l = S ,
or if m-k-l<0 ,
  • alphak+1:m = C ,
  • alpham+1:k+l = 0 ,
  • betak+1:m = S ,
  • betam+1:k+l = 1 , and
  • alphak+l+1:n = 0 ,
  • betak+l+1:n = 0 .
The notation alphak:n above refers to consecutive elements alphai, for i=k,,n.
7:     uldu: – complex array
The first dimension, ldu, of the array u will be
  • if jobu='U', ldu= max1,m ;
  • otherwise ldu=1.
The second dimension of the array u will be max1,m if jobu='U' and 1 otherwise.
If jobu='U', u contains the m by m unitary matrix U.
If jobu='N', u is not referenced.
8:     vldv: – complex array
The first dimension, ldv, of the array v will be
  • if jobv='V', ldv= max1,p ;
  • otherwise ldv=1.
The second dimension of the array v will be max1,p if jobv='V' and 1 otherwise.
If jobv='V', v contains the p by p unitary matrix V.
If jobv='N', v is not referenced.
9:     qldq: – complex array
The first dimension, ldq, of the array q will be
  • if jobq='Q', ldq= max1,n ;
  • otherwise ldq=1.
The second dimension of the array q will be max1,n if jobq='Q' and 1 otherwise.
If jobq='Q', q contains the n by n unitary matrix Q.
If jobq='N', q is not referenced.
10:   iworkn int64int32nag_int array
Stores the sorting information. More precisely, the following loop will sort alpha
 for i=k+1:min(m,k+l)
 % swap alpha(i) and alpha(iwork(i))
 end 
such that alpha1alpha2alphan.
11:   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: n, 6: p, 7: k, 8: l, 9: a, 10: lda, 11: b, 12: ldb, 13: alpha, 14: beta, 15: u, 16: ldu, 17: v, 18: ldv, 19: q, 20: ldq, 21: work, 22: rwork, 23: iwork, 24: 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
If info=1, the Jacobi-type procedure failed to converge.

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 diagonal elements of the matrix R are real.
The real analogue of this function is nag_lapack_dggsvd (f08va).

Example

This example finds the generalized singular value decomposition
A = U Σ1 0R QH ,   B = V Σ2 0R QH ,  
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 ,  
together with estimates for the condition number of R and the error bound for the computed generalized singular values.
The example program assumes that mn, and would need slight modification if this is not the case.
function f08vn_example


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

% Generalized SVD of (A, B), where
m = 6;
n = 4;
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]);

% Compute the GSVD of (a, b)
  % (a,b) = (U*D1,v*D2)*(0 R)*(Q^H), m>=n
jobu = 'U';
jobv = 'V';
jobq = 'Q';
[k, l, DR, b, alpha, beta, U, V, Q, iwork, info] = ...
  f08vn( ...
	 jobu, jobv, jobq, a, b);

irank = k + l;
fprintf('\nInfinite generalized singular values = %5d\n', k);
fprintf('  Finite generalized singular values = %5d\n', l);
fprintf('  Numerical rank of (a'' b'')''         = %5d\n', irank);

fprintf('\nFinite generalized singular values\n');
w = alpha(k+1:irank)./beta(k+1:irank);
disp(transpose(w));

fprintf('Orthogonal matrix u\n');
disp(U);

fprintf('Orthogonal matrix v\n');
disp(V);

fprintf('Orthogonal matrix q\n');
disp(Q);

fprintf('Non singular upper triangular matrix R\n');
R = DR(1:n,n-(irank-1):n);
disp(R);

% Estimate the reciprocal condition number of R
[rcond, info] = f07tu( ...
		       'Infinity-norm','Upper','Non-unit', R);

fprintf('Estimate of reciprocal condition number for R\n%11.1e\n', rcond);

% If irank = n, compute error bound for generalized singular values
if (irank==n)
  serrbd = x02aj/rcond;
  fprintf('\nError bound for the generalized singular values\n%11.1e\n', ...
          serrbd);
else
  fprintf('\n(A^H B^H)^H is not of full rank\n');
end


f08vn example results


Infinite generalized singular values =     2
  Finite generalized singular values =     2
  Numerical rank of (a' b')'         =     4

Finite generalized singular values
    2.0720    1.1058

Orthogonal matrix u
  -0.0130 - 0.3259i  -0.1404 - 0.2617i   0.2518 - 0.7979i  -0.0510 - 0.2175i  -0.0459 + 0.0001i  -0.0528 - 0.2249i
   0.4276 - 0.6258i   0.0863 - 0.0382i  -0.3219 + 0.1611i   0.1198 + 0.1632i  -0.0803 - 0.4360i  -0.0381 - 0.2191i
  -0.3259 + 0.1643i   0.3816 - 0.1822i   0.1323 - 0.0146i  -0.5067 + 0.1862i   0.0597 - 0.5897i  -0.1385 - 0.0909i
   0.1591 - 0.0052i  -0.2821 + 0.1973i   0.2160 + 0.1881i  -0.4016 + 0.2679i  -0.0464 + 0.3086i  -0.3735 - 0.5515i
  -0.1721 - 0.0130i  -0.5094 - 0.5032i   0.0365 + 0.2032i   0.1927 + 0.1557i   0.5784 - 0.1244i  -0.0188 - 0.0557i
  -0.2634 - 0.2477i  -0.1086 + 0.2847i   0.1091 - 0.1271i  -0.0882 + 0.5617i   0.0158 + 0.0471i   0.6501 + 0.0049i

Orthogonal matrix v
   0.9893 - 0.0000i  -0.1146 + 0.0903i
  -0.1146 - 0.0903i  -0.9893 - 0.0000i

Orthogonal matrix q
   0.7071 + 0.0000i   0.0000 + 0.0000i   0.6995 - 0.0000i   0.0810 - 0.0638i
   0.0000 + 0.0000i   0.7071 + 0.0000i  -0.0810 - 0.0638i   0.6995 + 0.0000i
   0.7071 + 0.0000i   0.0000 + 0.0000i  -0.6995 + 0.0000i  -0.0810 + 0.0638i
   0.0000 + 0.0000i   0.7071 + 0.0000i   0.0810 + 0.0638i  -0.6995 - 0.0000i

Non singular upper triangular matrix R
  -2.7118 + 0.0000i  -1.4390 - 1.0315i  -0.0769 + 1.3613i  -0.2814 - 0.0324i
   0.0000 + 0.0000i  -1.8583 + 0.0000i  -1.0760 + 0.0310i   1.3292 + 0.3677i
   0.0000 + 0.0000i   0.0000 + 0.0000i   3.2537 + 0.0000i  -0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -2.1084 + 0.0000i

Estimate of reciprocal condition number for R
    1.3e-01

Error bound for the generalized singular values
    8.3e-16

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