1 SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, 2 $ LWORK, RWORK, IWORK, INFO ) 3* 4* -- LAPACK driver 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 JOBZ 11 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 12* .. 13* .. Array Arguments .. 14 INTEGER IWORK( * ) 15 REAL RWORK( * ), S( * ) 16 COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 17 $ WORK( * ) 18* .. 19* 20* Purpose 21* ======= 22* 23* CGESDD computes the singular value decomposition (SVD) of a complex 24* M-by-N matrix A, optionally computing the left and/or right singular 25* vectors, by using divide-and-conquer method. The SVD is written 26* 27* A = U * SIGMA * conjugate-transpose(V) 28* 29* where SIGMA is an M-by-N matrix which is zero except for its 30* min(m,n) diagonal elements, U is an M-by-M unitary matrix, and 31* V is an N-by-N unitary matrix. The diagonal elements of SIGMA 32* are the singular values of A; they are real and non-negative, and 33* are returned in descending order. The first min(m,n) columns of 34* U and V are the left and right singular vectors of A. 35* 36* Note that the routine returns VT = V**H, not V. 37* 38* The divide and conquer algorithm makes very mild assumptions about 39* floating point arithmetic. It will work on machines with a guard 40* digit in add/subtract, or on those binary machines without guard 41* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 42* Cray-2. It could conceivably fail on hexadecimal or decimal machines 43* without guard digits, but we know of none. 44* 45* Arguments 46* ========= 47* 48* JOBZ (input) CHARACTER*1 49* Specifies options for computing all or part of the matrix U: 50* = 'A': all M columns of U and all N rows of V**H are 51* returned in the arrays U and VT; 52* = 'S': the first min(M,N) columns of U and the first 53* min(M,N) rows of V**H are returned in the arrays U 54* and VT; 55* = 'O': If M >= N, the first N columns of U are overwritten 56* on the array A and all rows of V**H are returned in 57* the array VT; 58* otherwise, all columns of U are returned in the 59* array U and the first M rows of V**H are overwritten 60* in the array VT; 61* = 'N': no columns of U or rows of V**H are computed. 62* 63* M (input) INTEGER 64* The number of rows of the input matrix A. M >= 0. 65* 66* N (input) INTEGER 67* The number of columns of the input matrix A. N >= 0. 68* 69* A (input/output) COMPLEX array, dimension (LDA,N) 70* On entry, the M-by-N matrix A. 71* On exit, 72* if JOBZ = 'O', A is overwritten with the first N columns 73* of U (the left singular vectors, stored 74* columnwise) if M >= N; 75* A is overwritten with the first M rows 76* of V**H (the right singular vectors, stored 77* rowwise) otherwise. 78* if JOBZ .ne. 'O', the contents of A are destroyed. 79* 80* LDA (input) INTEGER 81* The leading dimension of the array A. LDA >= max(1,M). 82* 83* S (output) REAL array, dimension (min(M,N)) 84* The singular values of A, sorted so that S(i) >= S(i+1). 85* 86* U (output) COMPLEX array, dimension (LDU,UCOL) 87* UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; 88* UCOL = min(M,N) if JOBZ = 'S'. 89* If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M 90* unitary matrix U; 91* if JOBZ = 'S', U contains the first min(M,N) columns of U 92* (the left singular vectors, stored columnwise); 93* if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. 94* 95* LDU (input) INTEGER 96* The leading dimension of the array U. LDU >= 1; if 97* JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. 98* 99* VT (output) COMPLEX array, dimension (LDVT,N) 100* If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the 101* N-by-N unitary matrix V**H; 102* if JOBZ = 'S', VT contains the first min(M,N) rows of 103* V**H (the right singular vectors, stored rowwise); 104* if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. 105* 106* LDVT (input) INTEGER 107* The leading dimension of the array VT. LDVT >= 1; if 108* JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; 109* if JOBZ = 'S', LDVT >= min(M,N). 110* 111* WORK (workspace/output) COMPLEX array, dimension (LWORK) 112* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 113* 114* LWORK (input) INTEGER 115* The dimension of the array WORK. LWORK >= 1. 116* if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N). 117* if JOBZ = 'O', 118* LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N). 119* if JOBZ = 'S' or 'A', 120* LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N). 121* For good performance, LWORK should generally be larger. 122* If LWORK < 0 but other input arguments are legal, WORK(1) 123* returns the optimal LWORK. 124* 125* RWORK (workspace) REAL array, dimension (LRWORK) 126* If JOBZ = 'N', LRWORK >= 7*min(M,N). 127* Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 5*min(M,N) 128* 129* IWORK (workspace) INTEGER array, dimension (8*min(M,N)) 130* 131* INFO (output) INTEGER 132* = 0: successful exit. 133* < 0: if INFO = -i, the i-th argument had an illegal value. 134* > 0: The updating process of SBDSDC did not converge. 135* 136* Further Details 137* =============== 138* 139* Based on contributions by 140* Ming Gu and Huan Ren, Computer Science Division, University of 141* California at Berkeley, USA 142* 143* ===================================================================== 144* 145* .. Parameters .. 146 COMPLEX CZERO, CONE 147 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 148 $ CONE = ( 1.0E0, 0.0E0 ) ) 149 REAL ZERO, ONE 150 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 151* .. 152* .. Local Scalars .. 153 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS 154 INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT, 155 $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, 156 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, 157 $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL 158 REAL ANRM, BIGNUM, EPS, SMLNUM 159* .. 160* .. Local Arrays .. 161 INTEGER IDUM( 1 ) 162 REAL DUM( 1 ) 163* .. 164* .. External Subroutines .. 165 EXTERNAL CGEBRD, CGELQF, CGEMM, CGEQRF, CLACP2, CLACPY, 166 $ CLACRM, CLARCM, CLASCL, CLASET, CUNGBR, CUNGLQ, 167 $ CUNGQR, CUNMBR, SBDSDC, SLASCL, XERBLA 168* .. 169* .. External Functions .. 170 LOGICAL LSAME 171 INTEGER ILAENV 172 REAL CLANGE, SLAMCH 173 EXTERNAL CLANGE, ILAENV, LSAME, SLAMCH 174* .. 175* .. Intrinsic Functions .. 176 INTRINSIC INT, MAX, MIN, SQRT 177* .. 178* .. Executable Statements .. 179* 180* Test the input arguments 181* 182 INFO = 0 183 MINMN = MIN( M, N ) 184 MNTHR1 = INT( MINMN*17.0E0 / 9.0E0 ) 185 MNTHR2 = INT( MINMN*5.0E0 / 3.0E0 ) 186 WNTQA = LSAME( JOBZ, 'A' ) 187 WNTQS = LSAME( JOBZ, 'S' ) 188 WNTQAS = WNTQA .OR. WNTQS 189 WNTQO = LSAME( JOBZ, 'O' ) 190 WNTQN = LSAME( JOBZ, 'N' ) 191 MINWRK = 1 192 MAXWRK = 1 193 LQUERY = ( LWORK.EQ.-1 ) 194* 195 IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN 196 INFO = -1 197 ELSE IF( M.LT.0 ) THEN 198 INFO = -2 199 ELSE IF( N.LT.0 ) THEN 200 INFO = -3 201 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 202 INFO = -5 203 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR. 204 $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN 205 INFO = -8 206 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR. 207 $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR. 208 $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN 209 INFO = -10 210 END IF 211* 212* Compute workspace 213* (Note: Comments in the code beginning "Workspace:" describe the 214* minimal amount of workspace needed at that point in the code, 215* as well as the preferred amount for good performance. 216* CWorkspace refers to complex workspace, and RWorkspace to 217* real workspace. NB refers to the optimal block size for the 218* immediately following subroutine, as returned by ILAENV.) 219* 220 IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN 221 IF( M.GE.N ) THEN 222* 223* There is no complex work space needed for bidiagonal SVD 224* The real work space needed for bidiagonal SVD is BDSPAC, 225* BDSPAC = 3*N*N + 4*N 226* 227 IF( M.GE.MNTHR1 ) THEN 228 IF( WNTQN ) THEN 229* 230* Path 1 (M much larger than N, JOBZ='N') 231* 232 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, 233 $ -1 ) 234 WRKBL = MAX( WRKBL, 2*N+2*N* 235 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 236 MAXWRK = WRKBL 237 MINWRK = 3*N 238 ELSE IF( WNTQO ) THEN 239* 240* Path 2 (M much larger than N, JOBZ='O') 241* 242 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 243 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, 244 $ N, N, -1 ) ) 245 WRKBL = MAX( WRKBL, 2*N+2*N* 246 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 247 WRKBL = MAX( WRKBL, 2*N+N* 248 $ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) ) 249 WRKBL = MAX( WRKBL, 2*N+N* 250 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 251 MAXWRK = M*N + N*N + WRKBL 252 MINWRK = 2*N*N + 3*N 253 ELSE IF( WNTQS ) THEN 254* 255* Path 3 (M much larger than N, JOBZ='S') 256* 257 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 258 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M, 259 $ N, N, -1 ) ) 260 WRKBL = MAX( WRKBL, 2*N+2*N* 261 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 262 WRKBL = MAX( WRKBL, 2*N+N* 263 $ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) ) 264 WRKBL = MAX( WRKBL, 2*N+N* 265 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 266 MAXWRK = N*N + WRKBL 267 MINWRK = N*N + 3*N 268 ELSE IF( WNTQA ) THEN 269* 270* Path 4 (M much larger than N, JOBZ='A') 271* 272 WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 273 WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M, 274 $ M, N, -1 ) ) 275 WRKBL = MAX( WRKBL, 2*N+2*N* 276 $ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) ) 277 WRKBL = MAX( WRKBL, 2*N+N* 278 $ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) ) 279 WRKBL = MAX( WRKBL, 2*N+N* 280 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 281 MAXWRK = N*N + WRKBL 282 MINWRK = N*N + 2*N + M 283 END IF 284 ELSE IF( M.GE.MNTHR2 ) THEN 285* 286* Path 5 (M much larger than N, but not as much as MNTHR1) 287* 288 MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 289 $ -1, -1 ) 290 MINWRK = 2*N + M 291 IF( WNTQO ) THEN 292 MAXWRK = MAX( MAXWRK, 2*N+N* 293 $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) 294 MAXWRK = MAX( MAXWRK, 2*N+N* 295 $ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) ) 296 MAXWRK = MAXWRK + M*N 297 MINWRK = MINWRK + N*N 298 ELSE IF( WNTQS ) THEN 299 MAXWRK = MAX( MAXWRK, 2*N+N* 300 $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) 301 MAXWRK = MAX( MAXWRK, 2*N+N* 302 $ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) ) 303 ELSE IF( WNTQA ) THEN 304 MAXWRK = MAX( MAXWRK, 2*N+N* 305 $ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) ) 306 MAXWRK = MAX( MAXWRK, 2*N+M* 307 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 308 END IF 309 ELSE 310* 311* Path 6 (M at least N, but not much larger) 312* 313 MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 314 $ -1, -1 ) 315 MINWRK = 2*N + M 316 IF( WNTQO ) THEN 317 MAXWRK = MAX( MAXWRK, 2*N+N* 318 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 319 MAXWRK = MAX( MAXWRK, 2*N+N* 320 $ ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) ) 321 MAXWRK = MAXWRK + M*N 322 MINWRK = MINWRK + N*N 323 ELSE IF( WNTQS ) THEN 324 MAXWRK = MAX( MAXWRK, 2*N+N* 325 $ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) ) 326 MAXWRK = MAX( MAXWRK, 2*N+N* 327 $ ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) ) 328 ELSE IF( WNTQA ) THEN 329 MAXWRK = MAX( MAXWRK, 2*N+N* 330 $ ILAENV( 1, 'CUNGBR', 'PRC', N, N, N, -1 ) ) 331 MAXWRK = MAX( MAXWRK, 2*N+M* 332 $ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) ) 333 END IF 334 END IF 335 ELSE 336* 337* There is no complex work space needed for bidiagonal SVD 338* The real work space needed for bidiagonal SVD is BDSPAC, 339* BDSPAC = 3*M*M + 4*M 340* 341 IF( N.GE.MNTHR1 ) THEN 342 IF( WNTQN ) THEN 343* 344* Path 1t (N much larger than M, JOBZ='N') 345* 346 MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, 347 $ -1 ) 348 MAXWRK = MAX( MAXWRK, 2*M+2*M* 349 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 350 MINWRK = 3*M 351 ELSE IF( WNTQO ) THEN 352* 353* Path 2t (N much larger than M, JOBZ='O') 354* 355 WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) 356 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, 357 $ N, M, -1 ) ) 358 WRKBL = MAX( WRKBL, 2*M+2*M* 359 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 360 WRKBL = MAX( WRKBL, 2*M+M* 361 $ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) ) 362 WRKBL = MAX( WRKBL, 2*M+M* 363 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) ) 364 MAXWRK = M*N + M*M + WRKBL 365 MINWRK = 2*M*M + 3*M 366 ELSE IF( WNTQS ) THEN 367* 368* Path 3t (N much larger than M, JOBZ='S') 369* 370 WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) 371 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M, 372 $ N, M, -1 ) ) 373 WRKBL = MAX( WRKBL, 2*M+2*M* 374 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 375 WRKBL = MAX( WRKBL, 2*M+M* 376 $ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) ) 377 WRKBL = MAX( WRKBL, 2*M+M* 378 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) ) 379 MAXWRK = M*M + WRKBL 380 MINWRK = M*M + 3*M 381 ELSE IF( WNTQA ) THEN 382* 383* Path 4t (N much larger than M, JOBZ='A') 384* 385 WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 ) 386 WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N, 387 $ N, M, -1 ) ) 388 WRKBL = MAX( WRKBL, 2*M+2*M* 389 $ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) ) 390 WRKBL = MAX( WRKBL, 2*M+M* 391 $ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) ) 392 WRKBL = MAX( WRKBL, 2*M+M* 393 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) ) 394 MAXWRK = M*M + WRKBL 395 MINWRK = M*M + 2*M + N 396 END IF 397 ELSE IF( N.GE.MNTHR2 ) THEN 398* 399* Path 5t (N much larger than M, but not as much as MNTHR1) 400* 401 MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 402 $ -1, -1 ) 403 MINWRK = 2*M + N 404 IF( WNTQO ) THEN 405 MAXWRK = MAX( MAXWRK, 2*M+M* 406 $ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) ) 407 MAXWRK = MAX( MAXWRK, 2*M+M* 408 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 409 MAXWRK = MAXWRK + M*N 410 MINWRK = MINWRK + M*M 411 ELSE IF( WNTQS ) THEN 412 MAXWRK = MAX( MAXWRK, 2*M+M* 413 $ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) ) 414 MAXWRK = MAX( MAXWRK, 2*M+M* 415 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 416 ELSE IF( WNTQA ) THEN 417 MAXWRK = MAX( MAXWRK, 2*M+N* 418 $ ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) ) 419 MAXWRK = MAX( MAXWRK, 2*M+M* 420 $ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) ) 421 END IF 422 ELSE 423* 424* Path 6t (N greater than M, but not much larger) 425* 426 MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N, 427 $ -1, -1 ) 428 MINWRK = 2*M + N 429 IF( WNTQO ) THEN 430 MAXWRK = MAX( MAXWRK, 2*M+M* 431 $ ILAENV( 1, 'CUNMBR', 'PRC', M, N, M, -1 ) ) 432 MAXWRK = MAX( MAXWRK, 2*M+M* 433 $ ILAENV( 1, 'CUNMBR', 'QLN', M, M, N, -1 ) ) 434 MAXWRK = MAXWRK + M*N 435 MINWRK = MINWRK + M*M 436 ELSE IF( WNTQS ) THEN 437 MAXWRK = MAX( MAXWRK, 2*M+M* 438 $ ILAENV( 1, 'CUNGBR', 'PRC', M, N, M, -1 ) ) 439 MAXWRK = MAX( MAXWRK, 2*M+M* 440 $ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) ) 441 ELSE IF( WNTQA ) THEN 442 MAXWRK = MAX( MAXWRK, 2*M+N* 443 $ ILAENV( 1, 'CUNGBR', 'PRC', N, N, M, -1 ) ) 444 MAXWRK = MAX( MAXWRK, 2*M+M* 445 $ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) ) 446 END IF 447 END IF 448 END IF 449 MAXWRK = MAX( MAXWRK, MINWRK ) 450 WORK( 1 ) = MAXWRK 451 END IF 452* 453 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 454 INFO = -13 455 END IF 456 IF( INFO.NE.0 ) THEN 457 CALL XERBLA( 'CGESDD', -INFO ) 458 RETURN 459 ELSE IF( LQUERY ) THEN 460 RETURN 461 END IF 462* 463* Quick return if possible 464* 465 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 466 IF( LWORK.GE.1 ) 467 $ WORK( 1 ) = ONE 468 RETURN 469 END IF 470* 471* Get machine constants 472* 473 EPS = SLAMCH( 'P' ) 474 SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS 475 BIGNUM = ONE / SMLNUM 476* 477* Scale A if max element outside range [SMLNUM,BIGNUM] 478* 479 ANRM = CLANGE( 'M', M, N, A, LDA, DUM ) 480 ISCL = 0 481 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 482 ISCL = 1 483 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 484 ELSE IF( ANRM.GT.BIGNUM ) THEN 485 ISCL = 1 486 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 487 END IF 488* 489 IF( M.GE.N ) THEN 490* 491* A has at least as many rows as columns. If A has sufficiently 492* more rows than columns, first reduce using the QR 493* decomposition (if sufficient workspace available) 494* 495 IF( M.GE.MNTHR1 ) THEN 496* 497 IF( WNTQN ) THEN 498* 499* Path 1 (M much larger than N, JOBZ='N') 500* No singular vectors to be computed 501* 502 ITAU = 1 503 NWORK = ITAU + N 504* 505* Compute A=Q*R 506* (CWorkspace: need 2*N, prefer N+N*NB) 507* (RWorkspace: need 0) 508* 509 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 510 $ LWORK-NWORK+1, IERR ) 511* 512* Zero out below R 513* 514 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 515 $ LDA ) 516 IE = 1 517 ITAUQ = 1 518 ITAUP = ITAUQ + N 519 NWORK = ITAUP + N 520* 521* Bidiagonalize R in A 522* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 523* (RWorkspace: need N) 524* 525 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 526 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 527 $ IERR ) 528 NRWORK = IE + N 529* 530* Perform bidiagonal SVD, compute singular values only 531* (CWorkspace: 0) 532* (RWorkspace: need BDSPAC) 533* 534 CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 535 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 536* 537 ELSE IF( WNTQO ) THEN 538* 539* Path 2 (M much larger than N, JOBZ='O') 540* N left singular vectors to be overwritten on A and 541* N right singular vectors to be computed in VT 542* 543 IU = 1 544* 545* WORK(IU) is N by N 546* 547 LDWRKU = N 548 IR = IU + LDWRKU*N 549 IF( LWORK.GE.M*N+N*N+3*N ) THEN 550* 551* WORK(IR) is M by N 552* 553 LDWRKR = M 554 ELSE 555 LDWRKR = ( LWORK-N*N-3*N ) / N 556 END IF 557 ITAU = IR + LDWRKR*N 558 NWORK = ITAU + N 559* 560* Compute A=Q*R 561* (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB) 562* (RWorkspace: 0) 563* 564 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 565 $ LWORK-NWORK+1, IERR ) 566* 567* Copy R to WORK( IR ), zeroing out below it 568* 569 CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 570 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ), 571 $ LDWRKR ) 572* 573* Generate Q in A 574* (CWorkspace: need 2*N, prefer N+N*NB) 575* (RWorkspace: 0) 576* 577 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 578 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 579 IE = 1 580 ITAUQ = ITAU 581 ITAUP = ITAUQ + N 582 NWORK = ITAUP + N 583* 584* Bidiagonalize R in WORK(IR) 585* (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB) 586* (RWorkspace: need N) 587* 588 CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 589 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 590 $ LWORK-NWORK+1, IERR ) 591* 592* Perform bidiagonal SVD, computing left singular vectors 593* of R in WORK(IRU) and computing right singular vectors 594* of R in WORK(IRVT) 595* (CWorkspace: need 0) 596* (RWorkspace: need BDSPAC) 597* 598 IRU = IE + N 599 IRVT = IRU + N*N 600 NRWORK = IRVT + N*N 601 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 602 $ N, RWORK( IRVT ), N, DUM, IDUM, 603 $ RWORK( NRWORK ), IWORK, INFO ) 604* 605* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 606* Overwrite WORK(IU) by the left singular vectors of R 607* (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB) 608* (RWorkspace: 0) 609* 610 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 611 $ LDWRKU ) 612 CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 613 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 614 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 615* 616* Copy real matrix RWORK(IRVT) to complex matrix VT 617* Overwrite VT by the right singular vectors of R 618* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 619* (RWorkspace: 0) 620* 621 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 622 CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, 623 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 624 $ LWORK-NWORK+1, IERR ) 625* 626* Multiply Q in A by left singular vectors of R in 627* WORK(IU), storing result in WORK(IR) and copying to A 628* (CWorkspace: need 2*N*N, prefer N*N+M*N) 629* (RWorkspace: 0) 630* 631 DO 10 I = 1, M, LDWRKR 632 CHUNK = MIN( M-I+1, LDWRKR ) 633 CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 634 $ LDA, WORK( IU ), LDWRKU, CZERO, 635 $ WORK( IR ), LDWRKR ) 636 CALL CLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 637 $ A( I, 1 ), LDA ) 638 10 CONTINUE 639* 640 ELSE IF( WNTQS ) THEN 641* 642* Path 3 (M much larger than N, JOBZ='S') 643* N left singular vectors to be computed in U and 644* N right singular vectors to be computed in VT 645* 646 IR = 1 647* 648* WORK(IR) is N by N 649* 650 LDWRKR = N 651 ITAU = IR + LDWRKR*N 652 NWORK = ITAU + N 653* 654* Compute A=Q*R 655* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 656* (RWorkspace: 0) 657* 658 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 659 $ LWORK-NWORK+1, IERR ) 660* 661* Copy R to WORK(IR), zeroing out below it 662* 663 CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 664 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ), 665 $ LDWRKR ) 666* 667* Generate Q in A 668* (CWorkspace: need 2*N, prefer N+N*NB) 669* (RWorkspace: 0) 670* 671 CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ), 672 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 673 IE = 1 674 ITAUQ = ITAU 675 ITAUP = ITAUQ + N 676 NWORK = ITAUP + N 677* 678* Bidiagonalize R in WORK(IR) 679* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 680* (RWorkspace: need N) 681* 682 CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 683 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 684 $ LWORK-NWORK+1, IERR ) 685* 686* Perform bidiagonal SVD, computing left singular vectors 687* of bidiagonal matrix in RWORK(IRU) and computing right 688* singular vectors of bidiagonal matrix in RWORK(IRVT) 689* (CWorkspace: need 0) 690* (RWorkspace: need BDSPAC) 691* 692 IRU = IE + N 693 IRVT = IRU + N*N 694 NRWORK = IRVT + N*N 695 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 696 $ N, RWORK( IRVT ), N, DUM, IDUM, 697 $ RWORK( NRWORK ), IWORK, INFO ) 698* 699* Copy real matrix RWORK(IRU) to complex matrix U 700* Overwrite U by left singular vectors of R 701* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 702* (RWorkspace: 0) 703* 704 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 705 CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 706 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 707 $ LWORK-NWORK+1, IERR ) 708* 709* Copy real matrix RWORK(IRVT) to complex matrix VT 710* Overwrite VT by right singular vectors of R 711* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 712* (RWorkspace: 0) 713* 714 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 715 CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, 716 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 717 $ LWORK-NWORK+1, IERR ) 718* 719* Multiply Q in A by left singular vectors of R in 720* WORK(IR), storing result in U 721* (CWorkspace: need N*N) 722* (RWorkspace: 0) 723* 724 CALL CLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) 725 CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ), 726 $ LDWRKR, CZERO, U, LDU ) 727* 728 ELSE IF( WNTQA ) THEN 729* 730* Path 4 (M much larger than N, JOBZ='A') 731* M left singular vectors to be computed in U and 732* N right singular vectors to be computed in VT 733* 734 IU = 1 735* 736* WORK(IU) is N by N 737* 738 LDWRKU = N 739 ITAU = IU + LDWRKU*N 740 NWORK = ITAU + N 741* 742* Compute A=Q*R, copying result to U 743* (CWorkspace: need 2*N, prefer N+N*NB) 744* (RWorkspace: 0) 745* 746 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 747 $ LWORK-NWORK+1, IERR ) 748 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 749* 750* Generate Q in U 751* (CWorkspace: need N+M, prefer N+M*NB) 752* (RWorkspace: 0) 753* 754 CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ), 755 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 756* 757* Produce R in A, zeroing out below it 758* 759 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 760 $ LDA ) 761 IE = 1 762 ITAUQ = ITAU 763 ITAUP = ITAUQ + N 764 NWORK = ITAUP + N 765* 766* Bidiagonalize R in A 767* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 768* (RWorkspace: need N) 769* 770 CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 771 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 772 $ IERR ) 773 IRU = IE + N 774 IRVT = IRU + N*N 775 NRWORK = IRVT + N*N 776* 777* Perform bidiagonal SVD, computing left singular vectors 778* of bidiagonal matrix in RWORK(IRU) and computing right 779* singular vectors of bidiagonal matrix in RWORK(IRVT) 780* (CWorkspace: need 0) 781* (RWorkspace: need BDSPAC) 782* 783 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 784 $ N, RWORK( IRVT ), N, DUM, IDUM, 785 $ RWORK( NRWORK ), IWORK, INFO ) 786* 787* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 788* Overwrite WORK(IU) by left singular vectors of R 789* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 790* (RWorkspace: 0) 791* 792 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 793 $ LDWRKU ) 794 CALL CUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA, 795 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 796 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 797* 798* Copy real matrix RWORK(IRVT) to complex matrix VT 799* Overwrite VT by right singular vectors of R 800* (CWorkspace: need 3*N, prefer 2*N+N*NB) 801* (RWorkspace: 0) 802* 803 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 804 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 805 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 806 $ LWORK-NWORK+1, IERR ) 807* 808* Multiply Q in U by left singular vectors of R in 809* WORK(IU), storing result in A 810* (CWorkspace: need N*N) 811* (RWorkspace: 0) 812* 813 CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ), 814 $ LDWRKU, CZERO, A, LDA ) 815* 816* Copy left singular vectors of A from A to U 817* 818 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 819* 820 END IF 821* 822 ELSE IF( M.GE.MNTHR2 ) THEN 823* 824* MNTHR2 <= M < MNTHR1 825* 826* Path 5 (M much larger than N, but not as much as MNTHR1) 827* Reduce to bidiagonal form without QR decomposition, use 828* CUNGBR and matrix multiplication to compute singular vectors 829* 830 IE = 1 831 NRWORK = IE + N 832 ITAUQ = 1 833 ITAUP = ITAUQ + N 834 NWORK = ITAUP + N 835* 836* Bidiagonalize A 837* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 838* (RWorkspace: need N) 839* 840 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 841 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 842 $ IERR ) 843 IF( WNTQN ) THEN 844* 845* Compute singular values only 846* (Cworkspace: 0) 847* (Rworkspace: need BDSPAC) 848* 849 CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 850 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 851 ELSE IF( WNTQO ) THEN 852 IU = NWORK 853 IRU = NRWORK 854 IRVT = IRU + N*N 855 NRWORK = IRVT + N*N 856* 857* Copy A to VT, generate P**H 858* (Cworkspace: need 2*N, prefer N+N*NB) 859* (Rworkspace: 0) 860* 861 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 862 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 863 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 864* 865* Generate Q in A 866* (CWorkspace: need 2*N, prefer N+N*NB) 867* (RWorkspace: 0) 868* 869 CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 870 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 871* 872 IF( LWORK.GE.M*N+3*N ) THEN 873* 874* WORK( IU ) is M by N 875* 876 LDWRKU = M 877 ELSE 878* 879* WORK(IU) is LDWRKU by N 880* 881 LDWRKU = ( LWORK-3*N ) / N 882 END IF 883 NWORK = IU + LDWRKU*N 884* 885* Perform bidiagonal SVD, computing left singular vectors 886* of bidiagonal matrix in RWORK(IRU) and computing right 887* singular vectors of bidiagonal matrix in RWORK(IRVT) 888* (CWorkspace: need 0) 889* (RWorkspace: need BDSPAC) 890* 891 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 892 $ N, RWORK( IRVT ), N, DUM, IDUM, 893 $ RWORK( NRWORK ), IWORK, INFO ) 894* 895* Multiply real matrix RWORK(IRVT) by P**H in VT, 896* storing the result in WORK(IU), copying to VT 897* (Cworkspace: need 0) 898* (Rworkspace: need 3*N*N) 899* 900 CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, 901 $ WORK( IU ), LDWRKU, RWORK( NRWORK ) ) 902 CALL CLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT ) 903* 904* Multiply Q in A by real matrix RWORK(IRU), storing the 905* result in WORK(IU), copying to A 906* (CWorkspace: need N*N, prefer M*N) 907* (Rworkspace: need 3*N*N, prefer N*N+2*M*N) 908* 909 NRWORK = IRVT 910 DO 20 I = 1, M, LDWRKU 911 CHUNK = MIN( M-I+1, LDWRKU ) 912 CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ), 913 $ N, WORK( IU ), LDWRKU, RWORK( NRWORK ) ) 914 CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 915 $ A( I, 1 ), LDA ) 916 20 CONTINUE 917* 918 ELSE IF( WNTQS ) THEN 919* 920* Copy A to VT, generate P**H 921* (Cworkspace: need 2*N, prefer N+N*NB) 922* (Rworkspace: 0) 923* 924 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 925 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 926 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 927* 928* Copy A to U, generate Q 929* (Cworkspace: need 2*N, prefer N+N*NB) 930* (Rworkspace: 0) 931* 932 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 933 CALL CUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ), 934 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 935* 936* Perform bidiagonal SVD, computing left singular vectors 937* of bidiagonal matrix in RWORK(IRU) and computing right 938* singular vectors of bidiagonal matrix in RWORK(IRVT) 939* (CWorkspace: need 0) 940* (RWorkspace: need BDSPAC) 941* 942 IRU = NRWORK 943 IRVT = IRU + N*N 944 NRWORK = IRVT + N*N 945 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 946 $ N, RWORK( IRVT ), N, DUM, IDUM, 947 $ RWORK( NRWORK ), IWORK, INFO ) 948* 949* Multiply real matrix RWORK(IRVT) by P**H in VT, 950* storing the result in A, copying to VT 951* (Cworkspace: need 0) 952* (Rworkspace: need 3*N*N) 953* 954 CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, 955 $ RWORK( NRWORK ) ) 956 CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT ) 957* 958* Multiply Q in U by real matrix RWORK(IRU), storing the 959* result in A, copying to U 960* (CWorkspace: need 0) 961* (Rworkspace: need N*N+2*M*N) 962* 963 NRWORK = IRVT 964 CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, 965 $ RWORK( NRWORK ) ) 966 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 967 ELSE 968* 969* Copy A to VT, generate P**H 970* (Cworkspace: need 2*N, prefer N+N*NB) 971* (Rworkspace: 0) 972* 973 CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT ) 974 CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 975 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 976* 977* Copy A to U, generate Q 978* (Cworkspace: need 2*N, prefer N+N*NB) 979* (Rworkspace: 0) 980* 981 CALL CLACPY( 'L', M, N, A, LDA, U, LDU ) 982 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 983 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 984* 985* Perform bidiagonal SVD, computing left singular vectors 986* of bidiagonal matrix in RWORK(IRU) and computing right 987* singular vectors of bidiagonal matrix in RWORK(IRVT) 988* (CWorkspace: need 0) 989* (RWorkspace: need BDSPAC) 990* 991 IRU = NRWORK 992 IRVT = IRU + N*N 993 NRWORK = IRVT + N*N 994 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 995 $ N, RWORK( IRVT ), N, DUM, IDUM, 996 $ RWORK( NRWORK ), IWORK, INFO ) 997* 998* Multiply real matrix RWORK(IRVT) by P**H in VT, 999* storing the result in A, copying to VT 1000* (Cworkspace: need 0) 1001* (Rworkspace: need 3*N*N) 1002* 1003 CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, 1004 $ RWORK( NRWORK ) ) 1005 CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT ) 1006* 1007* Multiply Q in U by real matrix RWORK(IRU), storing the 1008* result in A, copying to U 1009* (CWorkspace: 0) 1010* (Rworkspace: need 3*N*N) 1011* 1012 NRWORK = IRVT 1013 CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, 1014 $ RWORK( NRWORK ) ) 1015 CALL CLACPY( 'F', M, N, A, LDA, U, LDU ) 1016 END IF 1017* 1018 ELSE 1019* 1020* M .LT. MNTHR2 1021* 1022* Path 6 (M at least N, but not much larger) 1023* Reduce to bidiagonal form without QR decomposition 1024* Use CUNMBR to compute singular vectors 1025* 1026 IE = 1 1027 NRWORK = IE + N 1028 ITAUQ = 1 1029 ITAUP = ITAUQ + N 1030 NWORK = ITAUP + N 1031* 1032* Bidiagonalize A 1033* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 1034* (RWorkspace: need N) 1035* 1036 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 1037 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1038 $ IERR ) 1039 IF( WNTQN ) THEN 1040* 1041* Compute singular values only 1042* (Cworkspace: 0) 1043* (Rworkspace: need BDSPAC) 1044* 1045 CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 1046 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 1047 ELSE IF( WNTQO ) THEN 1048 IU = NWORK 1049 IRU = NRWORK 1050 IRVT = IRU + N*N 1051 NRWORK = IRVT + N*N 1052 IF( LWORK.GE.M*N+3*N ) THEN 1053* 1054* WORK( IU ) is M by N 1055* 1056 LDWRKU = M 1057 ELSE 1058* 1059* WORK( IU ) is LDWRKU by N 1060* 1061 LDWRKU = ( LWORK-3*N ) / N 1062 END IF 1063 NWORK = IU + LDWRKU*N 1064* 1065* Perform bidiagonal SVD, computing left singular vectors 1066* of bidiagonal matrix in RWORK(IRU) and computing right 1067* singular vectors of bidiagonal matrix in RWORK(IRVT) 1068* (CWorkspace: need 0) 1069* (RWorkspace: need BDSPAC) 1070* 1071 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 1072 $ N, RWORK( IRVT ), N, DUM, IDUM, 1073 $ RWORK( NRWORK ), IWORK, INFO ) 1074* 1075* Copy real matrix RWORK(IRVT) to complex matrix VT 1076* Overwrite VT by right singular vectors of A 1077* (Cworkspace: need 2*N, prefer N+N*NB) 1078* (Rworkspace: need 0) 1079* 1080 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 1081 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 1082 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1083 $ LWORK-NWORK+1, IERR ) 1084* 1085 IF( LWORK.GE.M*N+3*N ) THEN 1086* 1087* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 1088* Overwrite WORK(IU) by left singular vectors of A, copying 1089* to A 1090* (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB) 1091* (Rworkspace: need 0) 1092* 1093 CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IU ), 1094 $ LDWRKU ) 1095 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 1096 $ LDWRKU ) 1097 CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 1098 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 1099 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1100 CALL CLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) 1101 ELSE 1102* 1103* Generate Q in A 1104* (Cworkspace: need 2*N, prefer N+N*NB) 1105* (Rworkspace: need 0) 1106* 1107 CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 1108 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1109* 1110* Multiply Q in A by real matrix RWORK(IRU), storing the 1111* result in WORK(IU), copying to A 1112* (CWorkspace: need N*N, prefer M*N) 1113* (Rworkspace: need 3*N*N, prefer N*N+2*M*N) 1114* 1115 NRWORK = IRVT 1116 DO 30 I = 1, M, LDWRKU 1117 CHUNK = MIN( M-I+1, LDWRKU ) 1118 CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, 1119 $ RWORK( IRU ), N, WORK( IU ), LDWRKU, 1120 $ RWORK( NRWORK ) ) 1121 CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 1122 $ A( I, 1 ), LDA ) 1123 30 CONTINUE 1124 END IF 1125* 1126 ELSE IF( WNTQS ) THEN 1127* 1128* Perform bidiagonal SVD, computing left singular vectors 1129* of bidiagonal matrix in RWORK(IRU) and computing right 1130* singular vectors of bidiagonal matrix in RWORK(IRVT) 1131* (CWorkspace: need 0) 1132* (RWorkspace: need BDSPAC) 1133* 1134 IRU = NRWORK 1135 IRVT = IRU + N*N 1136 NRWORK = IRVT + N*N 1137 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 1138 $ N, RWORK( IRVT ), N, DUM, IDUM, 1139 $ RWORK( NRWORK ), IWORK, INFO ) 1140* 1141* Copy real matrix RWORK(IRU) to complex matrix U 1142* Overwrite U by left singular vectors of A 1143* (CWorkspace: need 3*N, prefer 2*N+N*NB) 1144* (RWorkspace: 0) 1145* 1146 CALL CLASET( 'F', M, N, CZERO, CZERO, U, LDU ) 1147 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 1148 CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 1149 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1150 $ LWORK-NWORK+1, IERR ) 1151* 1152* Copy real matrix RWORK(IRVT) to complex matrix VT 1153* Overwrite VT by right singular vectors of A 1154* (CWorkspace: need 3*N, prefer 2*N+N*NB) 1155* (RWorkspace: 0) 1156* 1157 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 1158 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 1159 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1160 $ LWORK-NWORK+1, IERR ) 1161 ELSE 1162* 1163* Perform bidiagonal SVD, computing left singular vectors 1164* of bidiagonal matrix in RWORK(IRU) and computing right 1165* singular vectors of bidiagonal matrix in RWORK(IRVT) 1166* (CWorkspace: need 0) 1167* (RWorkspace: need BDSPAC) 1168* 1169 IRU = NRWORK 1170 IRVT = IRU + N*N 1171 NRWORK = IRVT + N*N 1172 CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 1173 $ N, RWORK( IRVT ), N, DUM, IDUM, 1174 $ RWORK( NRWORK ), IWORK, INFO ) 1175* 1176* Set the right corner of U to identity matrix 1177* 1178 CALL CLASET( 'F', M, M, CZERO, CZERO, U, LDU ) 1179 CALL CLASET( 'F', M-N, M-N, CZERO, CONE, U( N+1, N+1 ), 1180 $ LDU ) 1181* 1182* Copy real matrix RWORK(IRU) to complex matrix U 1183* Overwrite U by left singular vectors of A 1184* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1185* (RWorkspace: 0) 1186* 1187 CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 1188 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1189 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1190 $ LWORK-NWORK+1, IERR ) 1191* 1192* Copy real matrix RWORK(IRVT) to complex matrix VT 1193* Overwrite VT by right singular vectors of A 1194* (CWorkspace: need 3*N, prefer 2*N+N*NB) 1195* (RWorkspace: 0) 1196* 1197 CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 1198 CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 1199 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1200 $ LWORK-NWORK+1, IERR ) 1201 END IF 1202* 1203 END IF 1204* 1205 ELSE 1206* 1207* A has more columns than rows. If A has sufficiently more 1208* columns than rows, first reduce using the LQ decomposition 1209* (if sufficient workspace available) 1210* 1211 IF( N.GE.MNTHR1 ) THEN 1212* 1213 IF( WNTQN ) THEN 1214* 1215* Path 1t (N much larger than M, JOBZ='N') 1216* No singular vectors to be computed 1217* 1218 ITAU = 1 1219 NWORK = ITAU + M 1220* 1221* Compute A=L*Q 1222* (CWorkspace: need 2*M, prefer M+M*NB) 1223* (RWorkspace: 0) 1224* 1225 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1226 $ LWORK-NWORK+1, IERR ) 1227* 1228* Zero out above L 1229* 1230 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 1231 $ LDA ) 1232 IE = 1 1233 ITAUQ = 1 1234 ITAUP = ITAUQ + M 1235 NWORK = ITAUP + M 1236* 1237* Bidiagonalize L in A 1238* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 1239* (RWorkspace: need M) 1240* 1241 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 1242 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1243 $ IERR ) 1244 NRWORK = IE + M 1245* 1246* Perform bidiagonal SVD, compute singular values only 1247* (CWorkspace: 0) 1248* (RWorkspace: need BDSPAC) 1249* 1250 CALL SBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 1251 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 1252* 1253 ELSE IF( WNTQO ) THEN 1254* 1255* Path 2t (N much larger than M, JOBZ='O') 1256* M right singular vectors to be overwritten on A and 1257* M left singular vectors to be computed in U 1258* 1259 IVT = 1 1260 LDWKVT = M 1261* 1262* WORK(IVT) is M by M 1263* 1264 IL = IVT + LDWKVT*M 1265 IF( LWORK.GE.M*N+M*M+3*M ) THEN 1266* 1267* WORK(IL) M by N 1268* 1269 LDWRKL = M 1270 CHUNK = N 1271 ELSE 1272* 1273* WORK(IL) is M by CHUNK 1274* 1275 LDWRKL = M 1276 CHUNK = ( LWORK-M*M-3*M ) / M 1277 END IF 1278 ITAU = IL + LDWRKL*CHUNK 1279 NWORK = ITAU + M 1280* 1281* Compute A=L*Q 1282* (CWorkspace: need 2*M, prefer M+M*NB) 1283* (RWorkspace: 0) 1284* 1285 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1286 $ LWORK-NWORK+1, IERR ) 1287* 1288* Copy L to WORK(IL), zeroing about above it 1289* 1290 CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 1291 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 1292 $ WORK( IL+LDWRKL ), LDWRKL ) 1293* 1294* Generate Q in A 1295* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 1296* (RWorkspace: 0) 1297* 1298 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 1299 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1300 IE = 1 1301 ITAUQ = ITAU 1302 ITAUP = ITAUQ + M 1303 NWORK = ITAUP + M 1304* 1305* Bidiagonalize L in WORK(IL) 1306* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 1307* (RWorkspace: need M) 1308* 1309 CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), 1310 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 1311 $ LWORK-NWORK+1, IERR ) 1312* 1313* Perform bidiagonal SVD, computing left singular vectors 1314* of bidiagonal matrix in RWORK(IRU) and computing right 1315* singular vectors of bidiagonal matrix in RWORK(IRVT) 1316* (CWorkspace: need 0) 1317* (RWorkspace: need BDSPAC) 1318* 1319 IRU = IE + M 1320 IRVT = IRU + M*M 1321 NRWORK = IRVT + M*M 1322 CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1323 $ M, RWORK( IRVT ), M, DUM, IDUM, 1324 $ RWORK( NRWORK ), IWORK, INFO ) 1325* 1326* Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 1327* Overwrite WORK(IU) by the left singular vectors of L 1328* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 1329* (RWorkspace: 0) 1330* 1331 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 1332 CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 1333 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1334 $ LWORK-NWORK+1, IERR ) 1335* 1336* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 1337* Overwrite WORK(IVT) by the right singular vectors of L 1338* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 1339* (RWorkspace: 0) 1340* 1341 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 1342 $ LDWKVT ) 1343 CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, 1344 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 1345 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1346* 1347* Multiply right singular vectors of L in WORK(IL) by Q 1348* in A, storing result in WORK(IL) and copying to A 1349* (CWorkspace: need 2*M*M, prefer M*M+M*N)) 1350* (RWorkspace: 0) 1351* 1352 DO 40 I = 1, N, CHUNK 1353 BLK = MIN( N-I+1, CHUNK ) 1354 CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M, 1355 $ A( 1, I ), LDA, CZERO, WORK( IL ), 1356 $ LDWRKL ) 1357 CALL CLACPY( 'F', M, BLK, WORK( IL ), LDWRKL, 1358 $ A( 1, I ), LDA ) 1359 40 CONTINUE 1360* 1361 ELSE IF( WNTQS ) THEN 1362* 1363* Path 3t (N much larger than M, JOBZ='S') 1364* M right singular vectors to be computed in VT and 1365* M left singular vectors to be computed in U 1366* 1367 IL = 1 1368* 1369* WORK(IL) is M by M 1370* 1371 LDWRKL = M 1372 ITAU = IL + LDWRKL*M 1373 NWORK = ITAU + M 1374* 1375* Compute A=L*Q 1376* (CWorkspace: need 2*M, prefer M+M*NB) 1377* (RWorkspace: 0) 1378* 1379 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1380 $ LWORK-NWORK+1, IERR ) 1381* 1382* Copy L to WORK(IL), zeroing out above it 1383* 1384 CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 1385 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, 1386 $ WORK( IL+LDWRKL ), LDWRKL ) 1387* 1388* Generate Q in A 1389* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 1390* (RWorkspace: 0) 1391* 1392 CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 1393 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1394 IE = 1 1395 ITAUQ = ITAU 1396 ITAUP = ITAUQ + M 1397 NWORK = ITAUP + M 1398* 1399* Bidiagonalize L in WORK(IL) 1400* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 1401* (RWorkspace: need M) 1402* 1403 CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), 1404 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 1405 $ LWORK-NWORK+1, IERR ) 1406* 1407* Perform bidiagonal SVD, computing left singular vectors 1408* of bidiagonal matrix in RWORK(IRU) and computing right 1409* singular vectors of bidiagonal matrix in RWORK(IRVT) 1410* (CWorkspace: need 0) 1411* (RWorkspace: need BDSPAC) 1412* 1413 IRU = IE + M 1414 IRVT = IRU + M*M 1415 NRWORK = IRVT + M*M 1416 CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1417 $ M, RWORK( IRVT ), M, DUM, IDUM, 1418 $ RWORK( NRWORK ), IWORK, INFO ) 1419* 1420* Copy real matrix RWORK(IRU) to complex matrix U 1421* Overwrite U by left singular vectors of L 1422* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 1423* (RWorkspace: 0) 1424* 1425 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 1426 CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 1427 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1428 $ LWORK-NWORK+1, IERR ) 1429* 1430* Copy real matrix RWORK(IRVT) to complex matrix VT 1431* Overwrite VT by left singular vectors of L 1432* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 1433* (RWorkspace: 0) 1434* 1435 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 1436 CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, 1437 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1438 $ LWORK-NWORK+1, IERR ) 1439* 1440* Copy VT to WORK(IL), multiply right singular vectors of L 1441* in WORK(IL) by Q in A, storing result in VT 1442* (CWorkspace: need M*M) 1443* (RWorkspace: 0) 1444* 1445 CALL CLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) 1446 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL, 1447 $ A, LDA, CZERO, VT, LDVT ) 1448* 1449 ELSE IF( WNTQA ) THEN 1450* 1451* Path 9t (N much larger than M, JOBZ='A') 1452* N right singular vectors to be computed in VT and 1453* M left singular vectors to be computed in U 1454* 1455 IVT = 1 1456* 1457* WORK(IVT) is M by M 1458* 1459 LDWKVT = M 1460 ITAU = IVT + LDWKVT*M 1461 NWORK = ITAU + M 1462* 1463* Compute A=L*Q, copying result to VT 1464* (CWorkspace: need 2*M, prefer M+M*NB) 1465* (RWorkspace: 0) 1466* 1467 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1468 $ LWORK-NWORK+1, IERR ) 1469 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 1470* 1471* Generate Q in VT 1472* (CWorkspace: need M+N, prefer M+N*NB) 1473* (RWorkspace: 0) 1474* 1475 CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 1476 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1477* 1478* Produce L in A, zeroing out above it 1479* 1480 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 1481 $ LDA ) 1482 IE = 1 1483 ITAUQ = ITAU 1484 ITAUP = ITAUQ + M 1485 NWORK = ITAUP + M 1486* 1487* Bidiagonalize L in A 1488* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 1489* (RWorkspace: need M) 1490* 1491 CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 1492 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1493 $ IERR ) 1494* 1495* Perform bidiagonal SVD, computing left singular vectors 1496* of bidiagonal matrix in RWORK(IRU) and computing right 1497* singular vectors of bidiagonal matrix in RWORK(IRVT) 1498* (CWorkspace: need 0) 1499* (RWorkspace: need BDSPAC) 1500* 1501 IRU = IE + M 1502 IRVT = IRU + M*M 1503 NRWORK = IRVT + M*M 1504 CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1505 $ M, RWORK( IRVT ), M, DUM, IDUM, 1506 $ RWORK( NRWORK ), IWORK, INFO ) 1507* 1508* Copy real matrix RWORK(IRU) to complex matrix U 1509* Overwrite U by left singular vectors of L 1510* (CWorkspace: need 3*M, prefer 2*M+M*NB) 1511* (RWorkspace: 0) 1512* 1513 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 1514 CALL CUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA, 1515 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1516 $ LWORK-NWORK+1, IERR ) 1517* 1518* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 1519* Overwrite WORK(IVT) by right singular vectors of L 1520* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 1521* (RWorkspace: 0) 1522* 1523 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 1524 $ LDWKVT ) 1525 CALL CUNMBR( 'P', 'R', 'C', M, M, M, A, LDA, 1526 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 1527 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1528* 1529* Multiply right singular vectors of L in WORK(IVT) by 1530* Q in VT, storing result in A 1531* (CWorkspace: need M*M) 1532* (RWorkspace: 0) 1533* 1534 CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT, 1535 $ VT, LDVT, CZERO, A, LDA ) 1536* 1537* Copy right singular vectors of A from A to VT 1538* 1539 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 1540* 1541 END IF 1542* 1543 ELSE IF( N.GE.MNTHR2 ) THEN 1544* 1545* MNTHR2 <= N < MNTHR1 1546* 1547* Path 5t (N much larger than M, but not as much as MNTHR1) 1548* Reduce to bidiagonal form without QR decomposition, use 1549* CUNGBR and matrix multiplication to compute singular vectors 1550* 1551* 1552 IE = 1 1553 NRWORK = IE + M 1554 ITAUQ = 1 1555 ITAUP = ITAUQ + M 1556 NWORK = ITAUP + M 1557* 1558* Bidiagonalize A 1559* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 1560* (RWorkspace: M) 1561* 1562 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 1563 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1564 $ IERR ) 1565* 1566 IF( WNTQN ) THEN 1567* 1568* Compute singular values only 1569* (Cworkspace: 0) 1570* (Rworkspace: need BDSPAC) 1571* 1572 CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 1573 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 1574 ELSE IF( WNTQO ) THEN 1575 IRVT = NRWORK 1576 IRU = IRVT + M*M 1577 NRWORK = IRU + M*M 1578 IVT = NWORK 1579* 1580* Copy A to U, generate Q 1581* (Cworkspace: need 2*M, prefer M+M*NB) 1582* (Rworkspace: 0) 1583* 1584 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 1585 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 1586 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1587* 1588* Generate P**H in A 1589* (Cworkspace: need 2*M, prefer M+M*NB) 1590* (Rworkspace: 0) 1591* 1592 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 1593 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1594* 1595 LDWKVT = M 1596 IF( LWORK.GE.M*N+3*M ) THEN 1597* 1598* WORK( IVT ) is M by N 1599* 1600 NWORK = IVT + LDWKVT*N 1601 CHUNK = N 1602 ELSE 1603* 1604* WORK( IVT ) is M by CHUNK 1605* 1606 CHUNK = ( LWORK-3*M ) / M 1607 NWORK = IVT + LDWKVT*CHUNK 1608 END IF 1609* 1610* Perform bidiagonal SVD, computing left singular vectors 1611* of bidiagonal matrix in RWORK(IRU) and computing right 1612* singular vectors of bidiagonal matrix in RWORK(IRVT) 1613* (CWorkspace: need 0) 1614* (RWorkspace: need BDSPAC) 1615* 1616 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1617 $ M, RWORK( IRVT ), M, DUM, IDUM, 1618 $ RWORK( NRWORK ), IWORK, INFO ) 1619* 1620* Multiply Q in U by real matrix RWORK(IRVT) 1621* storing the result in WORK(IVT), copying to U 1622* (Cworkspace: need 0) 1623* (Rworkspace: need 2*M*M) 1624* 1625 CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ), 1626 $ LDWKVT, RWORK( NRWORK ) ) 1627 CALL CLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU ) 1628* 1629* Multiply RWORK(IRVT) by P**H in A, storing the 1630* result in WORK(IVT), copying to A 1631* (CWorkspace: need M*M, prefer M*N) 1632* (Rworkspace: need 2*M*M, prefer 2*M*N) 1633* 1634 NRWORK = IRU 1635 DO 50 I = 1, N, CHUNK 1636 BLK = MIN( N-I+1, CHUNK ) 1637 CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA, 1638 $ WORK( IVT ), LDWKVT, RWORK( NRWORK ) ) 1639 CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT, 1640 $ A( 1, I ), LDA ) 1641 50 CONTINUE 1642 ELSE IF( WNTQS ) THEN 1643* 1644* Copy A to U, generate Q 1645* (Cworkspace: need 2*M, prefer M+M*NB) 1646* (Rworkspace: 0) 1647* 1648 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 1649 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 1650 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1651* 1652* Copy A to VT, generate P**H 1653* (Cworkspace: need 2*M, prefer M+M*NB) 1654* (Rworkspace: 0) 1655* 1656 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 1657 CALL CUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ), 1658 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1659* 1660* Perform bidiagonal SVD, computing left singular vectors 1661* of bidiagonal matrix in RWORK(IRU) and computing right 1662* singular vectors of bidiagonal matrix in RWORK(IRVT) 1663* (CWorkspace: need 0) 1664* (RWorkspace: need BDSPAC) 1665* 1666 IRVT = NRWORK 1667 IRU = IRVT + M*M 1668 NRWORK = IRU + M*M 1669 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1670 $ M, RWORK( IRVT ), M, DUM, IDUM, 1671 $ RWORK( NRWORK ), IWORK, INFO ) 1672* 1673* Multiply Q in U by real matrix RWORK(IRU), storing the 1674* result in A, copying to U 1675* (CWorkspace: need 0) 1676* (Rworkspace: need 3*M*M) 1677* 1678 CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, 1679 $ RWORK( NRWORK ) ) 1680 CALL CLACPY( 'F', M, M, A, LDA, U, LDU ) 1681* 1682* Multiply real matrix RWORK(IRVT) by P**H in VT, 1683* storing the result in A, copying to VT 1684* (Cworkspace: need 0) 1685* (Rworkspace: need M*M+2*M*N) 1686* 1687 NRWORK = IRU 1688 CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, 1689 $ RWORK( NRWORK ) ) 1690 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 1691 ELSE 1692* 1693* Copy A to U, generate Q 1694* (Cworkspace: need 2*M, prefer M+M*NB) 1695* (Rworkspace: 0) 1696* 1697 CALL CLACPY( 'L', M, M, A, LDA, U, LDU ) 1698 CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 1699 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1700* 1701* Copy A to VT, generate P**H 1702* (Cworkspace: need 2*M, prefer M+M*NB) 1703* (Rworkspace: 0) 1704* 1705 CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT ) 1706 CALL CUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ), 1707 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1708* 1709* Perform bidiagonal SVD, computing left singular vectors 1710* of bidiagonal matrix in RWORK(IRU) and computing right 1711* singular vectors of bidiagonal matrix in RWORK(IRVT) 1712* (CWorkspace: need 0) 1713* (RWorkspace: need BDSPAC) 1714* 1715 IRVT = NRWORK 1716 IRU = IRVT + M*M 1717 NRWORK = IRU + M*M 1718 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1719 $ M, RWORK( IRVT ), M, DUM, IDUM, 1720 $ RWORK( NRWORK ), IWORK, INFO ) 1721* 1722* Multiply Q in U by real matrix RWORK(IRU), storing the 1723* result in A, copying to U 1724* (CWorkspace: need 0) 1725* (Rworkspace: need 3*M*M) 1726* 1727 CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, 1728 $ RWORK( NRWORK ) ) 1729 CALL CLACPY( 'F', M, M, A, LDA, U, LDU ) 1730* 1731* Multiply real matrix RWORK(IRVT) by P**H in VT, 1732* storing the result in A, copying to VT 1733* (Cworkspace: need 0) 1734* (Rworkspace: need M*M+2*M*N) 1735* 1736 CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, 1737 $ RWORK( NRWORK ) ) 1738 CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT ) 1739 END IF 1740* 1741 ELSE 1742* 1743* N .LT. MNTHR2 1744* 1745* Path 6t (N greater than M, but not much larger) 1746* Reduce to bidiagonal form without LQ decomposition 1747* Use CUNMBR to compute singular vectors 1748* 1749 IE = 1 1750 NRWORK = IE + M 1751 ITAUQ = 1 1752 ITAUP = ITAUQ + M 1753 NWORK = ITAUP + M 1754* 1755* Bidiagonalize A 1756* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 1757* (RWorkspace: M) 1758* 1759 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 1760 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1761 $ IERR ) 1762 IF( WNTQN ) THEN 1763* 1764* Compute singular values only 1765* (Cworkspace: 0) 1766* (Rworkspace: need BDSPAC) 1767* 1768 CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 1769 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 1770 ELSE IF( WNTQO ) THEN 1771 LDWKVT = M 1772 IVT = NWORK 1773 IF( LWORK.GE.M*N+3*M ) THEN 1774* 1775* WORK( IVT ) is M by N 1776* 1777 CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ), 1778 $ LDWKVT ) 1779 NWORK = IVT + LDWKVT*N 1780 ELSE 1781* 1782* WORK( IVT ) is M by CHUNK 1783* 1784 CHUNK = ( LWORK-3*M ) / M 1785 NWORK = IVT + LDWKVT*CHUNK 1786 END IF 1787* 1788* Perform bidiagonal SVD, computing left singular vectors 1789* of bidiagonal matrix in RWORK(IRU) and computing right 1790* singular vectors of bidiagonal matrix in RWORK(IRVT) 1791* (CWorkspace: need 0) 1792* (RWorkspace: need BDSPAC) 1793* 1794 IRVT = NRWORK 1795 IRU = IRVT + M*M 1796 NRWORK = IRU + M*M 1797 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1798 $ M, RWORK( IRVT ), M, DUM, IDUM, 1799 $ RWORK( NRWORK ), IWORK, INFO ) 1800* 1801* Copy real matrix RWORK(IRU) to complex matrix U 1802* Overwrite U by left singular vectors of A 1803* (Cworkspace: need 2*M, prefer M+M*NB) 1804* (Rworkspace: need 0) 1805* 1806 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 1807 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1808 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1809 $ LWORK-NWORK+1, IERR ) 1810* 1811 IF( LWORK.GE.M*N+3*M ) THEN 1812* 1813* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 1814* Overwrite WORK(IVT) by right singular vectors of A, 1815* copying to A 1816* (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB) 1817* (Rworkspace: need 0) 1818* 1819 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 1820 $ LDWKVT ) 1821 CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA, 1822 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 1823 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1824 CALL CLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) 1825 ELSE 1826* 1827* Generate P**H in A 1828* (Cworkspace: need 2*M, prefer M+M*NB) 1829* (Rworkspace: need 0) 1830* 1831 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 1832 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 1833* 1834* Multiply Q in A by real matrix RWORK(IRU), storing the 1835* result in WORK(IU), copying to A 1836* (CWorkspace: need M*M, prefer M*N) 1837* (Rworkspace: need 3*M*M, prefer M*M+2*M*N) 1838* 1839 NRWORK = IRU 1840 DO 60 I = 1, N, CHUNK 1841 BLK = MIN( N-I+1, CHUNK ) 1842 CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), 1843 $ LDA, WORK( IVT ), LDWKVT, 1844 $ RWORK( NRWORK ) ) 1845 CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT, 1846 $ A( 1, I ), LDA ) 1847 60 CONTINUE 1848 END IF 1849 ELSE IF( WNTQS ) THEN 1850* 1851* Perform bidiagonal SVD, computing left singular vectors 1852* of bidiagonal matrix in RWORK(IRU) and computing right 1853* singular vectors of bidiagonal matrix in RWORK(IRVT) 1854* (CWorkspace: need 0) 1855* (RWorkspace: need BDSPAC) 1856* 1857 IRVT = NRWORK 1858 IRU = IRVT + M*M 1859 NRWORK = IRU + M*M 1860 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1861 $ M, RWORK( IRVT ), M, DUM, IDUM, 1862 $ RWORK( NRWORK ), IWORK, INFO ) 1863* 1864* Copy real matrix RWORK(IRU) to complex matrix U 1865* Overwrite U by left singular vectors of A 1866* (CWorkspace: need 3*M, prefer 2*M+M*NB) 1867* (RWorkspace: M*M) 1868* 1869 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 1870 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1871 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1872 $ LWORK-NWORK+1, IERR ) 1873* 1874* Copy real matrix RWORK(IRVT) to complex matrix VT 1875* Overwrite VT by right singular vectors of A 1876* (CWorkspace: need 3*M, prefer 2*M+M*NB) 1877* (RWorkspace: M*M) 1878* 1879 CALL CLASET( 'F', M, N, CZERO, CZERO, VT, LDVT ) 1880 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 1881 CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA, 1882 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1883 $ LWORK-NWORK+1, IERR ) 1884 ELSE 1885* 1886* Perform bidiagonal SVD, computing left singular vectors 1887* of bidiagonal matrix in RWORK(IRU) and computing right 1888* singular vectors of bidiagonal matrix in RWORK(IRVT) 1889* (CWorkspace: need 0) 1890* (RWorkspace: need BDSPAC) 1891* 1892 IRVT = NRWORK 1893 IRU = IRVT + M*M 1894 NRWORK = IRU + M*M 1895* 1896 CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 1897 $ M, RWORK( IRVT ), M, DUM, IDUM, 1898 $ RWORK( NRWORK ), IWORK, INFO ) 1899* 1900* Copy real matrix RWORK(IRU) to complex matrix U 1901* Overwrite U by left singular vectors of A 1902* (CWorkspace: need 3*M, prefer 2*M+M*NB) 1903* (RWorkspace: M*M) 1904* 1905 CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 1906 CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1907 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1908 $ LWORK-NWORK+1, IERR ) 1909* 1910* Set the right corner of VT to identity matrix 1911* 1912 CALL CLASET( 'F', N-M, N-M, CZERO, CONE, VT( M+1, M+1 ), 1913 $ LDVT ) 1914* 1915* Copy real matrix RWORK(IRVT) to complex matrix VT 1916* Overwrite VT by right singular vectors of A 1917* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 1918* (RWorkspace: M*M) 1919* 1920 CALL CLASET( 'F', N, N, CZERO, CZERO, VT, LDVT ) 1921 CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 1922 CALL CUNMBR( 'P', 'R', 'C', N, N, M, A, LDA, 1923 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1924 $ LWORK-NWORK+1, IERR ) 1925 END IF 1926* 1927 END IF 1928* 1929 END IF 1930* 1931* Undo scaling if necessary 1932* 1933 IF( ISCL.EQ.1 ) THEN 1934 IF( ANRM.GT.BIGNUM ) 1935 $ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 1936 $ IERR ) 1937 IF( ANRM.LT.SMLNUM ) 1938 $ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 1939 $ IERR ) 1940 END IF 1941* 1942* Return optimal workspace in WORK(1) 1943* 1944 WORK( 1 ) = MAXWRK 1945* 1946 RETURN 1947* 1948* End of CGESDD 1949* 1950 END 1951