nag_dorcsd (f08rac) (PDF version)
f08 Chapter Contents
f08 Chapter Introduction
NAG Library Manual

NAG Library Function Document

nag_dorcsd (f08rac)

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

nag_dorcsd (f08rac) computes the CS decomposition of a real m by m orthogonal matrix X, partitioned into a 2 by 2 array of submatrices.

2  Specification

#include <nag.h>
#include <nagf08.h>
void  nag_dorcsd (Nag_OrderType order, Nag_ComputeUType jobu1, Nag_ComputeUType jobu2, Nag_ComputeVTType jobv1t, Nag_ComputeVTType jobv2t, Nag_SignsType signs, Integer m, Integer p, Integer q, double x11[], Integer pdx11, double x12[], Integer pdx12, double x21[], Integer pdx21, double x22[], Integer pdx22, double theta[], double u1[], Integer pdu1, double u2[], Integer pdu2, double v1t[], Integer pdv1t, double v2t[], Integer pdv2t, NagError *fail)

3  Description

The m by m orthogonal 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 an orthogonal matrix containing the p by p orthogonal matrix U1 and the m-p by m-p orthogonal matrix U2;
V= V1 0 0 V2
is an orthogonal matrix containing the q by q orthogonal matrix V1 and the m-q by m-q orthogonal 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:     orderNag_OrderTypeInput
On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by order=Nag_RowMajor. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.
Constraint: order=Nag_RowMajor or Nag_ColMajor.
2:     jobu1Nag_ComputeUTypeInput
On entry:
  • if jobu1=Nag_AllU, U1 is computed;
  • if jobu1=Nag_NotU, U1 is not computed.
Constraint: jobu1=Nag_AllU or Nag_NotU.
3:     jobu2Nag_ComputeUTypeInput
On entry:
  • if jobu2=Nag_AllU, U2 is computed;
  • if jobu2=Nag_NotU, U2 is not computed.
Constraint: jobu2=Nag_AllU or Nag_NotU.
4:     jobv1tNag_ComputeVTTypeInput
On entry:
  • if jobv1t=Nag_AllVT, V1T is computed;
  • if jobv1t=Nag_NotVT, V1T is not computed.
Constraint: jobv1t=Nag_AllVT or Nag_NotVT.
5:     jobv2tNag_ComputeVTTypeInput
On entry:
  • if jobv2t=Nag_AllVT, V2T is computed;
  • if jobv2t=Nag_NotVT, V2T is not computed.
Constraint: jobv2t=Nag_AllVT or Nag_NotVT.
6:     signsNag_SignsTypeInput
On entry:
  • if signs=Nag_LowerMinus, the lower-left block is made nonpositive (the other convention);
  • if signs=Nag_UpperMinus, the upper-right block is made nonpositive (the default convention).
Constraint: signs=Nag_LowerMinus or Nag_UpperMinus.
7:     mIntegerInput
On entry: m, the number of rows and columns in the orthogonal matrix X.
Constraint: m0.
8:     pIntegerInput
On entry: p, the number of rows in X11 and X12.
Constraint: 0pm.
9:     qIntegerInput
On entry: q, the number of columns in X11 and X21.
Constraint: 0qm.
10:   x11[dim]doubleInput/Output
Note: the dimension, dim, of the array x11 must be at least
  • max1,pdx11×p when order=Nag_RowMajor;
  • max1,pdx11×q when order=Nag_ColMajor.
The i,jth element of the matrix is stored in
  • x11[j-1×pdx11+i-1] when order=Nag_ColMajor;
  • x11[i-1×pdx11+j-1] when order=Nag_RowMajor.
On entry: the upper left partition of the orthogonal matrix X whose CSD is desired.
On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
11:   pdx11IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x11.
Constraints:
  • if order=Nag_RowMajor, pdx11max1,q;
  • if order=Nag_ColMajor, pdx11max1,p.
12:   x12[dim]doubleInput/Output
Note: the dimension, dim, of the array x12 must be at least
  • max1,pdx12×p when order=Nag_RowMajor;
  • max1,pdx12×(m-q) when order=Nag_ColMajor.
The i,jth element of the matrix is stored in
  • x12[j-1×pdx12+i-1] when order=Nag_ColMajor;
  • x12[i-1×pdx12+j-1] when order=Nag_RowMajor.
On entry: the upper right partition of the orthogonal matrix X whose CSD is desired.
On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
13:   pdx12IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x12.
Constraints:
  • if order=Nag_RowMajor, pdx12max1,m-q;
  • if order=Nag_ColMajor, pdx12max1,p.
14:   x21[dim]doubleInput/Output
Note: the dimension, dim, of the array x21 must be at least
  • max1,pdx21×(m-p) when order=Nag_RowMajor;
  • max1,pdx21×q when order=Nag_ColMajor.
The i,jth element of the matrix is stored in
  • x21[j-1×pdx21+i-1] when order=Nag_ColMajor;
  • x21[i-1×pdx21+j-1] when order=Nag_RowMajor.
On entry: the lower left partition of the orthogonal matrix X whose CSD is desired.
On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
15:   pdx21IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x21.
Constraints:
  • if order=Nag_RowMajor, pdx21max1,q;
  • if order=Nag_ColMajor, pdx21max1,m-p.
16:   x22[dim]doubleInput/Output
Note: the dimension, dim, of the array x22 must be at least
  • max1,pdx22×(m-p) when order=Nag_RowMajor;
  • max1,pdx22×(m-q) when order=Nag_ColMajor.
The i,jth element of the matrix is stored in
  • x22[j-1×pdx22+i-1] when order=Nag_ColMajor;
  • x22[i-1×pdx22+j-1] when order=Nag_RowMajor.
On entry: the lower right partition of the orthogonal matrix X CSD is desired.
On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
17:   pdx22IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array x22.
Constraints:
  • if order=Nag_RowMajor, pdx22max1,m-q;
  • if order=Nag_ColMajor, pdx22max1,m-p.
18:   theta[dim]doubleOutput
Note: the dimension, dim, of the array theta must be at least minp,m-p,q,m-q.
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 theta[0] ,,cos theta[r-1]  and
  • S= diagsin theta[0] ,,sin theta[r-1] .
19:   u1[dim]doubleOutput
Note: the dimension, dim, of the array u1 must be at least
  • max1,pdu1×p when jobu1=Nag_AllU;
  • otherwise u1 may be NULL.
The i,jth element of the matrix is stored in
  • u1[j-1×pdu1+i-1] when order=Nag_ColMajor;
  • u1[i-1×pdu1+j-1] when order=Nag_RowMajor.
On exit: if jobu1=Nag_AllU, u1 contains the p by p orthogonal matrix U1.
20:   pdu1IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array u1.
Constraint: if jobu1=Nag_AllU, pdu1max1,p
21:   u2[dim]doubleOutput
Note: the dimension, dim, of the array u2 must be at least
  • max1,pdu2×(m-p) when jobu2=Nag_AllU;
  • otherwise u2 may be NULL.
The i,jth element of the matrix is stored in
  • u2[j-1×pdu2+i-1] when order=Nag_ColMajor;
  • u2[i-1×pdu2+j-1] when order=Nag_RowMajor.
On exit: if jobu2=Nag_AllU, u2 contains the m-p by m-p orthogonal matrix U2.
22:   pdu2IntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array u2.
Constraint: if jobu2=Nag_AllU, pdu2max1,m-p
23:   v1t[dim]doubleOutput
Note: the dimension, dim, of the array v1t must be at least
  • max1,pdv1t×q when jobv1t=Nag_AllVT;
  • otherwise v1t may be NULL.
The i,jth element of the matrix is stored in
  • v1t[j-1×pdv1t+i-1] when order=Nag_ColMajor;
  • v1t[i-1×pdv1t+j-1] when order=Nag_RowMajor.
On exit: if jobv1t=Nag_AllVT, v1t contains the q by q orthogonal matrix V1T.
24:   pdv1tIntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array v1t.
Constraint: if jobv1t=Nag_AllVT, pdv1tmax1,q
25:   v2t[dim]doubleOutput
Note: the dimension, dim, of the array v2t must be at least
  • max1,pdv2t×(m-q) when jobv2t=Nag_AllVT;
  • otherwise v2t may be NULL.
The i,jth element of the matrix is stored in
  • v2t[j-1×pdv2t+i-1] when order=Nag_ColMajor;
  • v2t[i-1×pdv2t+j-1] when order=Nag_RowMajor.
On exit: if jobv2t=Nag_AllVT, v2t contains the m-q by m-q orthogonal matrix V2T.
26:   pdv2tIntegerInput
On entry: the stride separating row or column elements (depending on the value of order) in the array v2t.
Constraint: if jobv2t=Nag_AllVT, pdv2tmax1,m-q
27:   failNagError *Input/Output
The NAG error argument (see Section 3.6 in the Essential Introduction).

6  Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
NE_BAD_PARAM
On entry, argument value had an illegal value.
NE_CONVERGENCE
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 fail.errnum gives the number of converged values.
NE_ENUM_INT_2
On entry, jobu1=value, pdu1=value and p=value.
Constraint: if jobu1=Nag_AllU, pdu1max1,p.
On entry, jobv1t=value, pdv1t=value and q=value.
Constraint: if jobv1t=Nag_AllVT, pdv1tmax1,q.
NE_ENUM_INT_3
On entry, jobu2=value, pdu2=value, m=value and p=value.
Constraint: if jobu2=Nag_AllU, pdu2max1,m-p.
On entry, jobv2t=value, pdv2t=value, m=value and q=value.
Constraint: if jobv2t=Nag_AllVT, pdv2tmax1,m-q.
On entry, pdx11=value, p=value, q=value and order=value.
Constraint: if order=Nag_RowMajor, pdx11max1,q;
if order=Nag_ColMajor, pdx11max1,p.
NE_ENUM_INT_4
On entry, pdx12=value, m=value, p=value, q=value and order=value.
Constraint: if order=Nag_RowMajor, pdx12max1,m-q;
if order=Nag_ColMajor, pdx12max1,p.
On entry, pdx21=value, m=value, p=value, q=value and order=value.
Constraint: if order=Nag_RowMajor, pdx21max1,q;
if order=Nag_ColMajor, pdx21max1,m-p.
On entry, pdx22=value, m=value, p=value, q=value and order=value.
Constraint: if order=Nag_RowMajor, pdx22max1,m-q;
if order=Nag_ColMajor, pdx22max1,m-p.
NE_INT
On entry, m=value.
Constraint: m0.
NE_INT_2
On entry, m=value and p=value.
Constraint: 0pm.
On entry, m=value and q=value.
Constraint: 0qm.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.

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

nag_dorcsd (f08rac) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
nag_dorcsd (f08rac) 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 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 complex analogue of this function is nag_zuncsd (f08rnc).

10  Example

This example finds the full CS decomposition of
X = -0.13484 0.52524 -0.20924 0.81373 0.67420 -0.52213 -0.38886 0.34874 0.26968 0.52757 -0.65782 -0.46499 0.67420 0.41615 0.61014 0.00000
partitioned in 2 by 2 blocks.
The decomposition is performed both on submatrices of the orthogonal matrix X and on separated partition matrices. Code is also provided to perform a recombining check if required.

10.1  Program Text

Program Text (f08race.c)

10.2  Program Data

Program Data (f08race.d)

10.3  Program Results

Program Results (f08race.r)


nag_dorcsd (f08rac) (PDF version)
f08 Chapter Contents
f08 Chapter Introduction
NAG Library Manual

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