SRC/dbdsdc.f

   1 *> \brief \b DBDSDC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DBDSDC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
  22 *                          WORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          COMPQ, UPLO
  26 *       INTEGER            INFO, LDU, LDVT, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IQ( * ), IWORK( * )
  30 *       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
  31 *      $                   VT( LDVT, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DBDSDC computes the singular value decomposition (SVD) of a real
  41 *> N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
  42 *> using a divide and conquer method, where S is a diagonal matrix
  43 *> with non-negative diagonal elements (the singular values of B), and
  44 *> U and VT are orthogonal matrices of left and right singular vectors,
  45 *> respectively. DBDSDC can be used to compute all singular values,
  46 *> and optionally, singular vectors or singular vectors in compact form.
  47 *>
  48 *> This code makes very mild assumptions about floating point
  49 *> arithmetic. It will work on machines with a guard digit in
  50 *> add/subtract, or on those binary machines without guard digits
  51 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
  52 *> It could conceivably fail on hexadecimal or decimal machines
  53 *> without guard digits, but we know of none.  See DLASD3 for details.
  54 *>
  55 *> The code currently calls DLASDQ if singular values only are desired.
  56 *> However, it can be slightly modified to compute singular values
  57 *> using the divide and conquer method.
  58 *> \endverbatim
  59 *
  60 *  Arguments:
  61 *  ==========
  62 *
  63 *> \param[in] UPLO
  64 *> \verbatim
  65 *>          UPLO is CHARACTER*1
  66 *>          = 'U':  B is upper bidiagonal.
  67 *>          = 'L':  B is lower bidiagonal.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] COMPQ
  71 *> \verbatim
  72 *>          COMPQ is CHARACTER*1
  73 *>          Specifies whether singular vectors are to be computed
  74 *>          as follows:
  75 *>          = 'N':  Compute singular values only;
  76 *>          = 'P':  Compute singular values and compute singular
  77 *>                  vectors in compact form;
  78 *>          = 'I':  Compute singular values and singular vectors.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] N
  82 *> \verbatim
  83 *>          N is INTEGER
  84 *>          The order of the matrix B.  N >= 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in,out] D
  88 *> \verbatim
  89 *>          D is DOUBLE PRECISION array, dimension (N)
  90 *>          On entry, the n diagonal elements of the bidiagonal matrix B.
  91 *>          On exit, if INFO=0, the singular values of B.
  92 *> \endverbatim
  93 *>
  94 *> \param[in,out] E
  95 *> \verbatim
  96 *>          E is DOUBLE PRECISION array, dimension (N-1)
  97 *>          On entry, the elements of E contain the offdiagonal
  98 *>          elements of the bidiagonal matrix whose SVD is desired.
  99 *>          On exit, E has been destroyed.
 100 *> \endverbatim
 101 *>
 102 *> \param[out] U
 103 *> \verbatim
 104 *>          U is DOUBLE PRECISION array, dimension (LDU,N)
 105 *>          If  COMPQ = 'I', then:
 106 *>             On exit, if INFO = 0, U contains the left singular vectors
 107 *>             of the bidiagonal matrix.
 108 *>          For other values of COMPQ, U is not referenced.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDU
 112 *> \verbatim
 113 *>          LDU is INTEGER
 114 *>          The leading dimension of the array U.  LDU >= 1.
 115 *>          If singular vectors are desired, then LDU >= max( 1, N ).
 116 *> \endverbatim
 117 *>
 118 *> \param[out] VT
 119 *> \verbatim
 120 *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
 121 *>          If  COMPQ = 'I', then:
 122 *>             On exit, if INFO = 0, VT**T contains the right singular
 123 *>             vectors of the bidiagonal matrix.
 124 *>          For other values of COMPQ, VT is not referenced.
 125 *> \endverbatim
 126 *>
 127 *> \param[in] LDVT
 128 *> \verbatim
 129 *>          LDVT is INTEGER
 130 *>          The leading dimension of the array VT.  LDVT >= 1.
 131 *>          If singular vectors are desired, then LDVT >= max( 1, N ).
 132 *> \endverbatim
 133 *>
 134 *> \param[out] Q
 135 *> \verbatim
 136 *>          Q is DOUBLE PRECISION array, dimension (LDQ)
 137 *>          If  COMPQ = 'P', then:
 138 *>             On exit, if INFO = 0, Q and IQ contain the left
 139 *>             and right singular vectors in a compact form,
 140 *>             requiring O(N log N) space instead of 2*N**2.
 141 *>             In particular, Q contains all the DOUBLE PRECISION data in
 142 *>             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
 143 *>             words of memory, where SMLSIZ is returned by ILAENV and
 144 *>             is equal to the maximum size of the subproblems at the
 145 *>             bottom of the computation tree (usually about 25).
 146 *>          For other values of COMPQ, Q is not referenced.
 147 *> \endverbatim
 148 *>
 149 *> \param[out] IQ
 150 *> \verbatim
 151 *>          IQ is INTEGER array, dimension (LDIQ)
 152 *>          If  COMPQ = 'P', then:
 153 *>             On exit, if INFO = 0, Q and IQ contain the left
 154 *>             and right singular vectors in a compact form,
 155 *>             requiring O(N log N) space instead of 2*N**2.
 156 *>             In particular, IQ contains all INTEGER data in
 157 *>             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
 158 *>             words of memory, where SMLSIZ is returned by ILAENV and
 159 *>             is equal to the maximum size of the subproblems at the
 160 *>             bottom of the computation tree (usually about 25).
 161 *>          For other values of COMPQ, IQ is not referenced.
 162 *> \endverbatim
 163 *>
 164 *> \param[out] WORK
 165 *> \verbatim
 166 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 167 *>          If COMPQ = 'N' then LWORK >= (4 * N).
 168 *>          If COMPQ = 'P' then LWORK >= (6 * N).
 169 *>          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
 170 *> \endverbatim
 171 *>
 172 *> \param[out] IWORK
 173 *> \verbatim
 174 *>          IWORK is INTEGER array, dimension (8*N)
 175 *> \endverbatim
 176 *>
 177 *> \param[out] INFO
 178 *> \verbatim
 179 *>          INFO is INTEGER
 180 *>          = 0:  successful exit.
 181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 182 *>          > 0:  The algorithm failed to compute a singular value.
 183 *>                The update process of divide and conquer failed.
 184 *> \endverbatim
 185 *
 186 *  Authors:
 187 *  ========
 188 *
 189 *> \author Univ. of Tennessee
 190 *> \author Univ. of California Berkeley
 191 *> \author Univ. of Colorado Denver
 192 *> \author NAG Ltd.
 193 *
 194 *> \date June 2016
 195 *
 196 *> \ingroup auxOTHERcomputational
 197 *
 198 *> \par Contributors:
 199 *  ==================
 200 *>
 201 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 202 *>     California at Berkeley, USA
 203 *>
 204 *  =====================================================================
 205       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
 206      $                   WORK, IWORK, INFO )
 207 *
 208 *  -- LAPACK computational routine (version 3.6.1) --
 209 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 210 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 211 *     June 2016
 212 *
 213 *     .. Scalar Arguments ..
 214       CHARACTER          COMPQ, UPLO
 215       INTEGER            INFO, LDU, LDVT, N
 216 *     ..
 217 *     .. Array Arguments ..
 218       INTEGER            IQ( * ), IWORK( * )
 219       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
 220      $                   VT( LDVT, * ), WORK( * )
 221 *     ..
 222 *
 223 *  =====================================================================
 224 *  Changed dimension statement in comment describing E from (N) to
 225 *  (N-1).  Sven, 17 Feb 05.
 226 *  =====================================================================
 227 *
 228 *     .. Parameters ..
 229       DOUBLE PRECISION   ZERO, ONE, TWO
 230       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
 231 *     ..
 232 *     .. Local Scalars ..
 233       INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
 234      $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
 235      $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
 236      $                   SMLSZP, SQRE, START, WSTART, Z
 237       DOUBLE PRECISION   CS, EPS, ORGNRM, P, R, SN
 238 *     ..
 239 *     .. External Functions ..
 240       LOGICAL            LSAME
 241       INTEGER            ILAENV
 242       DOUBLE PRECISION   DLAMCH, DLANST
 243       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
 244 *     ..
 245 *     .. External Subroutines ..
 246       EXTERNAL           DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
 247      $                   DLASET, DLASR, DSWAP, XERBLA
 248 *     ..
 249 *     .. Intrinsic Functions ..
 250       INTRINSIC          ABS, DBLE, INT, LOG, SIGN
 251 *     ..
 252 *     .. Executable Statements ..
 253 *
 254 *     Test the input parameters.
 255 *
 256       INFO = 0
 257 *
 258       IUPLO = 0
 259       IF( LSAME( UPLO, 'U' ) )
 260      $   IUPLO = 1
 261       IF( LSAME( UPLO, 'L' ) )
 262      $   IUPLO = 2
 263       IF( LSAME( COMPQ, 'N' ) ) THEN
 264          ICOMPQ = 0
 265       ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
 266          ICOMPQ = 1
 267       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
 268          ICOMPQ = 2
 269       ELSE
 270          ICOMPQ = -1
 271       END IF
 272       IF( IUPLO.EQ.0 ) THEN
 273          INFO = -1
 274       ELSE IF( ICOMPQ.LT.0 ) THEN
 275          INFO = -2
 276       ELSE IF( N.LT.0 ) THEN
 277          INFO = -3
 278       ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
 279      $         N ) ) ) THEN
 280          INFO = -7
 281       ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
 282      $         N ) ) ) THEN
 283          INFO = -9
 284       END IF
 285       IF( INFO.NE.0 ) THEN
 286          CALL XERBLA( 'DBDSDC', -INFO )
 287          RETURN
 288       END IF
 289 *
 290 *     Quick return if possible
 291 *
 292       IF( N.EQ.0 )
 293      $   RETURN
 294       SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
 295       IF( N.EQ.1 ) THEN
 296          IF( ICOMPQ.EQ.1 ) THEN
 297             Q( 1 ) = SIGN( ONE, D( 1 ) )
 298             Q( 1+SMLSIZ*N ) = ONE
 299          ELSE IF( ICOMPQ.EQ.2 ) THEN
 300             U( 1, 1 ) = SIGN( ONE, D( 1 ) )
 301             VT( 1, 1 ) = ONE
 302          END IF
 303          D( 1 ) = ABS( D( 1 ) )
 304          RETURN
 305       END IF
 306       NM1 = N - 1
 307 *
 308 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
 309 *     by applying Givens rotations on the left
 310 *
 311       WSTART = 1
 312       QSTART = 3
 313       IF( ICOMPQ.EQ.1 ) THEN
 314          CALL DCOPY( N,   D, 1, Q( 1 ),   1 )
 315          CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
 316       END IF
 317       IF( IUPLO.EQ.2 ) THEN
 318          QSTART = 5
 319          WSTART = 2*N - 1
 320          DO 10 I = 1, N - 1
 321             CALL DLARTG( D( I ), E( I ), CS, SN, R )
 322             D( I ) = R
 323             E( I ) = SN*D( I+1 )
 324             D( I+1 ) = CS*D( I+1 )
 325             IF( ICOMPQ.EQ.1 ) THEN
 326                Q( I+2*N ) = CS
 327                Q( I+3*N ) = SN
 328             ELSE IF( ICOMPQ.EQ.2 ) THEN
 329                WORK( I ) = CS
 330                WORK( NM1+I ) = -SN
 331             END IF
 332    10    CONTINUE
 333       END IF
 334 *
 335 *     If ICOMPQ = 0, use DLASDQ to compute the singular values.
 336 *
 337       IF( ICOMPQ.EQ.0 ) THEN
 338 *        Ignore WSTART, instead using WORK( 1 ), since the two vectors
 339 *        for CS and -SN above are added only if ICOMPQ == 2,
 340 *        and adding them exceeds documented WORK size of 4*n.
 341          CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
 342      $                LDU, WORK( 1 ), INFO )
 343          GO TO 40
 344       END IF
 345 *
 346 *     If N is smaller than the minimum divide size SMLSIZ, then solve
 347 *     the problem with another solver.
 348 *
 349       IF( N.LE.SMLSIZ ) THEN
 350          IF( ICOMPQ.EQ.2 ) THEN
 351             CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
 352             CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
 353             CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
 354      $                   LDU, WORK( WSTART ), INFO )
 355          ELSE IF( ICOMPQ.EQ.1 ) THEN
 356             IU = 1
 357             IVT = IU + N
 358             CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
 359      $                   N )
 360             CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
 361      $                   N )
 362             CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
 363      $                   Q( IVT+( QSTART-1 )*N ), N,
 364      $                   Q( IU+( QSTART-1 )*N ), N,
 365      $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
 366      $                   INFO )
 367          END IF
 368          GO TO 40
 369       END IF
 370 *
 371       IF( ICOMPQ.EQ.2 ) THEN
 372          CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
 373          CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
 374       END IF
 375 *
 376 *     Scale.
 377 *
 378       ORGNRM = DLANST( 'M', N, D, E )
 379       IF( ORGNRM.EQ.ZERO )
 380      $   RETURN
 381       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
 382       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
 383 *
 384       EPS = (0.9D+0)*DLAMCH( 'Epsilon' )
 385 *
 386       MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
 387       SMLSZP = SMLSIZ + 1
 388 *
 389       IF( ICOMPQ.EQ.1 ) THEN
 390          IU = 1
 391          IVT = 1 + SMLSIZ
 392          DIFL = IVT + SMLSZP
 393          DIFR = DIFL + MLVL
 394          Z = DIFR + MLVL*2
 395          IC = Z + MLVL
 396          IS = IC + 1
 397          POLES = IS + 1
 398          GIVNUM = POLES + 2*MLVL
 399 *
 400          K = 1
 401          GIVPTR = 2
 402          PERM = 3
 403          GIVCOL = PERM + MLVL
 404       END IF
 405 *
 406       DO 20 I = 1, N
 407          IF( ABS( D( I ) ).LT.EPS ) THEN
 408             D( I ) = SIGN( EPS, D( I ) )
 409          END IF
 410    20 CONTINUE
 411 *
 412       START = 1
 413       SQRE = 0
 414 *
 415       DO 30 I = 1, NM1
 416          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
 417 *
 418 *           Subproblem found. First determine its size and then
 419 *           apply divide and conquer on it.
 420 *
 421             IF( I.LT.NM1 ) THEN
 422 *
 423 *              A subproblem with E(I) small for I < NM1.
 424 *
 425                NSIZE = I - START + 1
 426             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
 427 *
 428 *              A subproblem with E(NM1) not too small but I = NM1.
 429 *
 430                NSIZE = N - START + 1
 431             ELSE
 432 *
 433 *              A subproblem with E(NM1) small. This implies an
 434 *              1-by-1 subproblem at D(N). Solve this 1-by-1 problem
 435 *              first.
 436 *
 437                NSIZE = I - START + 1
 438                IF( ICOMPQ.EQ.2 ) THEN
 439                   U( N, N ) = SIGN( ONE, D( N ) )
 440                   VT( N, N ) = ONE
 441                ELSE IF( ICOMPQ.EQ.1 ) THEN
 442                   Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
 443                   Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
 444                END IF
 445                D( N ) = ABS( D( N ) )
 446             END IF
 447             IF( ICOMPQ.EQ.2 ) THEN
 448                CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
 449      $                      U( START, START ), LDU, VT( START, START ),
 450      $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
 451             ELSE
 452                CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
 453      $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
 454      $                      Q( START+( IVT+QSTART-2 )*N ),
 455      $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
 456      $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
 457      $                      Q( START+( Z+QSTART-2 )*N ),
 458      $                      Q( START+( POLES+QSTART-2 )*N ),
 459      $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
 460      $                      N, IQ( START+PERM*N ),
 461      $                      Q( START+( GIVNUM+QSTART-2 )*N ),
 462      $                      Q( START+( IC+QSTART-2 )*N ),
 463      $                      Q( START+( IS+QSTART-2 )*N ),
 464      $                      WORK( WSTART ), IWORK, INFO )
 465             END IF
 466             IF( INFO.NE.0 ) THEN
 467                RETURN
 468             END IF
 469             START = I + 1
 470          END IF
 471    30 CONTINUE
 472 *
 473 *     Unscale
 474 *
 475       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
 476    40 CONTINUE
 477 *
 478 *     Use Selection Sort to minimize swaps of singular vectors
 479 *
 480       DO 60 II = 2, N
 481          I = II - 1
 482          KK = I
 483          P = D( I )
 484          DO 50 J = II, N
 485             IF( D( J ).GT.P ) THEN
 486                KK = J
 487                P = D( J )
 488             END IF
 489    50    CONTINUE
 490          IF( KK.NE.I ) THEN
 491             D( KK ) = D( I )
 492             D( I ) = P
 493             IF( ICOMPQ.EQ.1 ) THEN
 494                IQ( I ) = KK
 495             ELSE IF( ICOMPQ.EQ.2 ) THEN
 496                CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
 497                CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
 498             END IF
 499          ELSE IF( ICOMPQ.EQ.1 ) THEN
 500             IQ( I ) = I
 501          END IF
 502    60 CONTINUE
 503 *
 504 *     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
 505 *
 506       IF( ICOMPQ.EQ.1 ) THEN
 507          IF( IUPLO.EQ.1 ) THEN
 508             IQ( N ) = 1
 509          ELSE
 510             IQ( N ) = 0
 511          END IF
 512       END IF
 513 *
 514 *     If B is lower bidiagonal, update U by those Givens rotations
 515 *     which rotated B to be upper bidiagonal
 516 *
 517       IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
 518      $   CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
 519 *
 520       RETURN
 521 *
 522 *     End of DBDSDC
 523 *
 524       END