# DOUBLE PRECISION routines for bidiagonal matrix

## dbdsdc

```USAGE:
u, vt, q, iq, info, d, e = NumRu::Lapack.dbdsdc( uplo, compq, d, e, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO )

*  Purpose
*  =======
*
*  DBDSDC computes the singular value decomposition (SVD) of a real
*  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
*  using a divide and conquer method, where S is a diagonal matrix
*  with non-negative diagonal elements (the singular values of B), and
*  U and VT are orthogonal matrices of left and right singular vectors,
*  respectively. DBDSDC can be used to compute all singular values,
*  and optionally, singular vectors or singular vectors in compact form.
*
*  This code makes very mild assumptions about floating point
*  arithmetic. It will work on machines with a guard digit in
*  add/subtract, or on those binary machines without guard digits
*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
*  It could conceivably fail on hexadecimal or decimal machines
*  without guard digits, but we know of none.  See DLASD3 for details.
*
*  The code currently calls DLASDQ if singular values only are desired.
*  However, it can be slightly modified to compute singular values
*  using the divide and conquer method.
*

*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  B is upper bidiagonal.
*          = 'L':  B is lower bidiagonal.
*
*  COMPQ   (input) CHARACTER*1
*          Specifies whether singular vectors are to be computed
*          as follows:
*          = 'N':  Compute singular values only;
*          = 'P':  Compute singular values and compute singular
*                  vectors in compact form;
*          = 'I':  Compute singular values and singular vectors.
*
*  N       (input) INTEGER
*          The order of the matrix B.  N >= 0.
*
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the n diagonal elements of the bidiagonal matrix B.
*          On exit, if INFO=0, the singular values of B.
*
*  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
*          On entry, the elements of E contain the offdiagonal
*          elements of the bidiagonal matrix whose SVD is desired.
*          On exit, E has been destroyed.
*
*  U       (output) DOUBLE PRECISION array, dimension (LDU,N)
*          If  COMPQ = 'I', then:
*             On exit, if INFO = 0, U contains the left singular vectors
*             of the bidiagonal matrix.
*          For other values of COMPQ, U is not referenced.
*
*  LDU     (input) INTEGER
*          The leading dimension of the array U.  LDU >= 1.
*          If singular vectors are desired, then LDU >= max( 1, N ).
*
*  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
*          If  COMPQ = 'I', then:
*             On exit, if INFO = 0, VT' contains the right singular
*             vectors of the bidiagonal matrix.
*          For other values of COMPQ, VT is not referenced.
*
*  LDVT    (input) INTEGER
*          The leading dimension of the array VT.  LDVT >= 1.
*          If singular vectors are desired, then LDVT >= max( 1, N ).
*
*  Q       (output) DOUBLE PRECISION array, dimension (LDQ)
*          If  COMPQ = 'P', then:
*             On exit, if INFO = 0, Q and IQ contain the left
*             and right singular vectors in a compact form,
*             requiring O(N log N) space instead of 2*N**2.
*             In particular, Q contains all the DOUBLE PRECISION data in
*             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
*             words of memory, where SMLSIZ is returned by ILAENV and
*             is equal to the maximum size of the subproblems at the
*             bottom of the computation tree (usually about 25).
*          For other values of COMPQ, Q is not referenced.
*
*  IQ      (output) INTEGER array, dimension (LDIQ)
*          If  COMPQ = 'P', then:
*             On exit, if INFO = 0, Q and IQ contain the left
*             and right singular vectors in a compact form,
*             requiring O(N log N) space instead of 2*N**2.
*             In particular, IQ contains all INTEGER data in
*             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
*             words of memory, where SMLSIZ is returned by ILAENV and
*             is equal to the maximum size of the subproblems at the
*             bottom of the computation tree (usually about 25).
*          For other values of COMPQ, IQ is not referenced.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*          If COMPQ = 'N' then LWORK >= (4 * N).
*          If COMPQ = 'P' then LWORK >= (6 * N).
*          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
*
*  IWORK   (workspace) INTEGER array, dimension (8*N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit.
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          > 0:  The algorithm failed to compute a singular value.
*                The update process of divide and conquer failed.
*

*  Further Details
*  ===============
*
*  Based on contributions by
*     Ming Gu and Huan Ren, Computer Science Division, University of
*     California at Berkeley, USA
*
*  =====================================================================
*  Changed dimension statement in comment describing E from (N) to
*  (N-1).  Sven, 17 Feb 05.
*  =====================================================================
*

```
go to the page top

## dbdsqr

```USAGE:
info, d, e, vt, u, c = NumRu::Lapack.dbdsqr( uplo, nru, d, e, vt, u, c, [:usage => usage, :help => help])

FORTRAN MANUAL
SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO )

*  Purpose
*  =======
*
*  DBDSQR 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 * P**T
*
*  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 P**T*VT instead of
*  P**T, for given real input matrices U and VT.  When U and VT are the
*  orthogonal matrices that reduce a general matrix A to bidiagonal
*  form:  A = U*B*VT, as computed by DGEBRD, then
*
*     A = (U*Q) * S * (P**T*VT)
*
*  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
*  for a given real input matrix C.
*
*  See "Computing  Small Singular Values of Bidiagonal Matrices With
*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
*  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
*  no. 5, pp. 873-912, Sept 1990) and
*  "Accurate singular values and differential qd algorithms," by
*  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
*  Department, University of California at Berkeley, July 1992
*  for a detailed description of the algorithm.
*

*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  B is upper bidiagonal;
*          = 'L':  B is lower bidiagonal.
*
*  N       (input) INTEGER
*          The order of the matrix B.  N >= 0.
*
*  NCVT    (input) INTEGER
*          The number of columns of the matrix VT. NCVT >= 0.
*
*  NRU     (input) INTEGER
*          The number of rows of the matrix U. NRU >= 0.
*
*  NCC     (input) INTEGER
*          The number of columns of the matrix C. NCC >= 0.
*
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the n diagonal elements of the bidiagonal matrix B.
*          On exit, if INFO=0, the singular values of B in decreasing
*          order.
*
*  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
*          On entry, the N-1 offdiagonal elements of the bidiagonal
*          matrix B.
*          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
*          will contain the diagonal and superdiagonal elements of a
*          bidiagonal matrix orthogonally equivalent to the one given
*          as input.
*
*  VT      (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
*          On entry, an N-by-NCVT matrix VT.
*          On exit, VT is overwritten by P**T * VT.
*          Not referenced if NCVT = 0.
*
*  LDVT    (input) INTEGER
*          The leading dimension of the array VT.
*          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
*
*  U       (input/output) DOUBLE PRECISION array, dimension (LDU, N)
*          On entry, an NRU-by-N matrix U.
*          On exit, U is overwritten by U * Q.
*          Not referenced if NRU = 0.
*
*  LDU     (input) INTEGER
*          The leading dimension of the array U.  LDU >= max(1,NRU).
*
*  C       (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
*          On entry, an N-by-NCC matrix C.
*          On exit, C is overwritten by Q**T * C.
*          Not referenced if NCC = 0.
*
*  LDC     (input) INTEGER
*          The leading dimension of the array C.
*          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
*
*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  If INFO = -i, the i-th argument had an illegal value
*          > 0:
*             if NCVT = NRU = NCC = 0,
*                = 1, a split was marked by a positive value in E
*                = 2, current block of Z not diagonalized after 30*N
*                     iterations (in inner while loop)
*                = 3, termination criterion of outer while loop not met
*                     (program created more than N unreduced blocks)
*             else NCVT = NRU = NCC = 0,
*                   the algorithm did not converge; D and E contain the
*                   elements of a bidiagonal matrix which is orthogonally
*                   similar to the input matrix B;  if INFO = i, i
*                   elements of E have not converged to zero.
*
*  Internal Parameters
*  ===================
*
*  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
*          TOLMUL controls the convergence criterion of the QR loop.
*          If it is positive, TOLMUL*EPS is the desired relative
*             precision in the computed singular values.
*          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
*             desired absolute accuracy in the computed singular
*             values (corresponds to relative accuracy
*             abs(TOLMUL*EPS) in the largest singular value.
*          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
*             between 10 (for fast convergence) and .1/EPS
*             (for there to be some accuracy in the results).
*          Default is to lose at either one eighth or 2 of the
*             available decimal digits in each computed singular value
*             (whichever is smaller).
*
*  MAXITR  INTEGER, default = 6
*          MAXITR controls the maximum number of passes of the
*          algorithm through its inner loop. The algorithms stops
*          (and so fails to converge) if the number of passes
*          through the inner loop exceeds MAXITR*N**2.
*

*  =====================================================================
*

```
go to the page top
back to matrix types
back to data types