!>
!> DBBCSD computes the CS decomposition of an orthogonal matrix in
!> bidiagonal-block form,
!>
!>
!>     [ B11 | B12 0  0 ]
!>     [  0  |  0 -I  0 ]
!> X = [----------------]
!>     [ B21 | B22 0  0 ]
!>     [  0  |  0  0  I ]
!>
!>                               [  C | -S  0  0 ]
!>                   [ U1 |    ] [  0 |  0 -I  0 ] [ V1 |    ]**T
!>                 = [---------] [---------------] [---------]   .
!>                   [    | U2 ] [  S |  C  0  0 ] [    | V2 ]
!>                               [  0 |  0  0  I ]
!>
!> X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
!> than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
!> transposed and/or permuted. This can be done in constant time using
!> the TRANS and SIGNS options. See DORCSD for details.)
!>
!> The bidiagonal matrices B11, B12, B21, and B22 are represented
!> implicitly by angles THETA(1:Q) and PHI(1:Q-1).
!>
!> The orthogonal matrices U1, U2, V1T, and V2T are input/output.
!> The input matrices are pre- or post-multiplied by the appropriate
!> singular vector matrices.
!> 

JOBU1

!>          JOBU1 is CHARACTER
!>          = 'Y':      U1 is updated;
!>          otherwise:  U1 is not updated.
!> 

JOBU2

!>          JOBU2 is CHARACTER
!>          = 'Y':      U2 is updated;
!>          otherwise:  U2 is not updated.
!> 

JOBV1T

!>          JOBV1T is CHARACTER
!>          = 'Y':      V1T is updated;
!>          otherwise:  V1T is not updated.
!> 

JOBV2T

!>          JOBV2T is CHARACTER
!>          = 'Y':      V2T is updated;
!>          otherwise:  V2T is not updated.
!> 

TRANS

!>          TRANS is CHARACTER
!>          = '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.
!> 

M

!>          M is INTEGER
!>          The number of rows and columns in X, the orthogonal matrix in
!>          bidiagonal-block form.
!> 

P

!>          P is INTEGER
!>          The number of rows in the top-left block of X. 0 <= P <= M.
!> 

Q

!>          Q is INTEGER
!>          The number of columns in the top-left block of X.
!>          0 <= Q <= MIN(P,M-P,M-Q).
!> 

THETA

!>          THETA is DOUBLE PRECISION array, dimension (Q)
!>          On entry, the angles THETA(1),...,THETA(Q) that, along with
!>          PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
!>          form. On exit, the angles whose cosines and sines define the
!>          diagonal blocks in the CS decomposition.
!> 

PHI

!>          PHI is DOUBLE PRECISION array, dimension (Q-1)
!>          The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
!>          THETA(Q), define the matrix in bidiagonal-block form.
!> 

U1

!>          U1 is DOUBLE PRECISION array, dimension (LDU1,P)
!>          On entry, a P-by-P matrix. On exit, U1 is postmultiplied
!>          by the left singular vector matrix common to [ B11 ; 0 ] and
!>          [ B12 0 0 ; 0 -I 0 0 ].
!> 

LDU1

!>          LDU1 is INTEGER
!>          The leading dimension of the array U1, LDU1 >= MAX(1,P).
!> 

U2

!>          U2 is DOUBLE PRECISION array, dimension (LDU2,M-P)
!>          On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is
!>          postmultiplied by the left singular vector matrix common to
!>          [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
!> 

LDU2

!>          LDU2 is INTEGER
!>          The leading dimension of the array U2, LDU2 >= MAX(1,M-P).
!> 

V1T

!>          V1T is DOUBLE PRECISION array, dimension (LDV1T,Q)
!>          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
!>          by the transpose of the right singular vector
!>          matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
!> 

LDV1T

!>          LDV1T is INTEGER
!>          The leading dimension of the array V1T, LDV1T >= MAX(1,Q).
!> 

V2T

!>          V2T is DOUBLE PRECISION array, dimension (LDV2T,M-Q)
!>          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
!>          premultiplied by the transpose of the right
!>          singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
!>          [ B22 0 0 ; 0 0 I ].
!> 

LDV2T

!>          LDV2T is INTEGER
!>          The leading dimension of the array V2T, LDV2T >= MAX(1,M-Q).
!> 

B11D

!>          B11D is DOUBLE PRECISION array, dimension (Q)
!>          When DBBCSD converges, B11D contains the cosines of THETA(1),
!>          ..., THETA(Q). If DBBCSD fails to converge, then B11D
!>          contains the diagonal of the partially reduced top-left
!>          block.
!> 

B11E

!>          B11E is DOUBLE PRECISION array, dimension (Q-1)
!>          When DBBCSD converges, B11E contains zeros. If DBBCSD fails
!>          to converge, then B11E contains the superdiagonal of the
!>          partially reduced top-left block.
!> 

B12D

!>          B12D is DOUBLE PRECISION array, dimension (Q)
!>          When DBBCSD converges, B12D contains the negative sines of
!>          THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then
!>          B12D contains the diagonal of the partially reduced top-right
!>          block.
!> 

B12E

!>          B12E is DOUBLE PRECISION array, dimension (Q-1)
!>          When DBBCSD converges, B12E contains zeros. If DBBCSD fails
!>          to converge, then B12E contains the subdiagonal of the
!>          partially reduced top-right block.
!> 

B21D

!>          B21D is DOUBLE PRECISION  array, dimension (Q)
!>          When DBBCSD converges, B21D contains the negative sines of
!>          THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then
!>          B21D contains the diagonal of the partially reduced bottom-left
!>          block.
!> 

B21E

!>          B21E is DOUBLE PRECISION  array, dimension (Q-1)
!>          When DBBCSD converges, B21E contains zeros. If DBBCSD fails
!>          to converge, then B21E contains the subdiagonal of the
!>          partially reduced bottom-left block.
!> 

B22D

!>          B22D is DOUBLE PRECISION  array, dimension (Q)
!>          When DBBCSD converges, B22D contains the negative sines of
!>          THETA(1), ..., THETA(Q). If DBBCSD fails to converge, then
!>          B22D contains the diagonal of the partially reduced bottom-right
!>          block.
!> 

B22E

!>          B22E is DOUBLE PRECISION  array, dimension (Q-1)
!>          When DBBCSD converges, B22E contains zeros. If DBBCSD fails
!>          to converge, then B22E contains the subdiagonal of the
!>          partially reduced bottom-right block.
!> 

WORK

!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!> 

LWORK

!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= MAX(1,8*Q).
!>
!>          If LWORK = -1, then 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 by XERBLA.
!> 

INFO

!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if DBBCSD did not converge, INFO specifies the number
!>                of nonzero entries in PHI, and B11D, B11E, etc.,
!>                contain the partially reduced matrix.
!> 

!>  TOLMUL  DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
!>          TOLMUL controls the convergence criterion of the QR loop.
!>          Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
!>          are within TOLMUL*EPS of either bound.
!> 

[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.