1 SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, 2 $ U, LDU, C, LDC, WORK, INFO ) 3* 4* -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* October 31, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER UPLO 11 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), 15 $ VT( LDVT, * ), WORK( * ) 16* .. 17* .. Common block to return operation count .. 18 COMMON / LATIME / OPS, ITCNT 19* .. 20* .. Scalars in Common .. 21 DOUBLE PRECISION ITCNT, OPS 22* .. 23* 24* Purpose 25* ======= 26* 27* DLASDQ computes the singular value decomposition (SVD) of a real 28* (upper or lower) bidiagonal matrix with diagonal D and offdiagonal 29* E, accumulating the transformations if desired. Letting B denote 30* the input bidiagonal matrix, the algorithm computes orthogonal 31* matrices Q and P such that B = Q * S * P' (P' denotes the transpose 32* of P). The singular values S are overwritten on D. 33* 34* The input matrix U is changed to U * Q if desired. 35* The input matrix VT is changed to P' * VT if desired. 36* The input matrix C is changed to Q' * C if desired. 37* 38* See "Computing Small Singular Values of Bidiagonal Matrices With 39* Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 40* LAPACK Working Note #3, for a detailed description of the algorithm. 41* 42* Arguments 43* ========= 44* 45* UPLO (input) CHARACTER*1 46* On entry, UPLO specifies whether the input bidiagonal matrix 47* is upper or lower bidiagonal, and wether it is square are 48* not. 49* UPLO = 'U' or 'u' B is upper bidiagonal. 50* UPLO = 'L' or 'l' B is lower bidiagonal. 51* 52* SQRE (input) INTEGER 53* = 0: then the input matrix is N-by-N. 54* = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and 55* (N+1)-by-N if UPLU = 'L'. 56* 57* The bidiagonal matrix has 58* N = NL + NR + 1 rows and 59* M = N + SQRE >= N columns. 60* 61* N (input) INTEGER 62* On entry, N specifies the number of rows and columns 63* in the matrix. N must be at least 0. 64* 65* NCVT (input) INTEGER 66* On entry, NCVT specifies the number of columns of 67* the matrix VT. NCVT must be at least 0. 68* 69* NRU (input) INTEGER 70* On entry, NRU specifies the number of rows of 71* the matrix U. NRU must be at least 0. 72* 73* NCC (input) INTEGER 74* On entry, NCC specifies the number of columns of 75* the matrix C. NCC must be at least 0. 76* 77* D (input/output) DOUBLE PRECISION array, dimension (N) 78* On entry, D contains the diagonal entries of the 79* bidiagonal matrix whose SVD is desired. On normal exit, 80* D contains the singular values in ascending order. 81* 82* E (input/output) DOUBLE PRECISION array. 83* dimension is (N-1) if SQRE = 0 and N if SQRE = 1. 84* On entry, the entries of E contain the offdiagonal entries 85* of the bidiagonal matrix whose SVD is desired. On normal 86* exit, E will contain 0. If the algorithm does not converge, 87* D and E will contain the diagonal and superdiagonal entries 88* of a bidiagonal matrix orthogonally equivalent to the one 89* given as input. 90* 91* VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) 92* On entry, contains a matrix which on exit has been 93* premultiplied by P', dimension N-by-NCVT if SQRE = 0 94* and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). 95* 96* LDVT (input) INTEGER 97* On entry, LDVT specifies the leading dimension of VT as 98* declared in the calling (sub) program. LDVT must be at 99* least 1. If NCVT is nonzero LDVT must also be at least N. 100* 101* U (input/output) DOUBLE PRECISION array, dimension (LDU, N) 102* On entry, contains a matrix which on exit has been 103* postmultiplied by Q, dimension NRU-by-N if SQRE = 0 104* and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). 105* 106* LDU (input) INTEGER 107* On entry, LDU specifies the leading dimension of U as 108* declared in the calling (sub) program. LDU must be at 109* least max( 1, NRU ) . 110* 111* C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) 112* On entry, contains an N-by-NCC matrix which on exit 113* has been premultiplied by Q' dimension N-by-NCC if SQRE = 0 114* and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). 115* 116* LDC (input) INTEGER 117* On entry, LDC specifies the leading dimension of C as 118* declared in the calling (sub) program. LDC must be at 119* least 1. If NCC is nonzero, LDC must also be at least N. 120* 121* WORK (workspace) DOUBLE PRECISION array, dimension (MAX( 1, 4*N )) 122* Workspace. Only referenced if one of NCVT, NRU, or NCC is 123* nonzero, and if N is at least 2. 124* 125* INFO (output) INTEGER 126* On exit, a value of 0 indicates a successful exit. 127* If INFO < 0, argument number -INFO is illegal. 128* If INFO > 0, the algorithm did not converge, and INFO 129* specifies how many superdiagonals did not converge. 130* 131* Further Details 132* =============== 133* 134* Based on contributions by 135* Ming Gu and Huan Ren, Computer Science Division, University of 136* California at Berkeley, USA 137* 138* ===================================================================== 139* 140* .. Parameters .. 141 DOUBLE PRECISION ZERO 142 PARAMETER ( ZERO = 0.0D0 ) 143* .. 144* .. Local Scalars .. 145 LOGICAL ROTATE 146 INTEGER I, ISUB, IUPLO, J, NP1, SQRE1 147 DOUBLE PRECISION CS, R, SMIN, SN 148* .. 149* .. External Subroutines .. 150 EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA 151* .. 152* .. External Functions .. 153 LOGICAL LSAME 154 EXTERNAL LSAME 155* .. 156* .. Intrinsic Functions .. 157 INTRINSIC DBLE, MAX 158* .. 159* .. Executable Statements .. 160* 161* Test the input parameters. 162* 163 INFO = 0 164 IUPLO = 0 165 IF( LSAME( UPLO, 'U' ) ) 166 $ IUPLO = 1 167 IF( LSAME( UPLO, 'L' ) ) 168 $ IUPLO = 2 169 IF( IUPLO.EQ.0 ) THEN 170 INFO = -1 171 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 172 INFO = -2 173 ELSE IF( N.LT.0 ) THEN 174 INFO = -3 175 ELSE IF( NCVT.LT.0 ) THEN 176 INFO = -4 177 ELSE IF( NRU.LT.0 ) THEN 178 INFO = -5 179 ELSE IF( NCC.LT.0 ) THEN 180 INFO = -6 181 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. 182 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN 183 INFO = -10 184 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN 185 INFO = -12 186 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. 187 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN 188 INFO = -14 189 END IF 190 IF( INFO.NE.0 ) THEN 191 CALL XERBLA( 'DBDSQR', -INFO ) 192 RETURN 193 END IF 194 IF( N.EQ.0 ) 195 $ RETURN 196* 197* ROTATE is true if any singular vectors desired, false otherwise 198* 199 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) 200 NP1 = N + 1 201 SQRE1 = SQRE 202* 203* If matrix non-square upper bidiagonal, rotate to be lower 204* bidiagonal. The rotations are on the right. 205* 206 IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN 207 OPS = OPS + DBLE( 8*( N-1 ) ) 208 DO 10 I = 1, N - 1 209 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 210 D( I ) = R 211 E( I ) = SN*D( I+1 ) 212 D( I+1 ) = CS*D( I+1 ) 213 IF( ROTATE ) THEN 214 WORK( I ) = CS 215 WORK( N+I ) = SN 216 END IF 217 10 CONTINUE 218 OPS = OPS + DBLE( 6 ) 219 CALL DLARTG( D( N ), E( N ), CS, SN, R ) 220 D( N ) = R 221 E( N ) = ZERO 222 IF( ROTATE ) THEN 223 WORK( N ) = CS 224 WORK( N+N ) = SN 225 END IF 226 IUPLO = 2 227 SQRE1 = 0 228* 229* Update singular vectors if desired. 230* 231 IF( NCVT.GT.0 ) THEN 232 OPS = OPS + DBLE( 6*( NP1-1 )*NCVT ) 233 CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ), 234 $ WORK( NP1 ), VT, LDVT ) 235 END IF 236 END IF 237* 238* If matrix lower bidiagonal, rotate to be upper bidiagonal 239* by applying Givens rotations on the left. 240* 241 IF( IUPLO.EQ.2 ) THEN 242 OPS = OPS + DBLE( 8*( N-1 ) ) 243 DO 20 I = 1, N - 1 244 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 245 D( I ) = R 246 E( I ) = SN*D( I+1 ) 247 D( I+1 ) = CS*D( I+1 ) 248 IF( ROTATE ) THEN 249 WORK( I ) = CS 250 WORK( N+I ) = SN 251 END IF 252 20 CONTINUE 253* 254* If matrix (N+1)-by-N lower bidiagonal, one additional 255* rotation is needed. 256* 257 IF( SQRE1.EQ.1 ) THEN 258 OPS = OPS + DBLE( 6 ) 259 CALL DLARTG( D( N ), E( N ), CS, SN, R ) 260 D( N ) = R 261 IF( ROTATE ) THEN 262 WORK( N ) = CS 263 WORK( N+N ) = SN 264 END IF 265 END IF 266* 267* Update singular vectors if desired. 268* 269 IF( NRU.GT.0 ) THEN 270 IF( SQRE1.EQ.0 ) THEN 271 OPS = OPS + DBLE( 6*( N-1 )*NRU ) 272 CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), 273 $ WORK( NP1 ), U, LDU ) 274 ELSE 275 OPS = OPS + DBLE( 6*N*NRU ) 276 CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ), 277 $ WORK( NP1 ), U, LDU ) 278 END IF 279 END IF 280 IF( NCC.GT.0 ) THEN 281 IF( SQRE1.EQ.0 ) THEN 282 OPS = OPS + DBLE( 6*( N-1 )*NCC ) 283 CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), 284 $ WORK( NP1 ), C, LDC ) 285 ELSE 286 OPS = OPS + DBLE( 6*N*NCC ) 287 CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ), 288 $ WORK( NP1 ), C, LDC ) 289 END IF 290 END IF 291 END IF 292* 293* Call DBDSQR to compute the SVD of the reduced real 294* N-by-N upper bidiagonal matrix. 295* 296 CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, 297 $ LDC, WORK, INFO ) 298* 299* Sort the singular values into ascending order (insertion sort on 300* singular values, but only one transposition per singular vector) 301* 302 DO 40 I = 1, N 303* 304* Scan for smallest D(I). 305* 306 ISUB = I 307 SMIN = D( I ) 308 DO 30 J = I + 1, N 309 IF( D( J ).LT.SMIN ) THEN 310 ISUB = J 311 SMIN = D( J ) 312 END IF 313 30 CONTINUE 314 IF( ISUB.NE.I ) THEN 315* 316* Swap singular values and vectors. 317* 318 D( ISUB ) = D( I ) 319 D( I ) = SMIN 320 IF( NCVT.GT.0 ) 321 $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT ) 322 IF( NRU.GT.0 ) 323 $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 ) 324 IF( NCC.GT.0 ) 325 $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC ) 326 END IF 327 40 CONTINUE 328* 329 RETURN 330* 331* End of DLASDQ 332* 333 END 334