F08RNF (ZUNCSD) (PDF version)
F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

F08RNF (ZUNCSD)

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

F08RNF (ZUNCSD) computes the CS decomposition of a complex m by m unitary matrix X, partitioned into a 2 by 2 array of submatrices.

2  Specification

SUBROUTINE F08RNF ( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
INTEGER  M, P, Q, LDX11, LDX12, LDX21, LDX22, LDU1, LDU2, LDV1T, LDV2T, LWORK, LRWORK, IWORK(M-min(P,M-P,Q,M-Q)), INFO
REAL (KIND=nag_wp)  THETA(min(P,M-P,Q,M-Q)), RWORK(max(1,LRWORK))
COMPLEX (KIND=nag_wp)  X11(LDX11,*), X12(LDX12,*), X21(LDX21,*), X22(LDX22,*), U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), V2T(LDV2T,*), WORK(max(1,LWORK))
CHARACTER(1)  JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS
The routine may be called by its LAPACK name zuncsd.

3  Description

The m by m unitary matrix X is partitioned as
X= X11 X12 X21 X22  
where X11 is a p by q submatrix and the dimensions of the other submatrices X12, X21 and X22 are such that X remains m by m.
The CS decomposition of X is X= UΣpVT  where U, V and Σp are m by m matrices, such that
U= U1 0 0 U2  
is a unitary matrix containing the p by p unitary matrix U1 and the m-p by m-p unitary matrix U2;
V= V1 0 0 V2  
is a unitary matrix containing the q by q unitary matrix V1 and the m-q by m-q unitary matrix V2; and
Σp= I11 0 0 0 C 0 0 -S 0 0 0 -I12 0 0 I22 0 0 S C 0 0 I21 0 0  
contains the r by r non-negative diagonal submatrices C and S satisfying C2+S2=I, where r=minp,m-p,q,m-q and the top left partition is p by q.
The identity matrix I11 is of order minp,q-r and vanishes if minp,q=r.
The identity matrix I12 is of order minp,m-q-r and vanishes if minp,m-q=r.
The identity matrix I21 is of order minm-p,q-r and vanishes if minm-p,q=r.
The identity matrix I22 is of order minm-p,m-q-r and vanishes if minm-p,m-q=r.
In each of the four cases r=p,q,m-p,m-q at least two of the identity matrices vanish.
The indicated zeros represent augmentations by additional rows or columns (but not both) to the square diagonal matrices formed by Iij and C or S.
Σp does not need to be stored in full; it is sufficient to return only the values θi for i=1,2,,r where Cii=cosθi and Sii=sinθi.
The algorithm used to perform the complete CS decomposition is described fully in Sutton (2009) including discussions of the stability and accuracy of the algorithm.

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 (2012) Matrix Computations (4th Edition) Johns Hopkins University Press, Baltimore
Sutton B D (2009) Computing the complete CS decomposition Numerical Algorithms (Volume 50) 1017–1398 Springer US 33–65 http://dx.doi.org/10.1007/s11075-008-9215-6

5  Arguments

1:     JOBU1 – CHARACTER(1)Input
On entry:
  • if JOBU1='Y', U1 is computed;
  • otherwise, U1 is not computed.
2:     JOBU2 – CHARACTER(1)Input
On entry:
  • if JOBU2='Y', U2 is computed;
  • otherwise, U2 is not computed.
3:     JOBV1T – CHARACTER(1)Input
On entry:
  • if JOBV1T='Y', V1T is computed;
  • otherwise, V1T is not computed.
4:     JOBV2T – CHARACTER(1)Input
On entry:
  • if JOBV2T='Y', V2T is computed;
  • otherwise, V2T is not computed.
5:     TRANS – CHARACTER(1)Input
On entry:
  • if TRANS='T', X, U1, U2, V1T and V2T are stored in row-major order;
  • otherwise, X, U1, U2, V1T and V2T are stored in column-major order.
6:     SIGNS – CHARACTER(1)Input
On entry:
  • if SIGNS='O', the lower-left block is made nonpositive (the other convention);
  • otherwise, the upper-right block is made nonpositive (the default convention).
7:     M – INTEGERInput
On entry: m, the number of rows and columns in the unitary matrix X.
Constraint: M0.
8:     P – INTEGERInput
On entry: p, the number of rows in X11 and X12.
Constraint: 0PM.
9:     Q – INTEGERInput
On entry: q, the number of columns in X11 and X21.
Constraint: 0QM.
10:   X11LDX11* – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array X11 must be at least max1,P if TRANS='T', and at least Q otherwise.
On entry: the upper left partition of the unitary matrix X whose CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
11:   LDX11 – INTEGERInput
On entry: the first dimension of the array X11 as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraints:
  • if TRANS='T', LDX11max1,Q;
  • otherwise LDX11max1,P.
12:   X12LDX12* – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array X12 must be at least max1,P if TRANS='T', and at least M-Q otherwise.
On entry: the upper right partition of the unitary matrix X whose CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
13:   LDX12 – INTEGERInput
On entry: the first dimension of the array X12 as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraints:
  • if TRANS='T', LDX12max1,M-Q;
  • otherwise LDX12max1,P.
14:   X21LDX21* – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array X21 must be at least max1,M-P if TRANS='T', and at least Q otherwise.
On entry: the lower left partition of the unitary matrix X whose CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
15:   LDX21 – INTEGERInput
On entry: the first dimension of the array X21 as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraints:
  • if TRANS='T', LDX21max1,Q;
  • otherwise LDX21max1,M-P.
16:   X22LDX22* – COMPLEX (KIND=nag_wp) arrayInput/Output
Note: the second dimension of the array X22 must be at least max1,M-P if TRANS='T', and at least M-Q otherwise.
On entry: the lower right partition of the unitary matrix X CSD is desired.
On exit: contains details of the unitary matrix used in a simultaneous bidiagonalization process.
17:   LDX22 – INTEGERInput
On entry: the first dimension of the array X22 as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraints:
  • if TRANS='T', LDX22max1,M-Q;
  • otherwise LDX22max1,M-P.
18:   THETAminP,M-P,Q,M-Q – REAL (KIND=nag_wp) arrayOutput
On exit: the values θi for i=1,2,,r where r=minp,m-p,q,m-q. The diagonal submatrices C and S of Σp are constructed from these values as
  • C= diagcos THETA1 ,,cos THETAr  and
  • S= diagsin THETA1 ,,sin THETAr .
19:   U1LDU1* – COMPLEX (KIND=nag_wp) arrayOutput
Note: the second dimension of the array U1 must be at least max1,P if JOBU1='Y', and at least 1 otherwise.
On exit: if JOBU1='Y', U1 contains the p by p unitary matrix U1.
20:   LDU1 – INTEGERInput
On entry: the first dimension of the array U1 as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraint: if JOBU1='Y', LDU1max1,P.
21:   U2LDU2* – COMPLEX (KIND=nag_wp) arrayOutput
Note: the second dimension of the array U2 must be at least max1,M-P if JOBU2='Y', and at least 1 otherwise.
On exit: if JOBU2='Y', U2 contains the m-p by m-p unitary matrix U2.
22:   LDU2 – INTEGERInput
On entry: the first dimension of the array U2 as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraint: if JOBU2='Y', LDU2max1,M-P.
23:   V1TLDV1T* – COMPLEX (KIND=nag_wp) arrayOutput
Note: the second dimension of the array V1T must be at least max1,Q if JOBV1T='Y', and at least 1 otherwise.
On exit: if JOBV1T='Y', V1T contains the q by q unitary matrix V1H.
24:   LDV1T – INTEGERInput
On entry: the first dimension of the array V1T as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraint: if JOBV1T='Y', LDV1Tmax1,Q.
25:   V2TLDV2T* – COMPLEX (KIND=nag_wp) arrayOutput
Note: the second dimension of the array V2T must be at least max1,M-Q if JOBV2T='Y', and at least 1 otherwise.
On exit: if JOBV2T='Y', V2T contains the m-q by m-q unitary matrix V2H.
26:   LDV2T – INTEGERInput
On entry: the first dimension of the array V2T as declared in the (sub)program from which F08RNF (ZUNCSD) is called.
Constraint: if JOBV2T='Y', LDV2Tmax1,M-Q.
27:   WORKmax1,LWORK – COMPLEX (KIND=nag_wp) arrayWorkspace
On exit: if INFO=0, WORK1 returns the optimal LWORK.
If INFO>0 on exit, WORK2:r contains the values PHI1,PHIr-1 that, together with xrefarg1,xrefargr, define the matrix in intermediate bidiagonal-block form remaining after nonconvergence. INFO specifies the number of nonzero PHI's.
28:   LWORK – INTEGERInput
On entry:
If LWORK=-1, a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued.
The minimum workspace required is max1,P + max1,M-P + max1,Q + max1,M-Q + max1,P,M-P,Q,M-Q + 1 ; the optimal amount of workspace depends on internal block sizes and the relative problem dimensions.
Constraint:
LWORK=-1 or LWORK max1,P + max1,M-P + max1,Q + max1,M-Q + max1,P,M-P,Q,M-Q + 1 .
29:   RWORKmax1,LRWORK – REAL (KIND=nag_wp) arrayWorkspace
30:   LRWORK – INTEGERInput
On entry:
If LRWORK=-1, a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LRWORK is issued. Otherwise the required workspace is 5 × max1,Q-1 + 4 × max1,Q + max1,8×Q +1  which equates to 11 for Q=0, 18 for Q=1 and 17×Q-4 when Q>1.
Constraint:
LRWORK=-1 ​ or ​ LRWORK5×max1,Q-1+4×max1,Q+max1,8×Q+1.
31:   IWORKM-minP,M-P,Q,M-Q – INTEGER arrayWorkspace
32:   INFO – INTEGEROutput
On exit: INFO=0 unless the routine detects an error (see Section 6).

6  Error Indicators and Warnings

INFO<0
If INFO=-i, argument i had an illegal value. An explanatory message is output, and execution of the program is terminated.
INFO>0
The Jacobi-type procedure failed to converge during an internal reduction to bidiagonal-block form. The process requires convergence to minP,M-P,Q,M-Q values, the value of INFO gives the number of converged values.

7  Accuracy

The computed CS decomposition is nearly the exact CS decomposition for the nearby matrix X+E , where
E2 = Oε ,  
and ε  is the machine precision.

8  Parallelism and Performance

F08RNF (ZUNCSD) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
F08RNF (ZUNCSD) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9  Further Comments

The total number of floating-point operations required to perform the full CS decomposition is approximately 2m3.
The real analogue of this routine is F08RAF (DORCSD).

10  Example

This example finds the full CS decomposition of a unitary 6 by 6 matrix X (see Section 10.2) partitioned so that the top left block is 2 by 4.
The decomposition is performed both on submatrices of the unitary matrix X and on separated partition matrices. Code is also provided to perform a recombining check if required.

10.1  Program Text

Program Text (f08rnfe.f90)

10.2  Program Data

Program Data (f08rnfe.d)

10.3  Program Results

Program Results (f08rnfe.r)


F08RNF (ZUNCSD) (PDF version)
F08 Chapter Contents
F08 Chapter Introduction
NAG Library Manual

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