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*> \ingroup OTHERauxiliary 201* 202*> \par Contributors: 203* ================== 204*> 205*> Ming Gu and Huan Ren, Computer Science Division, University of 206*> California at Berkeley, USA 207*> 208* ===================================================================== 209 SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, 210 $ U, LDU, C, LDC, WORK, INFO ) 211* 212* -- LAPACK auxiliary routine -- 213* -- LAPACK is a software package provided by Univ. of Tennessee, -- 214* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 215* 216* .. Scalar Arguments .. 217 CHARACTER UPLO 218 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE 219* .. 220* .. Array Arguments .. 221 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), 222 $ VT( LDVT, * ), WORK( * ) 223* .. 224* 225* ===================================================================== 226* 227* .. Parameters .. 228 DOUBLE PRECISION ZERO 229 PARAMETER ( ZERO = 0.0D+0 ) 230* .. 231* .. Local Scalars .. 232 LOGICAL ROTATE 233 INTEGER I, ISUB, IUPLO, J, NP1, SQRE1 234 DOUBLE PRECISION CS, R, SMIN, SN 235* .. 236* .. External Subroutines .. 237 EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA 238* .. 239* .. External Functions .. 240 LOGICAL LSAME 241 EXTERNAL LSAME 242* .. 243* .. Intrinsic Functions .. 244 INTRINSIC MAX 245* .. 246* .. Executable Statements .. 247* 248* Test the input parameters. 249* 250 INFO = 0 251 IUPLO = 0 252 IF( LSAME( UPLO, 'U' ) ) 253 $ IUPLO = 1 254 IF( LSAME( UPLO, 'L' ) ) 255 $ IUPLO = 2 256 IF( IUPLO.EQ.0 ) THEN 257 INFO = -1 258 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 259 INFO = -2 260 ELSE IF( N.LT.0 ) THEN 261 INFO = -3 262 ELSE IF( NCVT.LT.0 ) THEN 263 INFO = -4 264 ELSE IF( NRU.LT.0 ) THEN 265 INFO = -5 266 ELSE IF( NCC.LT.0 ) THEN 267 INFO = -6 268 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. 269 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN 270 INFO = -10 271 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN 272 INFO = -12 273 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. 274 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN 275 INFO = -14 276 END IF 277 IF( INFO.NE.0 ) THEN 278 CALL XERBLA( 'DLASDQ', -INFO ) 279 RETURN 280 END IF 281 IF( N.EQ.0 ) 282 $ RETURN 283* 284* ROTATE is true if any singular vectors desired, false otherwise 285* 286 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) 287 NP1 = N + 1 288 SQRE1 = SQRE 289* 290* If matrix non-square upper bidiagonal, rotate to be lower 291* bidiagonal. The rotations are on the right. 292* 293 IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN 294 DO 10 I = 1, N - 1 295 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 296 D( I ) = R 297 E( I ) = SN*D( I+1 ) 298 D( I+1 ) = CS*D( I+1 ) 299 IF( ROTATE ) THEN 300 WORK( I ) = CS 301 WORK( N+I ) = SN 302 END IF 303 10 CONTINUE 304 CALL DLARTG( D( N ), E( N ), CS, SN, R ) 305 D( N ) = R 306 E( N ) = ZERO 307 IF( ROTATE ) THEN 308 WORK( N ) = CS 309 WORK( N+N ) = SN 310 END IF 311 IUPLO = 2 312 SQRE1 = 0 313* 314* Update singular vectors if desired. 315* 316 IF( NCVT.GT.0 ) 317 $ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ), 318 $ WORK( NP1 ), VT, LDVT ) 319 END IF 320* 321* If matrix lower bidiagonal, rotate to be upper bidiagonal 322* by applying Givens rotations on the left. 323* 324 IF( IUPLO.EQ.2 ) THEN 325 DO 20 I = 1, N - 1 326 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 327 D( I ) = R 328 E( I ) = SN*D( I+1 ) 329 D( I+1 ) = CS*D( I+1 ) 330 IF( ROTATE ) THEN 331 WORK( I ) = CS 332 WORK( N+I ) = SN 333 END IF 334 20 CONTINUE 335* 336* If matrix (N+1)-by-N lower bidiagonal, one additional 337* rotation is needed. 338* 339 IF( SQRE1.EQ.1 ) THEN 340 CALL DLARTG( D( N ), E( N ), CS, SN, R ) 341 D( N ) = R 342 IF( ROTATE ) THEN 343 WORK( N ) = CS 344 WORK( N+N ) = SN 345 END IF 346 END IF 347* 348* Update singular vectors if desired. 349* 350 IF( NRU.GT.0 ) THEN 351 IF( SQRE1.EQ.0 ) THEN 352 CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), 353 $ WORK( NP1 ), U, LDU ) 354 ELSE 355 CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ), 356 $ WORK( NP1 ), U, LDU ) 357 END IF 358 END IF 359 IF( NCC.GT.0 ) THEN 360 IF( SQRE1.EQ.0 ) THEN 361 CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), 362 $ WORK( NP1 ), C, LDC ) 363 ELSE 364 CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ), 365 $ WORK( NP1 ), C, LDC ) 366 END IF 367 END IF 368 END IF 369* 370* Call DBDSQR to compute the SVD of the reduced real 371* N-by-N upper bidiagonal matrix. 372* 373 CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, 374 $ LDC, WORK, INFO ) 375* 376* Sort the singular values into ascending order (insertion sort on 377* singular values, but only one transposition per singular vector) 378* 379 DO 40 I = 1, N 380* 381* Scan for smallest D(I). 382* 383 ISUB = I 384 SMIN = D( I ) 385 DO 30 J = I + 1, N 386 IF( D( J ).LT.SMIN ) THEN 387 ISUB = J 388 SMIN = D( J ) 389 END IF 390 30 CONTINUE 391 IF( ISUB.NE.I ) THEN 392* 393* Swap singular values and vectors. 394* 395 D( ISUB ) = D( I ) 396 D( I ) = SMIN 397 IF( NCVT.GT.0 ) 398 $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT ) 399 IF( NRU.GT.0 ) 400 $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 ) 401 IF( NCC.GT.0 ) 402 $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC ) 403 END IF 404 40 CONTINUE 405* 406 RETURN 407* 408* End of DLASDQ 409* 410 END 411