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