# ?bdsqr

Computes the singular value decomposition of a general matrix that has been reduced to bidiagonal form.

## Syntax

FORTRAN 77:

call sbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)

call dbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)

call cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)

call zbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)

Fortran 95:

call rbdsqr(d, e [,vt] [,u] [,c] [,uplo] [,info])

call bdsqr(d, e [,vt] [,u] [,c] [,uplo] [,info])

## Description

This routine computes the singular values and, optionally, the right and/or left singular vectors from the Singular Value Decomposition (SVD) of a real n-by-n (upper or lower) bidiagonal matrix B using the implicit zero-shift QR algorithm. The SVD of B has the form `B = Q*S*PH` where S is the diagonal matrix of singular values, Q is an orthogonal matrix of left singular vectors, and P is an orthogonal matrix of right singular vectors. If left singular vectors are requested, this subroutine actually returns U *Q instead of Q, and, if right singular vectors are requested, this subroutine returns PH *VT instead of PH, for given real/complex input matrices U and VT. When U and VT are the orthogonal/unitary matrices that reduce a general matrix A to bidiagonal form: `A = U*B*VT`, as computed by ?gebrd, then

`A = (U*Q)*S*(PH*VT)`

is the SVD of A. Optionally, the subroutine may also compute QH *C for a given real/complex input matrix C.

See also ?lasq1, ?lasq2, ?lasq3, ?lasq4, ?lasq5, ?lasq6 used by this routine.

## Input Parameters

uplo

  CHARACTER*1. Must be 'U' or 'L'.

  If `uplo = 'U'`, B is an upper bidiagonal matrix.

  If `uplo = 'L'`, B is a lower bidiagonal matrix.

n

  INTEGER. The order of the matrix B (`n ≥ 0`).

ncvt

  INTEGER. The number of columns of the matrix VT, that is, the number of right singular vectors (`ncvt ≥ 0`).

  Set `ncvt = 0` if no right singular vectors are required.

nru

  INTEGER. The number of rows in U, that is, the number of left singular vectors (`nru ≥ 0`).

  Set `nru = 0` if no left singular vectors are required.

ncc

  INTEGER. The number of columns in the matrix C used for computing the product QH*C (`ncc ≥ 0`). Set `ncc = 0` if no matrix C is supplied.

d, e, work

  REAL for single-precision flavors

  DOUBLE PRECISION for double-precision flavors.

  Arrays:

  d`(*)` contains the diagonal elements of B.

  The dimension of d must be at least `max(1, n)`.

  e`(*)` contains the `(n-1)` off-diagonal elements of B.

  The dimension of e must be at least `max(1, n)`. `e(n)` is used for workspace.

  `work(*)` is a workspace array.

  The dimension of work must be at least `max(1, 2*n)` if `ncvt = nru` = `ncc = 0`;

  `max(1, 4*(n-1))` otherwise.

vt, u, c

  REAL for sbdsqr

  DOUBLE PRECISION for dbdsqr

  COMPLEX for cbdsqr

  DOUBLE COMPLEX for zbdsqr.

  Arrays:

  `vt(ldvt,*)` contains an n-by-ncvt matrix VT.

  The second dimension of vt must be at least `max(1, ncvt)`.

  vt is not referenced if `ncvt = 0`.

  `u(ldu,*)` contains an nru by n unit matrix U.

  The second dimension of u must be at least `max(1, n)`.

  u is not referenced if `nru = 0`.

  c(ldc,*) contains the matrix C for computing the product `QH*C`.

  The second dimension of c must be at least `max(1, ncc)`. The array is not referenced if `ncc = 0`.

ldvt

  INTEGER. The first dimension of vt. Constraints:

  `ldvt ≥ max(1, n)` if `ncvt > 0`;

  `ldvt ≥ 1` if `ncvt = 0`.

ldu

  INTEGER. The first dimension of u. Constraint:

  `ldu ≥ max(1, nru)`.

ldc

  INTEGER. The first dimension of c. Constraints:

  `ldc ≥ max(1, n)` if `ncc > 0;ldc ≥ 1` otherwise.

## Output Parameters

d

  On exit, if `info = 0`, overwritten by the singular values in decreasing order (see info).

e

  On exit, if `info = 0`, e is destroyed. See also info below.

c

  Overwritten by the product QH*C.

vt

  On exit, this array is overwritten by PH *VT.

u

  On exit, this array is overwritten by U *Q .

info

  INTEGER.

  If `info = 0`, the execution is successful.

  If `info = -i`, the i-th parameter had an illegal value.

  If `info = i`, the algorithm failed to converge; i specifies how many off-diagonals did not converge.

  In this case, d and e contain on exit the diagonal and off-diagonal elements, respectively, of a bidiagonal matrix orthogonally equivalent to B.

## Fortran 95 Interface Notes

Routines in Fortran 95 interface have fewer arguments in the calling sequence than their FORTRAN 77 counterparts. For general conventions applied to skip redundant or restorable arguments, see Fortran 95 Interface Conventions.

Specific details for the routine bdsqr interface are the following:

d

  Holds the vector of length (n).

e

  Holds the vector of length (n).

vt

  Holds the matrix VT of size `(n, ncvt)`.

u

  Holds the matrix U of size (nru,n).

c

  Holds the matrix C of size (n,ncc).

uplo

  Must be 'U' or 'L'. The default value is 'U'.

ncvt

  If argument vt is present, then ncvt is equal to the number of columns in matrix VT; otherwise, ncvt is set to zero.

nru

  If argument u is present, then nru is equal to the number of rows in matrix U; otherwise, nru is set to zero.

ncc

  If argument c is present, then ncc is equal to the number of columns in matrix C; otherwise, ncc is set to zero.

Note that two variants of Fortran 95 interface for bdsqr routine are needed because of an ambiguous choice between real and complex cases appear when vt, u, and c are omitted. Thus, the name rbdsqr is used in real cases (single or double precision), and the name bdsqr is used in complex cases (single or double precision).

## Application Notes

Each singular value and singular vector is computed to high relative accuracy. However, the reduction to bidiagonal form (prior to calling the routine) may decrease the relative accuracy in the small singular values of the original matrix if its singular values vary widely in magnitude.

If si is an exact singular value of B, and si is the corresponding computed value, then

`|si - σi| ≤ p*(m,n)*ε*σi`

where p(m, n) is a modestly increasing function of m and n, and ε is the machine precision.

If only singular values are computed, they are computed more accurately than when some singular vectors are also computed (that is, the function p(m, n) is smaller).

If ui is the corresponding exact left singular vector of B, and wi is the corresponding computed left singular vector, then the angle θ(ui, wi) between them is bounded as follows:

`θ(ui, wi) ≤ p(m,n)*ε / min i≠j(|σi - σj|/|σi + σj|)`.

Here `mini≠j(|σi - σj|/|σi + σj|)` is the relative gap between σi and the other singular values. A similar error bound holds for the right singular vectors.

The total number of real floating-point operations is roughly proportional to n2 if only the singular values are computed. About `6n2*nru` additional operations (`12n2*nru` for complex flavors) are required to compute the left singular vectors and about `6n2*ncvt` operations (`12n2*ncvt` for complex flavors) to compute the right singular vectors.

Parent topic: Singular Value Decomposition