dbdsdc.f

*> \brief \b DBDSDC
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DBDSDC + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
*                          WORK, IWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          COMPQ, UPLO
*       INTEGER            INFO, LDU, LDVT, N
*       ..
*       .. Array Arguments ..
*       INTEGER            IQ( * ), IWORK( * )
*       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
*      $                   VT( LDVT, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  B is upper bidiagonal.
*>          = 'L':  B is lower bidiagonal.
*> \endverbatim
*>
*> \param[in] COMPQ
*> \verbatim
*>          COMPQ is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix B.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] D
*> \verbatim
*>          D is 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.
*> \endverbatim
*>
*> \param[in,out] E
*> \verbatim
*>          E is 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.
*> \endverbatim
*>
*> \param[out] U
*> \verbatim
*>          U is 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.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U.  LDU >= 1.
*>          If singular vectors are desired, then LDU >= max( 1, N ).
*> \endverbatim
*>
*> \param[out] VT
*> \verbatim
*>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
*>          If  COMPQ = 'I', then:
*>             On exit, if INFO = 0, VT**T contains the right singular
*>             vectors of the bidiagonal matrix.
*>          For other values of COMPQ, VT is not referenced.
*> \endverbatim
*>
*> \param[in] LDVT
*> \verbatim
*>          LDVT is INTEGER
*>          The leading dimension of the array VT.  LDVT >= 1.
*>          If singular vectors are desired, then LDVT >= max( 1, N ).
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is 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.
*> \endverbatim
*>
*> \param[out] IQ
*> \verbatim
*>          IQ is 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.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is 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).
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (8*N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is 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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup bdsdc
*
*> \par Contributors:
*  ==================
*>
*>     Ming Gu and Huan Ren, Computer Science Division, University of
*>     California at Berkeley, USA
*>
*  =====================================================================
      SUBROUTINE dbdsdc( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
     $                   WORK, IWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          COMPQ, UPLO
      INTEGER            INFO, LDU, LDVT, N
*     ..
*     .. Array Arguments ..
      INTEGER            IQ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( ldu, * ),
     $                   vt( ldvt, * ), work( * )
*     ..
*
*  =====================================================================
*  Changed dimension statement in comment describing E from (N) to
*  (N-1).  Sven, 17 Feb 05.
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO
      parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
     $                   icompq, ierr, ii, is, iu, iuplo, ivt, j, k, kk,
     $                   mlvl, nm1, nsize, perm, poles, qstart, smlsiz,
     $                   smlszp, sqre, start, wstart, z
      DOUBLE PRECISION   CS, EPS, ORGNRM, P, R, SN
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           lsame, ilaenv, dlamch, dlanst
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlartg, dlascl, dlasd0, dlasda, dlasdq,
     $                   dlaset, dlasr, dswap, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, int, log, sign
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      info = 0
*
      iuplo = 0
      IF( lsame( uplo, 'U' ) )
     $   iuplo = 1
      IF( lsame( uplo, 'L' ) )
     $   iuplo = 2
      IF( lsame( compq, 'N' ) ) THEN
         icompq = 0
      ELSE IF( lsame( compq, 'P' ) ) THEN
         icompq = 1
      ELSE IF( lsame( compq, 'I' ) ) THEN
         icompq = 2
      ELSE
         icompq = -1
      END IF
      IF( iuplo.EQ.0 ) THEN
         info = -1
      ELSE IF( icompq.LT.0 ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( ( ldu.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldu.LT.
     $         n ) ) ) THEN
         info = -7
      ELSE IF( ( ldvt.LT.1 ) .OR. ( ( icompq.EQ.2 ) .AND. ( ldvt.LT.
     $         n ) ) ) THEN
         info = -9
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DBDSDC', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.EQ.0 )
     $   RETURN
      smlsiz = ilaenv( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
      IF( n.EQ.1 ) THEN
         IF( icompq.EQ.1 ) THEN
            q( 1 ) = sign( one, d( 1 ) )
            q( 1+smlsiz*n ) = one
         ELSE IF( icompq.EQ.2 ) THEN
            u( 1, 1 ) = sign( one, d( 1 ) )
            vt( 1, 1 ) = one
         END IF
         d( 1 ) = abs( d( 1 ) )
         RETURN
      END IF
      nm1 = n - 1
*
*     If matrix lower bidiagonal, rotate to be upper bidiagonal
*     by applying Givens rotations on the left
*
      wstart = 1
      qstart = 3
      IF( icompq.EQ.1 ) THEN
         CALL dcopy( n, d, 1, q( 1 ), 1 )
         CALL dcopy( n-1, e, 1, q( n+1 ), 1 )
      END IF
      IF( iuplo.EQ.2 ) THEN
         qstart = 5
         IF( icompq .EQ. 2 ) wstart = 2*n - 1
         DO 10 i = 1, n - 1
            CALL dlartg( d( i ), e( i ), cs, sn, r )
            d( i ) = r
            e( i ) = sn*d( i+1 )
            d( i+1 ) = cs*d( i+1 )
            IF( icompq.EQ.1 ) THEN
               q( i+2*n ) = cs
               q( i+3*n ) = sn
            ELSE IF( icompq.EQ.2 ) THEN
               work( i ) = cs
               work( nm1+i ) = -sn
            END IF
   10    CONTINUE
      END IF
*
*     If ICOMPQ = 0, use DLASDQ to compute the singular values.
*
      IF( icompq.EQ.0 ) THEN
*        Ignore WSTART, instead using WORK( 1 ), since the two vectors
*        for CS and -SN above are added only if ICOMPQ == 2,
*        and adding them exceeds documented WORK size of 4*n.
         CALL dlasdq( 'U', 0, n, 0, 0, 0, d, e, vt, ldvt, u, ldu, u,
     $                ldu, work( 1 ), info )
         GO TO 40
      END IF
*
*     If N is smaller than the minimum divide size SMLSIZ, then solve
*     the problem with another solver.
*
      IF( n.LE.smlsiz ) THEN
         IF( icompq.EQ.2 ) THEN
            CALL dlaset( 'A', n, n, zero, one, u, ldu )
            CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
            CALL dlasdq( 'U', 0, n, n, n, 0, d, e, vt, ldvt, u, ldu, u,
     $                   ldu, work( wstart ), info )
         ELSE IF( icompq.EQ.1 ) THEN
            iu = 1
            ivt = iu + n
            CALL dlaset( 'A', n, n, zero, one, q( iu+( qstart-1 )*n ),
     $                   n )
            CALL dlaset( 'A', n, n, zero, one, q( ivt+( qstart-1 )*n ),
     $                   n )
            CALL dlasdq( 'U', 0, n, n, n, 0, d, e,
     $                   q( ivt+( qstart-1 )*n ), n,
     $                   q( iu+( qstart-1 )*n ), n,
     $                   q( iu+( qstart-1 )*n ), n, work( wstart ),
     $                   info )
         END IF
         GO TO 40
      END IF
*
      IF( icompq.EQ.2 ) THEN
         CALL dlaset( 'A', n, n, zero, one, u, ldu )
         CALL dlaset( 'A', n, n, zero, one, vt, ldvt )
      END IF
*
*     Scale.
*
      orgnrm = dlanst( 'M', n, d, e )
      IF( orgnrm.EQ.zero )
     $   RETURN
      CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, ierr )
      CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, ierr )
*
      eps = (0.9d+0)*dlamch( 'Epsilon' )
*
      mlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
      smlszp = smlsiz + 1
*
      IF( icompq.EQ.1 ) THEN
         iu = 1
         ivt = 1 + smlsiz
         difl = ivt + smlszp
         difr = difl + mlvl
         z = difr + mlvl*2
         ic = z + mlvl
         is = ic + 1
         poles = is + 1
         givnum = poles + 2*mlvl
*
         k = 1
         givptr = 2
         perm = 3
         givcol = perm + mlvl
      END IF
*
      DO 20 i = 1, n
         IF( abs( d( i ) ).LT.eps ) THEN
            d( i ) = sign( eps, d( i ) )
         END IF
   20 CONTINUE
*
      start = 1
      sqre = 0
*
      DO 30 i = 1, nm1
         IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
*
*           Subproblem found. First determine its size and then
*           apply divide and conquer on it.
*
            IF( i.LT.nm1 ) THEN
*
*              A subproblem with E(I) small for I < NM1.
*
               nsize = i - start + 1
            ELSE IF( abs( e( i ) ).GE.eps ) THEN
*
*              A subproblem with E(NM1) not too small but I = NM1.
*
               nsize = n - start + 1
            ELSE
*
*              A subproblem with E(NM1) small. This implies an
*              1-by-1 subproblem at D(N). Solve this 1-by-1 problem
*              first.
*
               nsize = i - start + 1
               IF( icompq.EQ.2 ) THEN
                  u( n, n ) = sign( one, d( n ) )
                  vt( n, n ) = one
               ELSE IF( icompq.EQ.1 ) THEN
                  q( n+( qstart-1 )*n ) = sign( one, d( n ) )
                  q( n+( smlsiz+qstart-1 )*n ) = one
               END IF
               d( n ) = abs( d( n ) )
            END IF
            IF( icompq.EQ.2 ) THEN
               CALL dlasd0( nsize, sqre, d( start ), e( start ),
     $                      u( start, start ), ldu, vt( start, start ),
     $                      ldvt, smlsiz, iwork, work( wstart ), info )
            ELSE
               CALL dlasda( icompq, smlsiz, nsize, sqre, d( start ),
     $                      e( start ), q( start+( iu+qstart-2 )*n ), n,
     $                      q( start+( ivt+qstart-2 )*n ),
     $                      iq( start+k*n ), q( start+( difl+qstart-2 )*
     $                      n ), q( start+( difr+qstart-2 )*n ),
     $                      q( start+( z+qstart-2 )*n ),
     $                      q( start+( poles+qstart-2 )*n ),
     $                      iq( start+givptr*n ), iq( start+givcol*n ),
     $                      n, iq( start+perm*n ),
     $                      q( start+( givnum+qstart-2 )*n ),
     $                      q( start+( ic+qstart-2 )*n ),
     $                      q( start+( is+qstart-2 )*n ),
     $                      work( wstart ), iwork, info )
            END IF
            IF( info.NE.0 ) THEN
               RETURN
            END IF
            start = i + 1
         END IF
   30 CONTINUE
*
*     Unscale
*
      CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, ierr )
   40 CONTINUE
*
*     Use Selection Sort to minimize swaps of singular vectors
*
      DO 60 ii = 2, n
         i = ii - 1
         kk = i
         p = d( i )
         DO 50 j = ii, n
            IF( d( j ).GT.p ) THEN
               kk = j
               p = d( j )
            END IF
   50    CONTINUE
         IF( kk.NE.i ) THEN
            d( kk ) = d( i )
            d( i ) = p
            IF( icompq.EQ.1 ) THEN
               iq( i ) = kk
            ELSE IF( icompq.EQ.2 ) THEN
               CALL dswap( n, u( 1, i ), 1, u( 1, kk ), 1 )
               CALL dswap( n, vt( i, 1 ), ldvt, vt( kk, 1 ), ldvt )
            END IF
         ELSE IF( icompq.EQ.1 ) THEN
            iq( i ) = i
         END IF
   60 CONTINUE
*
*     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
*
      IF( icompq.EQ.1 ) THEN
         IF( iuplo.EQ.1 ) THEN
            iq( n ) = 1
         ELSE
            iq( n ) = 0
         END IF
      END IF
*
*     If B is lower bidiagonal, update U by those Givens rotations
*     which rotated B to be upper bidiagonal
*
      IF( ( iuplo.EQ.2 ) .AND. ( icompq.EQ.2 ) )
     $   CALL dlasr( 'L', 'V', 'B', n, n, work( 1 ), work( n ), u, ldu )
*
      RETURN
*
*     End of DBDSDC
*
      END