SRC/dlasdq.f

   1 *> \brief \b DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLASDQ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasdq.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasdq.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasdq.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
  22 *                          U, LDU, C, LDC, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
  30 *      $                   VT( LDVT, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DLASDQ computes the singular value decomposition (SVD) of a real
  40 *> (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
  41 *> E, accumulating the transformations if desired. Letting B denote
  42 *> the input bidiagonal matrix, the algorithm computes orthogonal
  43 *> matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose
  44 *> of P). The singular values S are overwritten on D.
  45 *>
  46 *> The input matrix U  is changed to U  * Q  if desired.
  47 *> The input matrix VT is changed to P**T * VT if desired.
  48 *> The input matrix C  is changed to Q**T * C  if desired.
  49 *>
  50 *> See "Computing  Small Singular Values of Bidiagonal Matrices With
  51 *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
  52 *> LAPACK Working Note #3, for a detailed description of the algorithm.
  53 *> \endverbatim
  54 *
  55 *  Arguments:
  56 *  ==========
  57 *
  58 *> \param[in] UPLO
  59 *> \verbatim
  60 *>          UPLO is CHARACTER*1
  61 *>        On entry, UPLO specifies whether the input bidiagonal matrix
  62 *>        is upper or lower bidiagonal, and whether it is square are
  63 *>        not.
  64 *>           UPLO = 'U' or 'u'   B is upper bidiagonal.
  65 *>           UPLO = 'L' or 'l'   B is lower bidiagonal.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] SQRE
  69 *> \verbatim
  70 *>          SQRE is INTEGER
  71 *>        = 0: then the input matrix is N-by-N.
  72 *>        = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
  73 *>             (N+1)-by-N if UPLU = 'L'.
  74 *>
  75 *>        The bidiagonal matrix has
  76 *>        N = NL + NR + 1 rows and
  77 *>        M = N + SQRE >= N columns.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] N
  81 *> \verbatim
  82 *>          N is INTEGER
  83 *>        On entry, N specifies the number of rows and columns
  84 *>        in the matrix. N must be at least 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] NCVT
  88 *> \verbatim
  89 *>          NCVT is INTEGER
  90 *>        On entry, NCVT specifies the number of columns of
  91 *>        the matrix VT. NCVT must be at least 0.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] NRU
  95 *> \verbatim
  96 *>          NRU is INTEGER
  97 *>        On entry, NRU specifies the number of rows of
  98 *>        the matrix U. NRU must be at least 0.
  99 *> \endverbatim
 100 *>
 101 *> \param[in] NCC
 102 *> \verbatim
 103 *>          NCC is INTEGER
 104 *>        On entry, NCC specifies the number of columns of
 105 *>        the matrix C. NCC must be at least 0.
 106 *> \endverbatim
 107 *>
 108 *> \param[in,out] D
 109 *> \verbatim
 110 *>          D is DOUBLE PRECISION array, dimension (N)
 111 *>        On entry, D contains the diagonal entries of the
 112 *>        bidiagonal matrix whose SVD is desired. On normal exit,
 113 *>        D contains the singular values in ascending order.
 114 *> \endverbatim
 115 *>
 116 *> \param[in,out] E
 117 *> \verbatim
 118 *>          E is DOUBLE PRECISION array.
 119 *>        dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
 120 *>        On entry, the entries of E contain the offdiagonal entries
 121 *>        of the bidiagonal matrix whose SVD is desired. On normal
 122 *>        exit, E will contain 0. If the algorithm does not converge,
 123 *>        D and E will contain the diagonal and superdiagonal entries
 124 *>        of a bidiagonal matrix orthogonally equivalent to the one
 125 *>        given as input.
 126 *> \endverbatim
 127 *>
 128 *> \param[in,out] VT
 129 *> \verbatim
 130 *>          VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
 131 *>        On entry, contains a matrix which on exit has been
 132 *>        premultiplied by P**T, dimension N-by-NCVT if SQRE = 0
 133 *>        and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
 134 *> \endverbatim
 135 *>
 136 *> \param[in] LDVT
 137 *> \verbatim
 138 *>          LDVT is INTEGER
 139 *>        On entry, LDVT specifies the leading dimension of VT as
 140 *>        declared in the calling (sub) program. LDVT must be at
 141 *>        least 1. If NCVT is nonzero LDVT must also be at least N.
 142 *> \endverbatim
 143 *>
 144 *> \param[in,out] U
 145 *> \verbatim
 146 *>          U is DOUBLE PRECISION array, dimension (LDU, N)
 147 *>        On entry, contains a  matrix which on exit has been
 148 *>        postmultiplied by Q, dimension NRU-by-N if SQRE = 0
 149 *>        and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
 150 *> \endverbatim
 151 *>
 152 *> \param[in] LDU
 153 *> \verbatim
 154 *>          LDU is INTEGER
 155 *>        On entry, LDU  specifies the leading dimension of U as
 156 *>        declared in the calling (sub) program. LDU must be at
 157 *>        least max( 1, NRU ) .
 158 *> \endverbatim
 159 *>
 160 *> \param[in,out] C
 161 *> \verbatim
 162 *>          C is DOUBLE PRECISION array, dimension (LDC, NCC)
 163 *>        On entry, contains an N-by-NCC matrix which on exit
 164 *>        has been premultiplied by Q**T  dimension N-by-NCC if SQRE = 0
 165 *>        and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
 166 *> \endverbatim
 167 *>
 168 *> \param[in] LDC
 169 *> \verbatim
 170 *>          LDC is INTEGER
 171 *>        On entry, LDC  specifies the leading dimension of C as
 172 *>        declared in the calling (sub) program. LDC must be at
 173 *>        least 1. If NCC is nonzero, LDC must also be at least N.
 174 *> \endverbatim
 175 *>
 176 *> \param[out] WORK
 177 *> \verbatim
 178 *>          WORK is DOUBLE PRECISION array, dimension (4*N)
 179 *>        Workspace. Only referenced if one of NCVT, NRU, or NCC is
 180 *>        nonzero, and if N is at least 2.
 181 *> \endverbatim
 182 *>
 183 *> \param[out] INFO
 184 *> \verbatim
 185 *>          INFO is INTEGER
 186 *>        On exit, a value of 0 indicates a successful exit.
 187 *>        If INFO < 0, argument number -INFO is illegal.
 188 *>        If INFO > 0, the algorithm did not converge, and INFO
 189 *>        specifies how many superdiagonals did not converge.
 190 *> \endverbatim
 191 *
 192 *  Authors:
 193 *  ========
 194 *
 195 *> \author Univ. of Tennessee
 196 *> \author Univ. of California Berkeley
 197 *> \author Univ. of Colorado Denver
 198 *> \author NAG Ltd.
 199 *
 200 *> \date June 2016
 201 *
 202 *> \ingroup auxOTHERauxiliary
 203 *
 204 *> \par Contributors:
 205 *  ==================
 206 *>
 207 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 208 *>     California at Berkeley, USA
 209 *>
 210 *  =====================================================================
 211       SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
 212      $                   U, LDU, C, LDC, WORK, INFO )
 213 *
 214 *  -- LAPACK auxiliary routine (version 3.6.1) --
 215 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 216 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 217 *     June 2016
 218 *
 219 *     .. Scalar Arguments ..
 220       CHARACTER          UPLO
 221       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
 222 *     ..
 223 *     .. Array Arguments ..
 224       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
 225      $                   VT( LDVT, * ), WORK( * )
 226 *     ..
 227 *
 228 *  =====================================================================
 229 *
 230 *     .. Parameters ..
 231       DOUBLE PRECISION   ZERO
 232       PARAMETER          ( ZERO = 0.0D+0 )
 233 *     ..
 234 *     .. Local Scalars ..
 235       LOGICAL            ROTATE
 236       INTEGER            I, ISUB, IUPLO, J, NP1, SQRE1
 237       DOUBLE PRECISION   CS, R, SMIN, SN
 238 *     ..
 239 *     .. External Subroutines ..
 240       EXTERNAL           DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
 241 *     ..
 242 *     .. External Functions ..
 243       LOGICAL            LSAME
 244       EXTERNAL           LSAME
 245 *     ..
 246 *     .. Intrinsic Functions ..
 247       INTRINSIC          MAX
 248 *     ..
 249 *     .. Executable Statements ..
 250 *
 251 *     Test the input parameters.
 252 *
 253       INFO = 0
 254       IUPLO = 0
 255       IF( LSAME( UPLO, 'U' ) )
 256      $   IUPLO = 1
 257       IF( LSAME( UPLO, 'L' ) )
 258      $   IUPLO = 2
 259       IF( IUPLO.EQ.0 ) THEN
 260          INFO = -1
 261       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
 262          INFO = -2
 263       ELSE IF( N.LT.0 ) THEN
 264          INFO = -3
 265       ELSE IF( NCVT.LT.0 ) THEN
 266          INFO = -4
 267       ELSE IF( NRU.LT.0 ) THEN
 268          INFO = -5
 269       ELSE IF( NCC.LT.0 ) THEN
 270          INFO = -6
 271       ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
 272      $         ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
 273          INFO = -10
 274       ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
 275          INFO = -12
 276       ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
 277      $         ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
 278          INFO = -14
 279       END IF
 280       IF( INFO.NE.0 ) THEN
 281          CALL XERBLA( 'DLASDQ', -INFO )
 282          RETURN
 283       END IF
 284       IF( N.EQ.0 )
 285      $   RETURN
 286 *
 287 *     ROTATE is true if any singular vectors desired, false otherwise
 288 *
 289       ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
 290       NP1 = N + 1
 291       SQRE1 = SQRE
 292 *
 293 *     If matrix non-square upper bidiagonal, rotate to be lower
 294 *     bidiagonal.  The rotations are on the right.
 295 *
 296       IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
 297          DO 10 I = 1, N - 1
 298             CALL DLARTG( D( I ), E( I ), CS, SN, R )
 299             D( I ) = R
 300             E( I ) = SN*D( I+1 )
 301             D( I+1 ) = CS*D( I+1 )
 302             IF( ROTATE ) THEN
 303                WORK( I ) = CS
 304                WORK( N+I ) = SN
 305             END IF
 306    10    CONTINUE
 307          CALL DLARTG( D( N ), E( N ), CS, SN, R )
 308          D( N ) = R
 309          E( N ) = ZERO
 310          IF( ROTATE ) THEN
 311             WORK( N ) = CS
 312             WORK( N+N ) = SN
 313          END IF
 314          IUPLO = 2
 315          SQRE1 = 0
 316 *
 317 *        Update singular vectors if desired.
 318 *
 319          IF( NCVT.GT.0 )
 320      $      CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
 321      $                  WORK( NP1 ), VT, LDVT )
 322       END IF
 323 *
 324 *     If matrix lower bidiagonal, rotate to be upper bidiagonal
 325 *     by applying Givens rotations on the left.
 326 *
 327       IF( IUPLO.EQ.2 ) THEN
 328          DO 20 I = 1, N - 1
 329             CALL DLARTG( D( I ), E( I ), CS, SN, R )
 330             D( I ) = R
 331             E( I ) = SN*D( I+1 )
 332             D( I+1 ) = CS*D( I+1 )
 333             IF( ROTATE ) THEN
 334                WORK( I ) = CS
 335                WORK( N+I ) = SN
 336             END IF
 337    20    CONTINUE
 338 *
 339 *        If matrix (N+1)-by-N lower bidiagonal, one additional
 340 *        rotation is needed.
 341 *
 342          IF( SQRE1.EQ.1 ) THEN
 343             CALL DLARTG( D( N ), E( N ), CS, SN, R )
 344             D( N ) = R
 345             IF( ROTATE ) THEN
 346                WORK( N ) = CS
 347                WORK( N+N ) = SN
 348             END IF
 349          END IF
 350 *
 351 *        Update singular vectors if desired.
 352 *
 353          IF( NRU.GT.0 ) THEN
 354             IF( SQRE1.EQ.0 ) THEN
 355                CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
 356      $                     WORK( NP1 ), U, LDU )
 357             ELSE
 358                CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
 359      $                     WORK( NP1 ), U, LDU )
 360             END IF
 361          END IF
 362          IF( NCC.GT.0 ) THEN
 363             IF( SQRE1.EQ.0 ) THEN
 364                CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
 365      $                     WORK( NP1 ), C, LDC )
 366             ELSE
 367                CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
 368      $                     WORK( NP1 ), C, LDC )
 369             END IF
 370          END IF
 371       END IF
 372 *
 373 *     Call DBDSQR to compute the SVD of the reduced real
 374 *     N-by-N upper bidiagonal matrix.
 375 *
 376       CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
 377      $             LDC, WORK, INFO )
 378 *
 379 *     Sort the singular values into ascending order (insertion sort on
 380 *     singular values, but only one transposition per singular vector)
 381 *
 382       DO 40 I = 1, N
 383 *
 384 *        Scan for smallest D(I).
 385 *
 386          ISUB = I
 387          SMIN = D( I )
 388          DO 30 J = I + 1, N
 389             IF( D( J ).LT.SMIN ) THEN
 390                ISUB = J
 391                SMIN = D( J )
 392             END IF
 393    30    CONTINUE
 394          IF( ISUB.NE.I ) THEN
 395 *
 396 *           Swap singular values and vectors.
 397 *
 398             D( ISUB ) = D( I )
 399             D( I ) = SMIN
 400             IF( NCVT.GT.0 )
 401      $         CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
 402             IF( NRU.GT.0 )
 403      $         CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
 404             IF( NCC.GT.0 )
 405      $         CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
 406          END IF
 407    40 CONTINUE
 408 *
 409       RETURN
 410 *
 411 *     End of DLASDQ
 412 *
 413       END