1*> \brief <b> ZGESVD computes the singular value decomposition (SVD) for GE matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZGESVD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, 22* WORK, LWORK, RWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBU, JOBVT 26* INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION RWORK( * ), S( * ) 30* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 31* $ WORK( * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> ZGESVD computes the singular value decomposition (SVD) of a complex 41*> M-by-N matrix A, optionally computing the left and/or right singular 42*> vectors. The SVD is written 43*> 44*> A = U * SIGMA * conjugate-transpose(V) 45*> 46*> where SIGMA is an M-by-N matrix which is zero except for its 47*> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and 48*> V is an N-by-N unitary matrix. The diagonal elements of SIGMA 49*> are the singular values of A; they are real and non-negative, and 50*> are returned in descending order. The first min(m,n) columns of 51*> U and V are the left and right singular vectors of A. 52*> 53*> Note that the routine returns V**H, not V. 54*> \endverbatim 55* 56* Arguments: 57* ========== 58* 59*> \param[in] JOBU 60*> \verbatim 61*> JOBU is CHARACTER*1 62*> Specifies options for computing all or part of the matrix U: 63*> = 'A': all M columns of U are returned in array U: 64*> = 'S': the first min(m,n) columns of U (the left singular 65*> vectors) are returned in the array U; 66*> = 'O': the first min(m,n) columns of U (the left singular 67*> vectors) are overwritten on the array A; 68*> = 'N': no columns of U (no left singular vectors) are 69*> computed. 70*> \endverbatim 71*> 72*> \param[in] JOBVT 73*> \verbatim 74*> JOBVT is CHARACTER*1 75*> Specifies options for computing all or part of the matrix 76*> V**H: 77*> = 'A': all N rows of V**H are returned in the array VT; 78*> = 'S': the first min(m,n) rows of V**H (the right singular 79*> vectors) are returned in the array VT; 80*> = 'O': the first min(m,n) rows of V**H (the right singular 81*> vectors) are overwritten on the array A; 82*> = 'N': no rows of V**H (no right singular vectors) are 83*> computed. 84*> 85*> JOBVT and JOBU cannot both be 'O'. 86*> \endverbatim 87*> 88*> \param[in] M 89*> \verbatim 90*> M is INTEGER 91*> The number of rows of the input matrix A. M >= 0. 92*> \endverbatim 93*> 94*> \param[in] N 95*> \verbatim 96*> N is INTEGER 97*> The number of columns of the input matrix A. N >= 0. 98*> \endverbatim 99*> 100*> \param[in,out] A 101*> \verbatim 102*> A is COMPLEX*16 array, dimension (LDA,N) 103*> On entry, the M-by-N matrix A. 104*> On exit, 105*> if JOBU = 'O', A is overwritten with the first min(m,n) 106*> columns of U (the left singular vectors, 107*> stored columnwise); 108*> if JOBVT = 'O', A is overwritten with the first min(m,n) 109*> rows of V**H (the right singular vectors, 110*> stored rowwise); 111*> if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A 112*> are destroyed. 113*> \endverbatim 114*> 115*> \param[in] LDA 116*> \verbatim 117*> LDA is INTEGER 118*> The leading dimension of the array A. LDA >= max(1,M). 119*> \endverbatim 120*> 121*> \param[out] S 122*> \verbatim 123*> S is DOUBLE PRECISION array, dimension (min(M,N)) 124*> The singular values of A, sorted so that S(i) >= S(i+1). 125*> \endverbatim 126*> 127*> \param[out] U 128*> \verbatim 129*> U is COMPLEX*16 array, dimension (LDU,UCOL) 130*> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. 131*> If JOBU = 'A', U contains the M-by-M unitary matrix U; 132*> if JOBU = 'S', U contains the first min(m,n) columns of U 133*> (the left singular vectors, stored columnwise); 134*> if JOBU = 'N' or 'O', U is not referenced. 135*> \endverbatim 136*> 137*> \param[in] LDU 138*> \verbatim 139*> LDU is INTEGER 140*> The leading dimension of the array U. LDU >= 1; if 141*> JOBU = 'S' or 'A', LDU >= M. 142*> \endverbatim 143*> 144*> \param[out] VT 145*> \verbatim 146*> VT is COMPLEX*16 array, dimension (LDVT,N) 147*> If JOBVT = 'A', VT contains the N-by-N unitary matrix 148*> V**H; 149*> if JOBVT = 'S', VT contains the first min(m,n) rows of 150*> V**H (the right singular vectors, stored rowwise); 151*> if JOBVT = 'N' or 'O', VT is not referenced. 152*> \endverbatim 153*> 154*> \param[in] LDVT 155*> \verbatim 156*> LDVT is INTEGER 157*> The leading dimension of the array VT. LDVT >= 1; if 158*> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). 159*> \endverbatim 160*> 161*> \param[out] WORK 162*> \verbatim 163*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 164*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 165*> \endverbatim 166*> 167*> \param[in] LWORK 168*> \verbatim 169*> LWORK is INTEGER 170*> The dimension of the array WORK. 171*> LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)). 172*> For good performance, LWORK should generally be larger. 173*> 174*> If LWORK = -1, then a workspace query is assumed; the routine 175*> only calculates the optimal size of the WORK array, returns 176*> this value as the first entry of the WORK array, and no error 177*> message related to LWORK is issued by XERBLA. 178*> \endverbatim 179*> 180*> \param[out] RWORK 181*> \verbatim 182*> RWORK is DOUBLE PRECISION array, dimension (5*min(M,N)) 183*> On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the 184*> unconverged superdiagonal elements of an upper bidiagonal 185*> matrix B whose diagonal is in S (not necessarily sorted). 186*> B satisfies A = U * B * VT, so it has the same singular 187*> values as A, and singular vectors related by U and VT. 188*> \endverbatim 189*> 190*> \param[out] INFO 191*> \verbatim 192*> INFO is INTEGER 193*> = 0: successful exit. 194*> < 0: if INFO = -i, the i-th argument had an illegal value. 195*> > 0: if ZBDSQR did not converge, INFO specifies how many 196*> superdiagonals of an intermediate bidiagonal form B 197*> did not converge to zero. See the description of RWORK 198*> above for details. 199*> \endverbatim 200* 201* Authors: 202* ======== 203* 204*> \author Univ. of Tennessee 205*> \author Univ. of California Berkeley 206*> \author Univ. of Colorado Denver 207*> \author NAG Ltd. 208* 209*> \date April 2012 210* 211*> \ingroup complex16GEsing 212* 213* ===================================================================== 214 SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, 215 $ VT, LDVT, WORK, LWORK, RWORK, INFO ) 216* 217* -- LAPACK driver routine (version 3.7.0) -- 218* -- LAPACK is a software package provided by Univ. of Tennessee, -- 219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 220* April 2012 221* 222* .. Scalar Arguments .. 223 CHARACTER JOBU, JOBVT 224 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 225* .. 226* .. Array Arguments .. 227 DOUBLE PRECISION RWORK( * ), S( * ) 228 COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 229 $ WORK( * ) 230* .. 231* 232* ===================================================================== 233* 234* .. Parameters .. 235 COMPLEX*16 CZERO, CONE 236 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), 237 $ CONE = ( 1.0D0, 0.0D0 ) ) 238 DOUBLE PRECISION ZERO, ONE 239 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 240* .. 241* .. Local Scalars .. 242 LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS, 243 $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS 244 INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL, 245 $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, 246 $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, 247 $ NRVT, WRKBL 248 INTEGER LWORK_ZGEQRF, LWORK_ZUNGQR_N, LWORK_ZUNGQR_M, 249 $ LWORK_ZGEBRD, LWORK_ZUNGBR_P, LWORK_ZUNGBR_Q, 250 $ LWORK_ZGELQF, LWORK_ZUNGLQ_N, LWORK_ZUNGLQ_M 251 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM 252* .. 253* .. Local Arrays .. 254 DOUBLE PRECISION DUM( 1 ) 255 COMPLEX*16 CDUM( 1 ) 256* .. 257* .. External Subroutines .. 258 EXTERNAL DLASCL, XERBLA, ZBDSQR, ZGEBRD, ZGELQF, ZGEMM, 259 $ ZGEQRF, ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNGLQ, 260 $ ZUNGQR, ZUNMBR 261* .. 262* .. External Functions .. 263 LOGICAL LSAME 264 INTEGER ILAENV 265 DOUBLE PRECISION DLAMCH, ZLANGE 266 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE 267* .. 268* .. Intrinsic Functions .. 269 INTRINSIC MAX, MIN, SQRT 270* .. 271* .. Executable Statements .. 272* 273* Test the input arguments 274* 275 INFO = 0 276 MINMN = MIN( M, N ) 277 WNTUA = LSAME( JOBU, 'A' ) 278 WNTUS = LSAME( JOBU, 'S' ) 279 WNTUAS = WNTUA .OR. WNTUS 280 WNTUO = LSAME( JOBU, 'O' ) 281 WNTUN = LSAME( JOBU, 'N' ) 282 WNTVA = LSAME( JOBVT, 'A' ) 283 WNTVS = LSAME( JOBVT, 'S' ) 284 WNTVAS = WNTVA .OR. WNTVS 285 WNTVO = LSAME( JOBVT, 'O' ) 286 WNTVN = LSAME( JOBVT, 'N' ) 287 LQUERY = ( LWORK.EQ.-1 ) 288* 289 IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN 290 INFO = -1 291 ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR. 292 $ ( WNTVO .AND. WNTUO ) ) THEN 293 INFO = -2 294 ELSE IF( M.LT.0 ) THEN 295 INFO = -3 296 ELSE IF( N.LT.0 ) THEN 297 INFO = -4 298 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 299 INFO = -6 300 ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN 301 INFO = -9 302 ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR. 303 $ ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN 304 INFO = -11 305 END IF 306* 307* Compute workspace 308* (Note: Comments in the code beginning "Workspace:" describe the 309* minimal amount of workspace needed at that point in the code, 310* as well as the preferred amount for good performance. 311* CWorkspace refers to complex workspace, and RWorkspace to 312* real workspace. NB refers to the optimal block size for the 313* immediately following subroutine, as returned by ILAENV.) 314* 315 IF( INFO.EQ.0 ) THEN 316 MINWRK = 1 317 MAXWRK = 1 318 IF( M.GE.N .AND. MINMN.GT.0 ) THEN 319* 320* Space needed for ZBDSQR is BDSPAC = 5*N 321* 322 MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) 323* Compute space needed for ZGEQRF 324 CALL ZGEQRF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR ) 325 LWORK_ZGEQRF = INT( CDUM(1) ) 326* Compute space needed for ZUNGQR 327 CALL ZUNGQR( M, N, N, A, LDA, CDUM(1), CDUM(1), -1, IERR ) 328 LWORK_ZUNGQR_N = INT( CDUM(1) ) 329 CALL ZUNGQR( M, M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR ) 330 LWORK_ZUNGQR_M = INT( CDUM(1) ) 331* Compute space needed for ZGEBRD 332 CALL ZGEBRD( N, N, A, LDA, S, DUM(1), CDUM(1), 333 $ CDUM(1), CDUM(1), -1, IERR ) 334 LWORK_ZGEBRD = INT( CDUM(1) ) 335* Compute space needed for ZUNGBR 336 CALL ZUNGBR( 'P', N, N, N, A, LDA, CDUM(1), 337 $ CDUM(1), -1, IERR ) 338 LWORK_ZUNGBR_P = INT( CDUM(1) ) 339 CALL ZUNGBR( 'Q', N, N, N, A, LDA, CDUM(1), 340 $ CDUM(1), -1, IERR ) 341 LWORK_ZUNGBR_Q = INT( CDUM(1) ) 342* 343 IF( M.GE.MNTHR ) THEN 344 IF( WNTUN ) THEN 345* 346* Path 1 (M much larger than N, JOBU='N') 347* 348 MAXWRK = N + LWORK_ZGEQRF 349 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZGEBRD ) 350 IF( WNTVO .OR. WNTVAS ) 351 $ MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P ) 352 MINWRK = 3*N 353 ELSE IF( WNTUO .AND. WNTVN ) THEN 354* 355* Path 2 (M much larger than N, JOBU='O', JOBVT='N') 356* 357 WRKBL = N + LWORK_ZGEQRF 358 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 359 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 360 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 361 MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) 362 MINWRK = 2*N + M 363 ELSE IF( WNTUO .AND. WNTVAS ) THEN 364* 365* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 366* 'A') 367* 368 WRKBL = N + LWORK_ZGEQRF 369 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 370 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 371 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 372 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 373 MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) 374 MINWRK = 2*N + M 375 ELSE IF( WNTUS .AND. WNTVN ) THEN 376* 377* Path 4 (M much larger than N, JOBU='S', JOBVT='N') 378* 379 WRKBL = N + LWORK_ZGEQRF 380 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 381 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 382 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 383 MAXWRK = N*N + WRKBL 384 MINWRK = 2*N + M 385 ELSE IF( WNTUS .AND. WNTVO ) THEN 386* 387* Path 5 (M much larger than N, JOBU='S', JOBVT='O') 388* 389 WRKBL = N + LWORK_ZGEQRF 390 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 391 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 392 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 393 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 394 MAXWRK = 2*N*N + WRKBL 395 MINWRK = 2*N + M 396 ELSE IF( WNTUS .AND. WNTVAS ) THEN 397* 398* Path 6 (M much larger than N, JOBU='S', JOBVT='S' or 399* 'A') 400* 401 WRKBL = N + LWORK_ZGEQRF 402 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 403 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 404 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 405 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 406 MAXWRK = N*N + WRKBL 407 MINWRK = 2*N + M 408 ELSE IF( WNTUA .AND. WNTVN ) THEN 409* 410* Path 7 (M much larger than N, JOBU='A', JOBVT='N') 411* 412 WRKBL = N + LWORK_ZGEQRF 413 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) 414 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 415 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 416 MAXWRK = N*N + WRKBL 417 MINWRK = 2*N + M 418 ELSE IF( WNTUA .AND. WNTVO ) THEN 419* 420* Path 8 (M much larger than N, JOBU='A', JOBVT='O') 421* 422 WRKBL = N + LWORK_ZGEQRF 423 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) 424 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 425 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 426 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 427 MAXWRK = 2*N*N + WRKBL 428 MINWRK = 2*N + M 429 ELSE IF( WNTUA .AND. WNTVAS ) THEN 430* 431* Path 9 (M much larger than N, JOBU='A', JOBVT='S' or 432* 'A') 433* 434 WRKBL = N + LWORK_ZGEQRF 435 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) 436 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 437 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 438 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 439 MAXWRK = N*N + WRKBL 440 MINWRK = 2*N + M 441 END IF 442 ELSE 443* 444* Path 10 (M at least N, but not much larger) 445* 446 CALL ZGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1), 447 $ CDUM(1), CDUM(1), -1, IERR ) 448 LWORK_ZGEBRD = INT( CDUM(1) ) 449 MAXWRK = 2*N + LWORK_ZGEBRD 450 IF( WNTUS .OR. WNTUO ) THEN 451 CALL ZUNGBR( 'Q', M, N, N, A, LDA, CDUM(1), 452 $ CDUM(1), -1, IERR ) 453 LWORK_ZUNGBR_Q = INT( CDUM(1) ) 454 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q ) 455 END IF 456 IF( WNTUA ) THEN 457 CALL ZUNGBR( 'Q', M, M, N, A, LDA, CDUM(1), 458 $ CDUM(1), -1, IERR ) 459 LWORK_ZUNGBR_Q = INT( CDUM(1) ) 460 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q ) 461 END IF 462 IF( .NOT.WNTVN ) THEN 463 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P ) 464 END IF 465 MINWRK = 2*N + M 466 END IF 467 ELSE IF( MINMN.GT.0 ) THEN 468* 469* Space needed for ZBDSQR is BDSPAC = 5*M 470* 471 MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) 472* Compute space needed for ZGELQF 473 CALL ZGELQF( M, N, A, LDA, CDUM(1), CDUM(1), -1, IERR ) 474 LWORK_ZGELQF = INT( CDUM(1) ) 475* Compute space needed for ZUNGLQ 476 CALL ZUNGLQ( N, N, M, CDUM(1), N, CDUM(1), CDUM(1), -1, 477 $ IERR ) 478 LWORK_ZUNGLQ_N = INT( CDUM(1) ) 479 CALL ZUNGLQ( M, N, M, A, LDA, CDUM(1), CDUM(1), -1, IERR ) 480 LWORK_ZUNGLQ_M = INT( CDUM(1) ) 481* Compute space needed for ZGEBRD 482 CALL ZGEBRD( M, M, A, LDA, S, DUM(1), CDUM(1), 483 $ CDUM(1), CDUM(1), -1, IERR ) 484 LWORK_ZGEBRD = INT( CDUM(1) ) 485* Compute space needed for ZUNGBR P 486 CALL ZUNGBR( 'P', M, M, M, A, N, CDUM(1), 487 $ CDUM(1), -1, IERR ) 488 LWORK_ZUNGBR_P = INT( CDUM(1) ) 489* Compute space needed for ZUNGBR Q 490 CALL ZUNGBR( 'Q', M, M, M, A, N, CDUM(1), 491 $ CDUM(1), -1, IERR ) 492 LWORK_ZUNGBR_Q = INT( CDUM(1) ) 493 IF( N.GE.MNTHR ) THEN 494 IF( WNTVN ) THEN 495* 496* Path 1t(N much larger than M, JOBVT='N') 497* 498 MAXWRK = M + LWORK_ZGELQF 499 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZGEBRD ) 500 IF( WNTUO .OR. WNTUAS ) 501 $ MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q ) 502 MINWRK = 3*M 503 ELSE IF( WNTVO .AND. WNTUN ) THEN 504* 505* Path 2t(N much larger than M, JOBU='N', JOBVT='O') 506* 507 WRKBL = M + LWORK_ZGELQF 508 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 509 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 510 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 511 MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) 512 MINWRK = 2*M + N 513 ELSE IF( WNTVO .AND. WNTUAS ) THEN 514* 515* Path 3t(N much larger than M, JOBU='S' or 'A', 516* JOBVT='O') 517* 518 WRKBL = M + LWORK_ZGELQF 519 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 520 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 521 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 522 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 523 MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) 524 MINWRK = 2*M + N 525 ELSE IF( WNTVS .AND. WNTUN ) THEN 526* 527* Path 4t(N much larger than M, JOBU='N', JOBVT='S') 528* 529 WRKBL = M + LWORK_ZGELQF 530 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 531 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 532 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 533 MAXWRK = M*M + WRKBL 534 MINWRK = 2*M + N 535 ELSE IF( WNTVS .AND. WNTUO ) THEN 536* 537* Path 5t(N much larger than M, JOBU='O', JOBVT='S') 538* 539 WRKBL = M + LWORK_ZGELQF 540 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 541 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 542 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 543 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 544 MAXWRK = 2*M*M + WRKBL 545 MINWRK = 2*M + N 546 ELSE IF( WNTVS .AND. WNTUAS ) THEN 547* 548* Path 6t(N much larger than M, JOBU='S' or 'A', 549* JOBVT='S') 550* 551 WRKBL = M + LWORK_ZGELQF 552 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 553 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 554 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 555 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 556 MAXWRK = M*M + WRKBL 557 MINWRK = 2*M + N 558 ELSE IF( WNTVA .AND. WNTUN ) THEN 559* 560* Path 7t(N much larger than M, JOBU='N', JOBVT='A') 561* 562 WRKBL = M + LWORK_ZGELQF 563 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) 564 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 565 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 566 MAXWRK = M*M + WRKBL 567 MINWRK = 2*M + N 568 ELSE IF( WNTVA .AND. WNTUO ) THEN 569* 570* Path 8t(N much larger than M, JOBU='O', JOBVT='A') 571* 572 WRKBL = M + LWORK_ZGELQF 573 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) 574 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 575 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 576 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 577 MAXWRK = 2*M*M + WRKBL 578 MINWRK = 2*M + N 579 ELSE IF( WNTVA .AND. WNTUAS ) THEN 580* 581* Path 9t(N much larger than M, JOBU='S' or 'A', 582* JOBVT='A') 583* 584 WRKBL = M + LWORK_ZGELQF 585 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) 586 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 587 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 588 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 589 MAXWRK = M*M + WRKBL 590 MINWRK = 2*M + N 591 END IF 592 ELSE 593* 594* Path 10t(N greater than M, but not much larger) 595* 596 CALL ZGEBRD( M, N, A, LDA, S, DUM(1), CDUM(1), 597 $ CDUM(1), CDUM(1), -1, IERR ) 598 LWORK_ZGEBRD = INT( CDUM(1) ) 599 MAXWRK = 2*M + LWORK_ZGEBRD 600 IF( WNTVS .OR. WNTVO ) THEN 601* Compute space needed for ZUNGBR P 602 CALL ZUNGBR( 'P', M, N, M, A, N, CDUM(1), 603 $ CDUM(1), -1, IERR ) 604 LWORK_ZUNGBR_P = INT( CDUM(1) ) 605 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P ) 606 END IF 607 IF( WNTVA ) THEN 608 CALL ZUNGBR( 'P', N, N, M, A, N, CDUM(1), 609 $ CDUM(1), -1, IERR ) 610 LWORK_ZUNGBR_P = INT( CDUM(1) ) 611 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P ) 612 END IF 613 IF( .NOT.WNTUN ) THEN 614 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q ) 615 END IF 616 MINWRK = 2*M + N 617 END IF 618 END IF 619 MAXWRK = MAX( MAXWRK, MINWRK ) 620 WORK( 1 ) = MAXWRK 621* 622 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 623 INFO = -13 624 END IF 625 END IF 626* 627 IF( INFO.NE.0 ) THEN 628 CALL XERBLA( 'ZGESVD', -INFO ) 629 RETURN 630 ELSE IF( LQUERY ) THEN 631 RETURN 632 END IF 633* 634* Quick return if possible 635* 636 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 637 RETURN 638 END IF 639* 640* Get machine constants 641* 642 EPS = DLAMCH( 'P' ) 643 SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS 644 BIGNUM = ONE / SMLNUM 645* 646* Scale A if max element outside range [SMLNUM,BIGNUM] 647* 648 ANRM = ZLANGE( 'M', M, N, A, LDA, DUM ) 649 ISCL = 0 650 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 651 ISCL = 1 652 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 653 ELSE IF( ANRM.GT.BIGNUM ) THEN 654 ISCL = 1 655 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 656 END IF 657* 658 IF( M.GE.N ) THEN 659* 660* A has at least as many rows as columns. If A has sufficiently 661* more rows than columns, first reduce using the QR 662* decomposition (if sufficient workspace available) 663* 664 IF( M.GE.MNTHR ) THEN 665* 666 IF( WNTUN ) THEN 667* 668* Path 1 (M much larger than N, JOBU='N') 669* No left singular vectors to be computed 670* 671 ITAU = 1 672 IWORK = ITAU + N 673* 674* Compute A=Q*R 675* (CWorkspace: need 2*N, prefer N+N*NB) 676* (RWorkspace: need 0) 677* 678 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 679 $ LWORK-IWORK+1, IERR ) 680* 681* Zero out below R 682* 683 IF( N .GT. 1 ) THEN 684 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 685 $ LDA ) 686 END IF 687 IE = 1 688 ITAUQ = 1 689 ITAUP = ITAUQ + N 690 IWORK = ITAUP + N 691* 692* Bidiagonalize R in A 693* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 694* (RWorkspace: need N) 695* 696 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 697 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 698 $ IERR ) 699 NCVT = 0 700 IF( WNTVO .OR. WNTVAS ) THEN 701* 702* If right singular vectors desired, generate P'. 703* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 704* (RWorkspace: 0) 705* 706 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 707 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 708 NCVT = N 709 END IF 710 IRWORK = IE + N 711* 712* Perform bidiagonal QR iteration, computing right 713* singular vectors of A in A if desired 714* (CWorkspace: 0) 715* (RWorkspace: need BDSPAC) 716* 717 CALL ZBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA, 718 $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO ) 719* 720* If right singular vectors desired in VT, copy them there 721* 722 IF( WNTVAS ) 723 $ CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT ) 724* 725 ELSE IF( WNTUO .AND. WNTVN ) THEN 726* 727* Path 2 (M much larger than N, JOBU='O', JOBVT='N') 728* N left singular vectors to be overwritten on A and 729* no right singular vectors to be computed 730* 731 IF( LWORK.GE.N*N+3*N ) THEN 732* 733* Sufficient workspace for a fast algorithm 734* 735 IR = 1 736 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN 737* 738* WORK(IU) is LDA by N, WORK(IR) is LDA by N 739* 740 LDWRKU = LDA 741 LDWRKR = LDA 742 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN 743* 744* WORK(IU) is LDA by N, WORK(IR) is N by N 745* 746 LDWRKU = LDA 747 LDWRKR = N 748 ELSE 749* 750* WORK(IU) is LDWRKU by N, WORK(IR) is N by N 751* 752 LDWRKU = ( LWORK-N*N ) / N 753 LDWRKR = N 754 END IF 755 ITAU = IR + LDWRKR*N 756 IWORK = ITAU + N 757* 758* Compute A=Q*R 759* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 760* (RWorkspace: 0) 761* 762 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 763 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 764* 765* Copy R to WORK(IR) and zero out below it 766* 767 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 768 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 769 $ WORK( IR+1 ), LDWRKR ) 770* 771* Generate Q in A 772* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 773* (RWorkspace: 0) 774* 775 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 776 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 777 IE = 1 778 ITAUQ = ITAU 779 ITAUP = ITAUQ + N 780 IWORK = ITAUP + N 781* 782* Bidiagonalize R in WORK(IR) 783* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 784* (RWorkspace: need N) 785* 786 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 787 $ WORK( ITAUQ ), WORK( ITAUP ), 788 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 789* 790* Generate left vectors bidiagonalizing R 791* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 792* (RWorkspace: need 0) 793* 794 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 795 $ WORK( ITAUQ ), WORK( IWORK ), 796 $ LWORK-IWORK+1, IERR ) 797 IRWORK = IE + N 798* 799* Perform bidiagonal QR iteration, computing left 800* singular vectors of R in WORK(IR) 801* (CWorkspace: need N*N) 802* (RWorkspace: need BDSPAC) 803* 804 CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1, 805 $ WORK( IR ), LDWRKR, CDUM, 1, 806 $ RWORK( IRWORK ), INFO ) 807 IU = ITAUQ 808* 809* Multiply Q in A by left singular vectors of R in 810* WORK(IR), storing result in WORK(IU) and copying to A 811* (CWorkspace: need N*N+N, prefer N*N+M*N) 812* (RWorkspace: 0) 813* 814 DO 10 I = 1, M, LDWRKU 815 CHUNK = MIN( M-I+1, LDWRKU ) 816 CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 817 $ LDA, WORK( IR ), LDWRKR, CZERO, 818 $ WORK( IU ), LDWRKU ) 819 CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 820 $ A( I, 1 ), LDA ) 821 10 CONTINUE 822* 823 ELSE 824* 825* Insufficient workspace for a fast algorithm 826* 827 IE = 1 828 ITAUQ = 1 829 ITAUP = ITAUQ + N 830 IWORK = ITAUP + N 831* 832* Bidiagonalize A 833* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 834* (RWorkspace: N) 835* 836 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), 837 $ WORK( ITAUQ ), WORK( ITAUP ), 838 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 839* 840* Generate left vectors bidiagonalizing A 841* (CWorkspace: need 3*N, prefer 2*N+N*NB) 842* (RWorkspace: 0) 843* 844 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 845 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 846 IRWORK = IE + N 847* 848* Perform bidiagonal QR iteration, computing left 849* singular vectors of A in A 850* (CWorkspace: need 0) 851* (RWorkspace: need BDSPAC) 852* 853 CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1, 854 $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO ) 855* 856 END IF 857* 858 ELSE IF( WNTUO .AND. WNTVAS ) THEN 859* 860* Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A') 861* N left singular vectors to be overwritten on A and 862* N right singular vectors to be computed in VT 863* 864 IF( LWORK.GE.N*N+3*N ) THEN 865* 866* Sufficient workspace for a fast algorithm 867* 868 IR = 1 869 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN 870* 871* WORK(IU) is LDA by N and WORK(IR) is LDA by N 872* 873 LDWRKU = LDA 874 LDWRKR = LDA 875 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN 876* 877* WORK(IU) is LDA by N and WORK(IR) is N by N 878* 879 LDWRKU = LDA 880 LDWRKR = N 881 ELSE 882* 883* WORK(IU) is LDWRKU by N and WORK(IR) is N by N 884* 885 LDWRKU = ( LWORK-N*N ) / N 886 LDWRKR = N 887 END IF 888 ITAU = IR + LDWRKR*N 889 IWORK = ITAU + N 890* 891* Compute A=Q*R 892* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 893* (RWorkspace: 0) 894* 895 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 896 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 897* 898* Copy R to VT, zeroing out below it 899* 900 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 901 IF( N.GT.1 ) 902 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 903 $ VT( 2, 1 ), LDVT ) 904* 905* Generate Q in A 906* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 907* (RWorkspace: 0) 908* 909 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 910 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 911 IE = 1 912 ITAUQ = ITAU 913 ITAUP = ITAUQ + N 914 IWORK = ITAUP + N 915* 916* Bidiagonalize R in VT, copying result to WORK(IR) 917* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 918* (RWorkspace: need N) 919* 920 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 921 $ WORK( ITAUQ ), WORK( ITAUP ), 922 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 923 CALL ZLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) 924* 925* Generate left vectors bidiagonalizing R in WORK(IR) 926* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 927* (RWorkspace: 0) 928* 929 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 930 $ WORK( ITAUQ ), WORK( IWORK ), 931 $ LWORK-IWORK+1, IERR ) 932* 933* Generate right vectors bidiagonalizing R in VT 934* (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB) 935* (RWorkspace: 0) 936* 937 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 938 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 939 IRWORK = IE + N 940* 941* Perform bidiagonal QR iteration, computing left 942* singular vectors of R in WORK(IR) and computing right 943* singular vectors of R in VT 944* (CWorkspace: need N*N) 945* (RWorkspace: need BDSPAC) 946* 947 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 948 $ LDVT, WORK( IR ), LDWRKR, CDUM, 1, 949 $ RWORK( IRWORK ), INFO ) 950 IU = ITAUQ 951* 952* Multiply Q in A by left singular vectors of R in 953* WORK(IR), storing result in WORK(IU) and copying to A 954* (CWorkspace: need N*N+N, prefer N*N+M*N) 955* (RWorkspace: 0) 956* 957 DO 20 I = 1, M, LDWRKU 958 CHUNK = MIN( M-I+1, LDWRKU ) 959 CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 960 $ LDA, WORK( IR ), LDWRKR, CZERO, 961 $ WORK( IU ), LDWRKU ) 962 CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 963 $ A( I, 1 ), LDA ) 964 20 CONTINUE 965* 966 ELSE 967* 968* Insufficient workspace for a fast algorithm 969* 970 ITAU = 1 971 IWORK = ITAU + N 972* 973* Compute A=Q*R 974* (CWorkspace: need 2*N, prefer N+N*NB) 975* (RWorkspace: 0) 976* 977 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 978 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 979* 980* Copy R to VT, zeroing out below it 981* 982 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 983 IF( N.GT.1 ) 984 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 985 $ VT( 2, 1 ), LDVT ) 986* 987* Generate Q in A 988* (CWorkspace: need 2*N, prefer N+N*NB) 989* (RWorkspace: 0) 990* 991 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 992 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 993 IE = 1 994 ITAUQ = ITAU 995 ITAUP = ITAUQ + N 996 IWORK = ITAUP + N 997* 998* Bidiagonalize R in VT 999* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 1000* (RWorkspace: N) 1001* 1002 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 1003 $ WORK( ITAUQ ), WORK( ITAUP ), 1004 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1005* 1006* Multiply Q in A by left vectors bidiagonalizing R 1007* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1008* (RWorkspace: 0) 1009* 1010 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 1011 $ WORK( ITAUQ ), A, LDA, WORK( IWORK ), 1012 $ LWORK-IWORK+1, IERR ) 1013* 1014* Generate right vectors bidiagonalizing R in VT 1015* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 1016* (RWorkspace: 0) 1017* 1018 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 1019 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1020 IRWORK = IE + N 1021* 1022* Perform bidiagonal QR iteration, computing left 1023* singular vectors of A in A and computing right 1024* singular vectors of A in VT 1025* (CWorkspace: 0) 1026* (RWorkspace: need BDSPAC) 1027* 1028 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 1029 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 1030 $ INFO ) 1031* 1032 END IF 1033* 1034 ELSE IF( WNTUS ) THEN 1035* 1036 IF( WNTVN ) THEN 1037* 1038* Path 4 (M much larger than N, JOBU='S', JOBVT='N') 1039* N left singular vectors to be computed in U and 1040* no right singular vectors to be computed 1041* 1042 IF( LWORK.GE.N*N+3*N ) THEN 1043* 1044* Sufficient workspace for a fast algorithm 1045* 1046 IR = 1 1047 IF( LWORK.GE.WRKBL+LDA*N ) THEN 1048* 1049* WORK(IR) is LDA by N 1050* 1051 LDWRKR = LDA 1052 ELSE 1053* 1054* WORK(IR) is N by N 1055* 1056 LDWRKR = N 1057 END IF 1058 ITAU = IR + LDWRKR*N 1059 IWORK = ITAU + N 1060* 1061* Compute A=Q*R 1062* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 1063* (RWorkspace: 0) 1064* 1065 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1066 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1067* 1068* Copy R to WORK(IR), zeroing out below it 1069* 1070 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), 1071 $ LDWRKR ) 1072 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1073 $ WORK( IR+1 ), LDWRKR ) 1074* 1075* Generate Q in A 1076* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 1077* (RWorkspace: 0) 1078* 1079 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 1080 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1081 IE = 1 1082 ITAUQ = ITAU 1083 ITAUP = ITAUQ + N 1084 IWORK = ITAUP + N 1085* 1086* Bidiagonalize R in WORK(IR) 1087* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 1088* (RWorkspace: need N) 1089* 1090 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, 1091 $ RWORK( IE ), WORK( ITAUQ ), 1092 $ WORK( ITAUP ), WORK( IWORK ), 1093 $ LWORK-IWORK+1, IERR ) 1094* 1095* Generate left vectors bidiagonalizing R in WORK(IR) 1096* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 1097* (RWorkspace: 0) 1098* 1099 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 1100 $ WORK( ITAUQ ), WORK( IWORK ), 1101 $ LWORK-IWORK+1, IERR ) 1102 IRWORK = IE + N 1103* 1104* Perform bidiagonal QR iteration, computing left 1105* singular vectors of R in WORK(IR) 1106* (CWorkspace: need N*N) 1107* (RWorkspace: need BDSPAC) 1108* 1109 CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1110 $ 1, WORK( IR ), LDWRKR, CDUM, 1, 1111 $ RWORK( IRWORK ), INFO ) 1112* 1113* Multiply Q in A by left singular vectors of R in 1114* WORK(IR), storing result in U 1115* (CWorkspace: need N*N) 1116* (RWorkspace: 0) 1117* 1118 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 1119 $ WORK( IR ), LDWRKR, CZERO, U, LDU ) 1120* 1121 ELSE 1122* 1123* Insufficient workspace for a fast algorithm 1124* 1125 ITAU = 1 1126 IWORK = ITAU + N 1127* 1128* Compute A=Q*R, copying result to U 1129* (CWorkspace: need 2*N, prefer N+N*NB) 1130* (RWorkspace: 0) 1131* 1132 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1133 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1134 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1135* 1136* Generate Q in U 1137* (CWorkspace: need 2*N, prefer N+N*NB) 1138* (RWorkspace: 0) 1139* 1140 CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ), 1141 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1142 IE = 1 1143 ITAUQ = ITAU 1144 ITAUP = ITAUQ + N 1145 IWORK = ITAUP + N 1146* 1147* Zero out below R in A 1148* 1149 IF( N .GT. 1 ) THEN 1150 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1151 $ A( 2, 1 ), LDA ) 1152 END IF 1153* 1154* Bidiagonalize R in A 1155* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 1156* (RWorkspace: need N) 1157* 1158 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 1159 $ WORK( ITAUQ ), WORK( ITAUP ), 1160 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1161* 1162* Multiply Q in U by left vectors bidiagonalizing R 1163* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1164* (RWorkspace: 0) 1165* 1166 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 1167 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 1168 $ LWORK-IWORK+1, IERR ) 1169 IRWORK = IE + N 1170* 1171* Perform bidiagonal QR iteration, computing left 1172* singular vectors of A in U 1173* (CWorkspace: 0) 1174* (RWorkspace: need BDSPAC) 1175* 1176 CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1177 $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ), 1178 $ INFO ) 1179* 1180 END IF 1181* 1182 ELSE IF( WNTVO ) THEN 1183* 1184* Path 5 (M much larger than N, JOBU='S', JOBVT='O') 1185* N left singular vectors to be computed in U and 1186* N right singular vectors to be overwritten on A 1187* 1188 IF( LWORK.GE.2*N*N+3*N ) THEN 1189* 1190* Sufficient workspace for a fast algorithm 1191* 1192 IU = 1 1193 IF( LWORK.GE.WRKBL+2*LDA*N ) THEN 1194* 1195* WORK(IU) is LDA by N and WORK(IR) is LDA by N 1196* 1197 LDWRKU = LDA 1198 IR = IU + LDWRKU*N 1199 LDWRKR = LDA 1200 ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN 1201* 1202* WORK(IU) is LDA by N and WORK(IR) is N by N 1203* 1204 LDWRKU = LDA 1205 IR = IU + LDWRKU*N 1206 LDWRKR = N 1207 ELSE 1208* 1209* WORK(IU) is N by N and WORK(IR) is N by N 1210* 1211 LDWRKU = N 1212 IR = IU + LDWRKU*N 1213 LDWRKR = N 1214 END IF 1215 ITAU = IR + LDWRKR*N 1216 IWORK = ITAU + N 1217* 1218* Compute A=Q*R 1219* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 1220* (RWorkspace: 0) 1221* 1222 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1223 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1224* 1225* Copy R to WORK(IU), zeroing out below it 1226* 1227 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 1228 $ LDWRKU ) 1229 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1230 $ WORK( IU+1 ), LDWRKU ) 1231* 1232* Generate Q in A 1233* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 1234* (RWorkspace: 0) 1235* 1236 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 1237 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1238 IE = 1 1239 ITAUQ = ITAU 1240 ITAUP = ITAUQ + N 1241 IWORK = ITAUP + N 1242* 1243* Bidiagonalize R in WORK(IU), copying result to 1244* WORK(IR) 1245* (CWorkspace: need 2*N*N+3*N, 1246* prefer 2*N*N+2*N+2*N*NB) 1247* (RWorkspace: need N) 1248* 1249 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 1250 $ RWORK( IE ), WORK( ITAUQ ), 1251 $ WORK( ITAUP ), WORK( IWORK ), 1252 $ LWORK-IWORK+1, IERR ) 1253 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, 1254 $ WORK( IR ), LDWRKR ) 1255* 1256* Generate left bidiagonalizing vectors in WORK(IU) 1257* (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 1258* (RWorkspace: 0) 1259* 1260 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 1261 $ WORK( ITAUQ ), WORK( IWORK ), 1262 $ LWORK-IWORK+1, IERR ) 1263* 1264* Generate right bidiagonalizing vectors in WORK(IR) 1265* (CWorkspace: need 2*N*N+3*N-1, 1266* prefer 2*N*N+2*N+(N-1)*NB) 1267* (RWorkspace: 0) 1268* 1269 CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR, 1270 $ WORK( ITAUP ), WORK( IWORK ), 1271 $ LWORK-IWORK+1, IERR ) 1272 IRWORK = IE + N 1273* 1274* Perform bidiagonal QR iteration, computing left 1275* singular vectors of R in WORK(IU) and computing 1276* right singular vectors of R in WORK(IR) 1277* (CWorkspace: need 2*N*N) 1278* (RWorkspace: need BDSPAC) 1279* 1280 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), 1281 $ WORK( IR ), LDWRKR, WORK( IU ), 1282 $ LDWRKU, CDUM, 1, RWORK( IRWORK ), 1283 $ INFO ) 1284* 1285* Multiply Q in A by left singular vectors of R in 1286* WORK(IU), storing result in U 1287* (CWorkspace: need N*N) 1288* (RWorkspace: 0) 1289* 1290 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 1291 $ WORK( IU ), LDWRKU, CZERO, U, LDU ) 1292* 1293* Copy right singular vectors of R to A 1294* (CWorkspace: need N*N) 1295* (RWorkspace: 0) 1296* 1297 CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A, 1298 $ LDA ) 1299* 1300 ELSE 1301* 1302* Insufficient workspace for a fast algorithm 1303* 1304 ITAU = 1 1305 IWORK = ITAU + N 1306* 1307* Compute A=Q*R, copying result to U 1308* (CWorkspace: need 2*N, prefer N+N*NB) 1309* (RWorkspace: 0) 1310* 1311 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1312 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1313 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1314* 1315* Generate Q in U 1316* (CWorkspace: need 2*N, prefer N+N*NB) 1317* (RWorkspace: 0) 1318* 1319 CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ), 1320 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1321 IE = 1 1322 ITAUQ = ITAU 1323 ITAUP = ITAUQ + N 1324 IWORK = ITAUP + N 1325* 1326* Zero out below R in A 1327* 1328 IF( N .GT. 1 ) THEN 1329 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1330 $ A( 2, 1 ), LDA ) 1331 END IF 1332* 1333* Bidiagonalize R in A 1334* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 1335* (RWorkspace: need N) 1336* 1337 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 1338 $ WORK( ITAUQ ), WORK( ITAUP ), 1339 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1340* 1341* Multiply Q in U by left vectors bidiagonalizing R 1342* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1343* (RWorkspace: 0) 1344* 1345 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 1346 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 1347 $ LWORK-IWORK+1, IERR ) 1348* 1349* Generate right vectors bidiagonalizing R in A 1350* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 1351* (RWorkspace: 0) 1352* 1353 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 1354 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1355 IRWORK = IE + N 1356* 1357* Perform bidiagonal QR iteration, computing left 1358* singular vectors of A in U and computing right 1359* singular vectors of A in A 1360* (CWorkspace: 0) 1361* (RWorkspace: need BDSPAC) 1362* 1363 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A, 1364 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 1365 $ INFO ) 1366* 1367 END IF 1368* 1369 ELSE IF( WNTVAS ) THEN 1370* 1371* Path 6 (M much larger than N, JOBU='S', JOBVT='S' 1372* or 'A') 1373* N left singular vectors to be computed in U and 1374* N right singular vectors to be computed in VT 1375* 1376 IF( LWORK.GE.N*N+3*N ) THEN 1377* 1378* Sufficient workspace for a fast algorithm 1379* 1380 IU = 1 1381 IF( LWORK.GE.WRKBL+LDA*N ) THEN 1382* 1383* WORK(IU) is LDA by N 1384* 1385 LDWRKU = LDA 1386 ELSE 1387* 1388* WORK(IU) is N by N 1389* 1390 LDWRKU = N 1391 END IF 1392 ITAU = IU + LDWRKU*N 1393 IWORK = ITAU + N 1394* 1395* Compute A=Q*R 1396* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 1397* (RWorkspace: 0) 1398* 1399 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1400 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1401* 1402* Copy R to WORK(IU), zeroing out below it 1403* 1404 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 1405 $ LDWRKU ) 1406 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1407 $ WORK( IU+1 ), LDWRKU ) 1408* 1409* Generate Q in A 1410* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 1411* (RWorkspace: 0) 1412* 1413 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 1414 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1415 IE = 1 1416 ITAUQ = ITAU 1417 ITAUP = ITAUQ + N 1418 IWORK = ITAUP + N 1419* 1420* Bidiagonalize R in WORK(IU), copying result to VT 1421* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 1422* (RWorkspace: need N) 1423* 1424 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 1425 $ RWORK( IE ), WORK( ITAUQ ), 1426 $ WORK( ITAUP ), WORK( IWORK ), 1427 $ LWORK-IWORK+1, IERR ) 1428 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT, 1429 $ LDVT ) 1430* 1431* Generate left bidiagonalizing vectors in WORK(IU) 1432* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 1433* (RWorkspace: 0) 1434* 1435 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 1436 $ WORK( ITAUQ ), WORK( IWORK ), 1437 $ LWORK-IWORK+1, IERR ) 1438* 1439* Generate right bidiagonalizing vectors in VT 1440* (CWorkspace: need N*N+3*N-1, 1441* prefer N*N+2*N+(N-1)*NB) 1442* (RWorkspace: 0) 1443* 1444 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 1445 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1446 IRWORK = IE + N 1447* 1448* Perform bidiagonal QR iteration, computing left 1449* singular vectors of R in WORK(IU) and computing 1450* right singular vectors of R in VT 1451* (CWorkspace: need N*N) 1452* (RWorkspace: need BDSPAC) 1453* 1454 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 1455 $ LDVT, WORK( IU ), LDWRKU, CDUM, 1, 1456 $ RWORK( IRWORK ), INFO ) 1457* 1458* Multiply Q in A by left singular vectors of R in 1459* WORK(IU), storing result in U 1460* (CWorkspace: need N*N) 1461* (RWorkspace: 0) 1462* 1463 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 1464 $ WORK( IU ), LDWRKU, CZERO, U, LDU ) 1465* 1466 ELSE 1467* 1468* Insufficient workspace for a fast algorithm 1469* 1470 ITAU = 1 1471 IWORK = ITAU + N 1472* 1473* Compute A=Q*R, copying result to U 1474* (CWorkspace: need 2*N, prefer N+N*NB) 1475* (RWorkspace: 0) 1476* 1477 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1478 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1479 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1480* 1481* Generate Q in U 1482* (CWorkspace: need 2*N, prefer N+N*NB) 1483* (RWorkspace: 0) 1484* 1485 CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ), 1486 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1487* 1488* Copy R to VT, zeroing out below it 1489* 1490 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 1491 IF( N.GT.1 ) 1492 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1493 $ VT( 2, 1 ), LDVT ) 1494 IE = 1 1495 ITAUQ = ITAU 1496 ITAUP = ITAUQ + N 1497 IWORK = ITAUP + N 1498* 1499* Bidiagonalize R in VT 1500* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 1501* (RWorkspace: need N) 1502* 1503 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 1504 $ WORK( ITAUQ ), WORK( ITAUP ), 1505 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1506* 1507* Multiply Q in U by left bidiagonalizing vectors 1508* in VT 1509* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1510* (RWorkspace: 0) 1511* 1512 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 1513 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 1514 $ LWORK-IWORK+1, IERR ) 1515* 1516* Generate right bidiagonalizing vectors in VT 1517* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 1518* (RWorkspace: 0) 1519* 1520 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 1521 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1522 IRWORK = IE + N 1523* 1524* Perform bidiagonal QR iteration, computing left 1525* singular vectors of A in U and computing right 1526* singular vectors of A in VT 1527* (CWorkspace: 0) 1528* (RWorkspace: need BDSPAC) 1529* 1530 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 1531 $ LDVT, U, LDU, CDUM, 1, 1532 $ RWORK( IRWORK ), INFO ) 1533* 1534 END IF 1535* 1536 END IF 1537* 1538 ELSE IF( WNTUA ) THEN 1539* 1540 IF( WNTVN ) THEN 1541* 1542* Path 7 (M much larger than N, JOBU='A', JOBVT='N') 1543* M left singular vectors to be computed in U and 1544* no right singular vectors to be computed 1545* 1546 IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN 1547* 1548* Sufficient workspace for a fast algorithm 1549* 1550 IR = 1 1551 IF( LWORK.GE.WRKBL+LDA*N ) THEN 1552* 1553* WORK(IR) is LDA by N 1554* 1555 LDWRKR = LDA 1556 ELSE 1557* 1558* WORK(IR) is N by N 1559* 1560 LDWRKR = N 1561 END IF 1562 ITAU = IR + LDWRKR*N 1563 IWORK = ITAU + N 1564* 1565* Compute A=Q*R, copying result to U 1566* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 1567* (RWorkspace: 0) 1568* 1569 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1570 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1571 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1572* 1573* Copy R to WORK(IR), zeroing out below it 1574* 1575 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), 1576 $ LDWRKR ) 1577 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1578 $ WORK( IR+1 ), LDWRKR ) 1579* 1580* Generate Q in U 1581* (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) 1582* (RWorkspace: 0) 1583* 1584 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 1585 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1586 IE = 1 1587 ITAUQ = ITAU 1588 ITAUP = ITAUQ + N 1589 IWORK = ITAUP + N 1590* 1591* Bidiagonalize R in WORK(IR) 1592* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 1593* (RWorkspace: need N) 1594* 1595 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, 1596 $ RWORK( IE ), WORK( ITAUQ ), 1597 $ WORK( ITAUP ), WORK( IWORK ), 1598 $ LWORK-IWORK+1, IERR ) 1599* 1600* Generate left bidiagonalizing vectors in WORK(IR) 1601* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 1602* (RWorkspace: 0) 1603* 1604 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 1605 $ WORK( ITAUQ ), WORK( IWORK ), 1606 $ LWORK-IWORK+1, IERR ) 1607 IRWORK = IE + N 1608* 1609* Perform bidiagonal QR iteration, computing left 1610* singular vectors of R in WORK(IR) 1611* (CWorkspace: need N*N) 1612* (RWorkspace: need BDSPAC) 1613* 1614 CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1615 $ 1, WORK( IR ), LDWRKR, CDUM, 1, 1616 $ RWORK( IRWORK ), INFO ) 1617* 1618* Multiply Q in U by left singular vectors of R in 1619* WORK(IR), storing result in A 1620* (CWorkspace: need N*N) 1621* (RWorkspace: 0) 1622* 1623 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 1624 $ WORK( IR ), LDWRKR, CZERO, A, LDA ) 1625* 1626* Copy left singular vectors of A from A to U 1627* 1628 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 1629* 1630 ELSE 1631* 1632* Insufficient workspace for a fast algorithm 1633* 1634 ITAU = 1 1635 IWORK = ITAU + N 1636* 1637* Compute A=Q*R, copying result to U 1638* (CWorkspace: need 2*N, prefer N+N*NB) 1639* (RWorkspace: 0) 1640* 1641 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1642 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1643 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1644* 1645* Generate Q in U 1646* (CWorkspace: need N+M, prefer N+M*NB) 1647* (RWorkspace: 0) 1648* 1649 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 1650 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1651 IE = 1 1652 ITAUQ = ITAU 1653 ITAUP = ITAUQ + N 1654 IWORK = ITAUP + N 1655* 1656* Zero out below R in A 1657* 1658 IF( N .GT. 1 ) THEN 1659 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1660 $ A( 2, 1 ), LDA ) 1661 END IF 1662* 1663* Bidiagonalize R in A 1664* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 1665* (RWorkspace: need N) 1666* 1667 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 1668 $ WORK( ITAUQ ), WORK( ITAUP ), 1669 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1670* 1671* Multiply Q in U by left bidiagonalizing vectors 1672* in A 1673* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1674* (RWorkspace: 0) 1675* 1676 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 1677 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 1678 $ LWORK-IWORK+1, IERR ) 1679 IRWORK = IE + N 1680* 1681* Perform bidiagonal QR iteration, computing left 1682* singular vectors of A in U 1683* (CWorkspace: 0) 1684* (RWorkspace: need BDSPAC) 1685* 1686 CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1687 $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ), 1688 $ INFO ) 1689* 1690 END IF 1691* 1692 ELSE IF( WNTVO ) THEN 1693* 1694* Path 8 (M much larger than N, JOBU='A', JOBVT='O') 1695* M left singular vectors to be computed in U and 1696* N right singular vectors to be overwritten on A 1697* 1698 IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN 1699* 1700* Sufficient workspace for a fast algorithm 1701* 1702 IU = 1 1703 IF( LWORK.GE.WRKBL+2*LDA*N ) THEN 1704* 1705* WORK(IU) is LDA by N and WORK(IR) is LDA by N 1706* 1707 LDWRKU = LDA 1708 IR = IU + LDWRKU*N 1709 LDWRKR = LDA 1710 ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN 1711* 1712* WORK(IU) is LDA by N and WORK(IR) is N by N 1713* 1714 LDWRKU = LDA 1715 IR = IU + LDWRKU*N 1716 LDWRKR = N 1717 ELSE 1718* 1719* WORK(IU) is N by N and WORK(IR) is N by N 1720* 1721 LDWRKU = N 1722 IR = IU + LDWRKU*N 1723 LDWRKR = N 1724 END IF 1725 ITAU = IR + LDWRKR*N 1726 IWORK = ITAU + N 1727* 1728* Compute A=Q*R, copying result to U 1729* (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 1730* (RWorkspace: 0) 1731* 1732 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1733 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1734 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1735* 1736* Generate Q in U 1737* (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) 1738* (RWorkspace: 0) 1739* 1740 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 1741 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1742* 1743* Copy R to WORK(IU), zeroing out below it 1744* 1745 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 1746 $ LDWRKU ) 1747 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1748 $ WORK( IU+1 ), LDWRKU ) 1749 IE = 1 1750 ITAUQ = ITAU 1751 ITAUP = ITAUQ + N 1752 IWORK = ITAUP + N 1753* 1754* Bidiagonalize R in WORK(IU), copying result to 1755* WORK(IR) 1756* (CWorkspace: need 2*N*N+3*N, 1757* prefer 2*N*N+2*N+2*N*NB) 1758* (RWorkspace: need N) 1759* 1760 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 1761 $ RWORK( IE ), WORK( ITAUQ ), 1762 $ WORK( ITAUP ), WORK( IWORK ), 1763 $ LWORK-IWORK+1, IERR ) 1764 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, 1765 $ WORK( IR ), LDWRKR ) 1766* 1767* Generate left bidiagonalizing vectors in WORK(IU) 1768* (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 1769* (RWorkspace: 0) 1770* 1771 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 1772 $ WORK( ITAUQ ), WORK( IWORK ), 1773 $ LWORK-IWORK+1, IERR ) 1774* 1775* Generate right bidiagonalizing vectors in WORK(IR) 1776* (CWorkspace: need 2*N*N+3*N-1, 1777* prefer 2*N*N+2*N+(N-1)*NB) 1778* (RWorkspace: 0) 1779* 1780 CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR, 1781 $ WORK( ITAUP ), WORK( IWORK ), 1782 $ LWORK-IWORK+1, IERR ) 1783 IRWORK = IE + N 1784* 1785* Perform bidiagonal QR iteration, computing left 1786* singular vectors of R in WORK(IU) and computing 1787* right singular vectors of R in WORK(IR) 1788* (CWorkspace: need 2*N*N) 1789* (RWorkspace: need BDSPAC) 1790* 1791 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), 1792 $ WORK( IR ), LDWRKR, WORK( IU ), 1793 $ LDWRKU, CDUM, 1, RWORK( IRWORK ), 1794 $ INFO ) 1795* 1796* Multiply Q in U by left singular vectors of R in 1797* WORK(IU), storing result in A 1798* (CWorkspace: need N*N) 1799* (RWorkspace: 0) 1800* 1801 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 1802 $ WORK( IU ), LDWRKU, CZERO, A, LDA ) 1803* 1804* Copy left singular vectors of A from A to U 1805* 1806 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 1807* 1808* Copy right singular vectors of R from WORK(IR) to A 1809* 1810 CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A, 1811 $ LDA ) 1812* 1813 ELSE 1814* 1815* Insufficient workspace for a fast algorithm 1816* 1817 ITAU = 1 1818 IWORK = ITAU + N 1819* 1820* Compute A=Q*R, copying result to U 1821* (CWorkspace: need 2*N, prefer N+N*NB) 1822* (RWorkspace: 0) 1823* 1824 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1825 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1826 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1827* 1828* Generate Q in U 1829* (CWorkspace: need N+M, prefer N+M*NB) 1830* (RWorkspace: 0) 1831* 1832 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 1833 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1834 IE = 1 1835 ITAUQ = ITAU 1836 ITAUP = ITAUQ + N 1837 IWORK = ITAUP + N 1838* 1839* Zero out below R in A 1840* 1841 IF( N .GT. 1 ) THEN 1842 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1843 $ A( 2, 1 ), LDA ) 1844 END IF 1845* 1846* Bidiagonalize R in A 1847* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 1848* (RWorkspace: need N) 1849* 1850 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 1851 $ WORK( ITAUQ ), WORK( ITAUP ), 1852 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1853* 1854* Multiply Q in U by left bidiagonalizing vectors 1855* in A 1856* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 1857* (RWorkspace: 0) 1858* 1859 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 1860 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 1861 $ LWORK-IWORK+1, IERR ) 1862* 1863* Generate right bidiagonalizing vectors in A 1864* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 1865* (RWorkspace: 0) 1866* 1867 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 1868 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1869 IRWORK = IE + N 1870* 1871* Perform bidiagonal QR iteration, computing left 1872* singular vectors of A in U and computing right 1873* singular vectors of A in A 1874* (CWorkspace: 0) 1875* (RWorkspace: need BDSPAC) 1876* 1877 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A, 1878 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 1879 $ INFO ) 1880* 1881 END IF 1882* 1883 ELSE IF( WNTVAS ) THEN 1884* 1885* Path 9 (M much larger than N, JOBU='A', JOBVT='S' 1886* or 'A') 1887* M left singular vectors to be computed in U and 1888* N right singular vectors to be computed in VT 1889* 1890 IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN 1891* 1892* Sufficient workspace for a fast algorithm 1893* 1894 IU = 1 1895 IF( LWORK.GE.WRKBL+LDA*N ) THEN 1896* 1897* WORK(IU) is LDA by N 1898* 1899 LDWRKU = LDA 1900 ELSE 1901* 1902* WORK(IU) is N by N 1903* 1904 LDWRKU = N 1905 END IF 1906 ITAU = IU + LDWRKU*N 1907 IWORK = ITAU + N 1908* 1909* Compute A=Q*R, copying result to U 1910* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 1911* (RWorkspace: 0) 1912* 1913 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1914 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1915 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1916* 1917* Generate Q in U 1918* (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) 1919* (RWorkspace: 0) 1920* 1921 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 1922 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1923* 1924* Copy R to WORK(IU), zeroing out below it 1925* 1926 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 1927 $ LDWRKU ) 1928 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 1929 $ WORK( IU+1 ), LDWRKU ) 1930 IE = 1 1931 ITAUQ = ITAU 1932 ITAUP = ITAUQ + N 1933 IWORK = ITAUP + N 1934* 1935* Bidiagonalize R in WORK(IU), copying result to VT 1936* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 1937* (RWorkspace: need N) 1938* 1939 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 1940 $ RWORK( IE ), WORK( ITAUQ ), 1941 $ WORK( ITAUP ), WORK( IWORK ), 1942 $ LWORK-IWORK+1, IERR ) 1943 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT, 1944 $ LDVT ) 1945* 1946* Generate left bidiagonalizing vectors in WORK(IU) 1947* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 1948* (RWorkspace: 0) 1949* 1950 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 1951 $ WORK( ITAUQ ), WORK( IWORK ), 1952 $ LWORK-IWORK+1, IERR ) 1953* 1954* Generate right bidiagonalizing vectors in VT 1955* (CWorkspace: need N*N+3*N-1, 1956* prefer N*N+2*N+(N-1)*NB) 1957* (RWorkspace: need 0) 1958* 1959 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 1960 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1961 IRWORK = IE + N 1962* 1963* Perform bidiagonal QR iteration, computing left 1964* singular vectors of R in WORK(IU) and computing 1965* right singular vectors of R in VT 1966* (CWorkspace: need N*N) 1967* (RWorkspace: need BDSPAC) 1968* 1969 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 1970 $ LDVT, WORK( IU ), LDWRKU, CDUM, 1, 1971 $ RWORK( IRWORK ), INFO ) 1972* 1973* Multiply Q in U by left singular vectors of R in 1974* WORK(IU), storing result in A 1975* (CWorkspace: need N*N) 1976* (RWorkspace: 0) 1977* 1978 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 1979 $ WORK( IU ), LDWRKU, CZERO, A, LDA ) 1980* 1981* Copy left singular vectors of A from A to U 1982* 1983 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 1984* 1985 ELSE 1986* 1987* Insufficient workspace for a fast algorithm 1988* 1989 ITAU = 1 1990 IWORK = ITAU + N 1991* 1992* Compute A=Q*R, copying result to U 1993* (CWorkspace: need 2*N, prefer N+N*NB) 1994* (RWorkspace: 0) 1995* 1996 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 1997 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 1998 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 1999* 2000* Generate Q in U 2001* (CWorkspace: need N+M, prefer N+M*NB) 2002* (RWorkspace: 0) 2003* 2004 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 2005 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2006* 2007* Copy R from A to VT, zeroing out below it 2008* 2009 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 2010 IF( N.GT.1 ) 2011 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 2012 $ VT( 2, 1 ), LDVT ) 2013 IE = 1 2014 ITAUQ = ITAU 2015 ITAUP = ITAUQ + N 2016 IWORK = ITAUP + N 2017* 2018* Bidiagonalize R in VT 2019* (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 2020* (RWorkspace: need N) 2021* 2022 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 2023 $ WORK( ITAUQ ), WORK( ITAUP ), 2024 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2025* 2026* Multiply Q in U by left bidiagonalizing vectors 2027* in VT 2028* (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 2029* (RWorkspace: 0) 2030* 2031 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 2032 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 2033 $ LWORK-IWORK+1, IERR ) 2034* 2035* Generate right bidiagonalizing vectors in VT 2036* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 2037* (RWorkspace: 0) 2038* 2039 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 2040 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2041 IRWORK = IE + N 2042* 2043* Perform bidiagonal QR iteration, computing left 2044* singular vectors of A in U and computing right 2045* singular vectors of A in VT 2046* (CWorkspace: 0) 2047* (RWorkspace: need BDSPAC) 2048* 2049 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 2050 $ LDVT, U, LDU, CDUM, 1, 2051 $ RWORK( IRWORK ), INFO ) 2052* 2053 END IF 2054* 2055 END IF 2056* 2057 END IF 2058* 2059 ELSE 2060* 2061* M .LT. MNTHR 2062* 2063* Path 10 (M at least N, but not much larger) 2064* Reduce to bidiagonal form without QR decomposition 2065* 2066 IE = 1 2067 ITAUQ = 1 2068 ITAUP = ITAUQ + N 2069 IWORK = ITAUP + N 2070* 2071* Bidiagonalize A 2072* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 2073* (RWorkspace: need N) 2074* 2075 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 2076 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 2077 $ IERR ) 2078 IF( WNTUAS ) THEN 2079* 2080* If left singular vectors desired in U, copy result to U 2081* and generate left bidiagonalizing vectors in U 2082* (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB) 2083* (RWorkspace: 0) 2084* 2085 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 2086 IF( WNTUS ) 2087 $ NCU = N 2088 IF( WNTUA ) 2089 $ NCU = M 2090 CALL ZUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ), 2091 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2092 END IF 2093 IF( WNTVAS ) THEN 2094* 2095* If right singular vectors desired in VT, copy result to 2096* VT and generate right bidiagonalizing vectors in VT 2097* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 2098* (RWorkspace: 0) 2099* 2100 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 2101 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 2102 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2103 END IF 2104 IF( WNTUO ) THEN 2105* 2106* If left singular vectors desired in A, generate left 2107* bidiagonalizing vectors in A 2108* (CWorkspace: need 3*N, prefer 2*N+N*NB) 2109* (RWorkspace: 0) 2110* 2111 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 2112 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2113 END IF 2114 IF( WNTVO ) THEN 2115* 2116* If right singular vectors desired in A, generate right 2117* bidiagonalizing vectors in A 2118* (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 2119* (RWorkspace: 0) 2120* 2121 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 2122 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2123 END IF 2124 IRWORK = IE + N 2125 IF( WNTUAS .OR. WNTUO ) 2126 $ NRU = M 2127 IF( WNTUN ) 2128 $ NRU = 0 2129 IF( WNTVAS .OR. WNTVO ) 2130 $ NCVT = N 2131 IF( WNTVN ) 2132 $ NCVT = 0 2133 IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN 2134* 2135* Perform bidiagonal QR iteration, if desired, computing 2136* left singular vectors in U and computing right singular 2137* vectors in VT 2138* (CWorkspace: 0) 2139* (RWorkspace: need BDSPAC) 2140* 2141 CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT, 2142 $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ), 2143 $ INFO ) 2144 ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN 2145* 2146* Perform bidiagonal QR iteration, if desired, computing 2147* left singular vectors in U and computing right singular 2148* vectors in A 2149* (CWorkspace: 0) 2150* (RWorkspace: need BDSPAC) 2151* 2152 CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A, 2153 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 2154 $ INFO ) 2155 ELSE 2156* 2157* Perform bidiagonal QR iteration, if desired, computing 2158* left singular vectors in A and computing right singular 2159* vectors in VT 2160* (CWorkspace: 0) 2161* (RWorkspace: need BDSPAC) 2162* 2163 CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT, 2164 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 2165 $ INFO ) 2166 END IF 2167* 2168 END IF 2169* 2170 ELSE 2171* 2172* A has more columns than rows. If A has sufficiently more 2173* columns than rows, first reduce using the LQ decomposition (if 2174* sufficient workspace available) 2175* 2176 IF( N.GE.MNTHR ) THEN 2177* 2178 IF( WNTVN ) THEN 2179* 2180* Path 1t(N much larger than M, JOBVT='N') 2181* No right singular vectors to be computed 2182* 2183 ITAU = 1 2184 IWORK = ITAU + M 2185* 2186* Compute A=L*Q 2187* (CWorkspace: need 2*M, prefer M+M*NB) 2188* (RWorkspace: 0) 2189* 2190 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 2191 $ LWORK-IWORK+1, IERR ) 2192* 2193* Zero out above L 2194* 2195 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 2196 $ LDA ) 2197 IE = 1 2198 ITAUQ = 1 2199 ITAUP = ITAUQ + M 2200 IWORK = ITAUP + M 2201* 2202* Bidiagonalize L in A 2203* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 2204* (RWorkspace: need M) 2205* 2206 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 2207 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 2208 $ IERR ) 2209 IF( WNTUO .OR. WNTUAS ) THEN 2210* 2211* If left singular vectors desired, generate Q 2212* (CWorkspace: need 3*M, prefer 2*M+M*NB) 2213* (RWorkspace: 0) 2214* 2215 CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 2216 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2217 END IF 2218 IRWORK = IE + M 2219 NRU = 0 2220 IF( WNTUO .OR. WNTUAS ) 2221 $ NRU = M 2222* 2223* Perform bidiagonal QR iteration, computing left singular 2224* vectors of A in A if desired 2225* (CWorkspace: 0) 2226* (RWorkspace: need BDSPAC) 2227* 2228 CALL ZBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1, 2229 $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO ) 2230* 2231* If left singular vectors desired in U, copy them there 2232* 2233 IF( WNTUAS ) 2234 $ CALL ZLACPY( 'F', M, M, A, LDA, U, LDU ) 2235* 2236 ELSE IF( WNTVO .AND. WNTUN ) THEN 2237* 2238* Path 2t(N much larger than M, JOBU='N', JOBVT='O') 2239* M right singular vectors to be overwritten on A and 2240* no left singular vectors to be computed 2241* 2242 IF( LWORK.GE.M*M+3*M ) THEN 2243* 2244* Sufficient workspace for a fast algorithm 2245* 2246 IR = 1 2247 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN 2248* 2249* WORK(IU) is LDA by N and WORK(IR) is LDA by M 2250* 2251 LDWRKU = LDA 2252 CHUNK = N 2253 LDWRKR = LDA 2254 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN 2255* 2256* WORK(IU) is LDA by N and WORK(IR) is M by M 2257* 2258 LDWRKU = LDA 2259 CHUNK = N 2260 LDWRKR = M 2261 ELSE 2262* 2263* WORK(IU) is M by CHUNK and WORK(IR) is M by M 2264* 2265 LDWRKU = M 2266 CHUNK = ( LWORK-M*M ) / M 2267 LDWRKR = M 2268 END IF 2269 ITAU = IR + LDWRKR*M 2270 IWORK = ITAU + M 2271* 2272* Compute A=L*Q 2273* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2274* (RWorkspace: 0) 2275* 2276 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2277 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2278* 2279* Copy L to WORK(IR) and zero out above it 2280* 2281 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR ) 2282 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 2283 $ WORK( IR+LDWRKR ), LDWRKR ) 2284* 2285* Generate Q in A 2286* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2287* (RWorkspace: 0) 2288* 2289 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 2290 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2291 IE = 1 2292 ITAUQ = ITAU 2293 ITAUP = ITAUQ + M 2294 IWORK = ITAUP + M 2295* 2296* Bidiagonalize L in WORK(IR) 2297* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 2298* (RWorkspace: need M) 2299* 2300 CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ), 2301 $ WORK( ITAUQ ), WORK( ITAUP ), 2302 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2303* 2304* Generate right vectors bidiagonalizing L 2305* (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) 2306* (RWorkspace: 0) 2307* 2308 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 2309 $ WORK( ITAUP ), WORK( IWORK ), 2310 $ LWORK-IWORK+1, IERR ) 2311 IRWORK = IE + M 2312* 2313* Perform bidiagonal QR iteration, computing right 2314* singular vectors of L in WORK(IR) 2315* (CWorkspace: need M*M) 2316* (RWorkspace: need BDSPAC) 2317* 2318 CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 2319 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 2320 $ RWORK( IRWORK ), INFO ) 2321 IU = ITAUQ 2322* 2323* Multiply right singular vectors of L in WORK(IR) by Q 2324* in A, storing result in WORK(IU) and copying to A 2325* (CWorkspace: need M*M+M, prefer M*M+M*N) 2326* (RWorkspace: 0) 2327* 2328 DO 30 I = 1, N, CHUNK 2329 BLK = MIN( N-I+1, CHUNK ) 2330 CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ), 2331 $ LDWRKR, A( 1, I ), LDA, CZERO, 2332 $ WORK( IU ), LDWRKU ) 2333 CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU, 2334 $ A( 1, I ), LDA ) 2335 30 CONTINUE 2336* 2337 ELSE 2338* 2339* Insufficient workspace for a fast algorithm 2340* 2341 IE = 1 2342 ITAUQ = 1 2343 ITAUP = ITAUQ + M 2344 IWORK = ITAUP + M 2345* 2346* Bidiagonalize A 2347* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 2348* (RWorkspace: need M) 2349* 2350 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), 2351 $ WORK( ITAUQ ), WORK( ITAUP ), 2352 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2353* 2354* Generate right vectors bidiagonalizing A 2355* (CWorkspace: need 3*M, prefer 2*M+M*NB) 2356* (RWorkspace: 0) 2357* 2358 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 2359 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2360 IRWORK = IE + M 2361* 2362* Perform bidiagonal QR iteration, computing right 2363* singular vectors of A in A 2364* (CWorkspace: 0) 2365* (RWorkspace: need BDSPAC) 2366* 2367 CALL ZBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA, 2368 $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO ) 2369* 2370 END IF 2371* 2372 ELSE IF( WNTVO .AND. WNTUAS ) THEN 2373* 2374* Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O') 2375* M right singular vectors to be overwritten on A and 2376* M left singular vectors to be computed in U 2377* 2378 IF( LWORK.GE.M*M+3*M ) THEN 2379* 2380* Sufficient workspace for a fast algorithm 2381* 2382 IR = 1 2383 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN 2384* 2385* WORK(IU) is LDA by N and WORK(IR) is LDA by M 2386* 2387 LDWRKU = LDA 2388 CHUNK = N 2389 LDWRKR = LDA 2390 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN 2391* 2392* WORK(IU) is LDA by N and WORK(IR) is M by M 2393* 2394 LDWRKU = LDA 2395 CHUNK = N 2396 LDWRKR = M 2397 ELSE 2398* 2399* WORK(IU) is M by CHUNK and WORK(IR) is M by M 2400* 2401 LDWRKU = M 2402 CHUNK = ( LWORK-M*M ) / M 2403 LDWRKR = M 2404 END IF 2405 ITAU = IR + LDWRKR*M 2406 IWORK = ITAU + M 2407* 2408* Compute A=L*Q 2409* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2410* (RWorkspace: 0) 2411* 2412 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2413 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2414* 2415* Copy L to U, zeroing about above it 2416* 2417 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 2418 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ), 2419 $ LDU ) 2420* 2421* Generate Q in A 2422* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2423* (RWorkspace: 0) 2424* 2425 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 2426 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2427 IE = 1 2428 ITAUQ = ITAU 2429 ITAUP = ITAUQ + M 2430 IWORK = ITAUP + M 2431* 2432* Bidiagonalize L in U, copying result to WORK(IR) 2433* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 2434* (RWorkspace: need M) 2435* 2436 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 2437 $ WORK( ITAUQ ), WORK( ITAUP ), 2438 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2439 CALL ZLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) 2440* 2441* Generate right vectors bidiagonalizing L in WORK(IR) 2442* (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) 2443* (RWorkspace: 0) 2444* 2445 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 2446 $ WORK( ITAUP ), WORK( IWORK ), 2447 $ LWORK-IWORK+1, IERR ) 2448* 2449* Generate left vectors bidiagonalizing L in U 2450* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 2451* (RWorkspace: 0) 2452* 2453 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 2454 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2455 IRWORK = IE + M 2456* 2457* Perform bidiagonal QR iteration, computing left 2458* singular vectors of L in U, and computing right 2459* singular vectors of L in WORK(IR) 2460* (CWorkspace: need M*M) 2461* (RWorkspace: need BDSPAC) 2462* 2463 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 2464 $ WORK( IR ), LDWRKR, U, LDU, CDUM, 1, 2465 $ RWORK( IRWORK ), INFO ) 2466 IU = ITAUQ 2467* 2468* Multiply right singular vectors of L in WORK(IR) by Q 2469* in A, storing result in WORK(IU) and copying to A 2470* (CWorkspace: need M*M+M, prefer M*M+M*N)) 2471* (RWorkspace: 0) 2472* 2473 DO 40 I = 1, N, CHUNK 2474 BLK = MIN( N-I+1, CHUNK ) 2475 CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ), 2476 $ LDWRKR, A( 1, I ), LDA, CZERO, 2477 $ WORK( IU ), LDWRKU ) 2478 CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU, 2479 $ A( 1, I ), LDA ) 2480 40 CONTINUE 2481* 2482 ELSE 2483* 2484* Insufficient workspace for a fast algorithm 2485* 2486 ITAU = 1 2487 IWORK = ITAU + M 2488* 2489* Compute A=L*Q 2490* (CWorkspace: need 2*M, prefer M+M*NB) 2491* (RWorkspace: 0) 2492* 2493 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2494 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2495* 2496* Copy L to U, zeroing out above it 2497* 2498 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 2499 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ), 2500 $ LDU ) 2501* 2502* Generate Q in A 2503* (CWorkspace: need 2*M, prefer M+M*NB) 2504* (RWorkspace: 0) 2505* 2506 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 2507 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2508 IE = 1 2509 ITAUQ = ITAU 2510 ITAUP = ITAUQ + M 2511 IWORK = ITAUP + M 2512* 2513* Bidiagonalize L in U 2514* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 2515* (RWorkspace: need M) 2516* 2517 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 2518 $ WORK( ITAUQ ), WORK( ITAUP ), 2519 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2520* 2521* Multiply right vectors bidiagonalizing L by Q in A 2522* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 2523* (RWorkspace: 0) 2524* 2525 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 2526 $ WORK( ITAUP ), A, LDA, WORK( IWORK ), 2527 $ LWORK-IWORK+1, IERR ) 2528* 2529* Generate left vectors bidiagonalizing L in U 2530* (CWorkspace: need 3*M, prefer 2*M+M*NB) 2531* (RWorkspace: 0) 2532* 2533 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 2534 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2535 IRWORK = IE + M 2536* 2537* Perform bidiagonal QR iteration, computing left 2538* singular vectors of A in U and computing right 2539* singular vectors of A in A 2540* (CWorkspace: 0) 2541* (RWorkspace: need BDSPAC) 2542* 2543 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA, 2544 $ U, LDU, CDUM, 1, RWORK( IRWORK ), INFO ) 2545* 2546 END IF 2547* 2548 ELSE IF( WNTVS ) THEN 2549* 2550 IF( WNTUN ) THEN 2551* 2552* Path 4t(N much larger than M, JOBU='N', JOBVT='S') 2553* M right singular vectors to be computed in VT and 2554* no left singular vectors to be computed 2555* 2556 IF( LWORK.GE.M*M+3*M ) THEN 2557* 2558* Sufficient workspace for a fast algorithm 2559* 2560 IR = 1 2561 IF( LWORK.GE.WRKBL+LDA*M ) THEN 2562* 2563* WORK(IR) is LDA by M 2564* 2565 LDWRKR = LDA 2566 ELSE 2567* 2568* WORK(IR) is M by M 2569* 2570 LDWRKR = M 2571 END IF 2572 ITAU = IR + LDWRKR*M 2573 IWORK = ITAU + M 2574* 2575* Compute A=L*Q 2576* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2577* (RWorkspace: 0) 2578* 2579 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2580 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2581* 2582* Copy L to WORK(IR), zeroing out above it 2583* 2584 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), 2585 $ LDWRKR ) 2586 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 2587 $ WORK( IR+LDWRKR ), LDWRKR ) 2588* 2589* Generate Q in A 2590* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2591* (RWorkspace: 0) 2592* 2593 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 2594 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2595 IE = 1 2596 ITAUQ = ITAU 2597 ITAUP = ITAUQ + M 2598 IWORK = ITAUP + M 2599* 2600* Bidiagonalize L in WORK(IR) 2601* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 2602* (RWorkspace: need M) 2603* 2604 CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, 2605 $ RWORK( IE ), WORK( ITAUQ ), 2606 $ WORK( ITAUP ), WORK( IWORK ), 2607 $ LWORK-IWORK+1, IERR ) 2608* 2609* Generate right vectors bidiagonalizing L in 2610* WORK(IR) 2611* (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) 2612* (RWorkspace: 0) 2613* 2614 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 2615 $ WORK( ITAUP ), WORK( IWORK ), 2616 $ LWORK-IWORK+1, IERR ) 2617 IRWORK = IE + M 2618* 2619* Perform bidiagonal QR iteration, computing right 2620* singular vectors of L in WORK(IR) 2621* (CWorkspace: need M*M) 2622* (RWorkspace: need BDSPAC) 2623* 2624 CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 2625 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 2626 $ RWORK( IRWORK ), INFO ) 2627* 2628* Multiply right singular vectors of L in WORK(IR) by 2629* Q in A, storing result in VT 2630* (CWorkspace: need M*M) 2631* (RWorkspace: 0) 2632* 2633 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ), 2634 $ LDWRKR, A, LDA, CZERO, VT, LDVT ) 2635* 2636 ELSE 2637* 2638* Insufficient workspace for a fast algorithm 2639* 2640 ITAU = 1 2641 IWORK = ITAU + M 2642* 2643* Compute A=L*Q 2644* (CWorkspace: need 2*M, prefer M+M*NB) 2645* (RWorkspace: 0) 2646* 2647 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2648 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2649* 2650* Copy result to VT 2651* 2652 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 2653* 2654* Generate Q in VT 2655* (CWorkspace: need 2*M, prefer M+M*NB) 2656* (RWorkspace: 0) 2657* 2658 CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 2659 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2660 IE = 1 2661 ITAUQ = ITAU 2662 ITAUP = ITAUQ + M 2663 IWORK = ITAUP + M 2664* 2665* Zero out above L in A 2666* 2667 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 2668 $ A( 1, 2 ), LDA ) 2669* 2670* Bidiagonalize L in A 2671* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 2672* (RWorkspace: need M) 2673* 2674 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 2675 $ WORK( ITAUQ ), WORK( ITAUP ), 2676 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2677* 2678* Multiply right vectors bidiagonalizing L by Q in VT 2679* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 2680* (RWorkspace: 0) 2681* 2682 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 2683 $ WORK( ITAUP ), VT, LDVT, 2684 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2685 IRWORK = IE + M 2686* 2687* Perform bidiagonal QR iteration, computing right 2688* singular vectors of A in VT 2689* (CWorkspace: 0) 2690* (RWorkspace: need BDSPAC) 2691* 2692 CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT, 2693 $ LDVT, CDUM, 1, CDUM, 1, 2694 $ RWORK( IRWORK ), INFO ) 2695* 2696 END IF 2697* 2698 ELSE IF( WNTUO ) THEN 2699* 2700* Path 5t(N much larger than M, JOBU='O', JOBVT='S') 2701* M right singular vectors to be computed in VT and 2702* M left singular vectors to be overwritten on A 2703* 2704 IF( LWORK.GE.2*M*M+3*M ) THEN 2705* 2706* Sufficient workspace for a fast algorithm 2707* 2708 IU = 1 2709 IF( LWORK.GE.WRKBL+2*LDA*M ) THEN 2710* 2711* WORK(IU) is LDA by M and WORK(IR) is LDA by M 2712* 2713 LDWRKU = LDA 2714 IR = IU + LDWRKU*M 2715 LDWRKR = LDA 2716 ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN 2717* 2718* WORK(IU) is LDA by M and WORK(IR) is M by M 2719* 2720 LDWRKU = LDA 2721 IR = IU + LDWRKU*M 2722 LDWRKR = M 2723 ELSE 2724* 2725* WORK(IU) is M by M and WORK(IR) is M by M 2726* 2727 LDWRKU = M 2728 IR = IU + LDWRKU*M 2729 LDWRKR = M 2730 END IF 2731 ITAU = IR + LDWRKR*M 2732 IWORK = ITAU + M 2733* 2734* Compute A=L*Q 2735* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 2736* (RWorkspace: 0) 2737* 2738 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2739 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2740* 2741* Copy L to WORK(IU), zeroing out below it 2742* 2743 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 2744 $ LDWRKU ) 2745 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 2746 $ WORK( IU+LDWRKU ), LDWRKU ) 2747* 2748* Generate Q in A 2749* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 2750* (RWorkspace: 0) 2751* 2752 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 2753 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2754 IE = 1 2755 ITAUQ = ITAU 2756 ITAUP = ITAUQ + M 2757 IWORK = ITAUP + M 2758* 2759* Bidiagonalize L in WORK(IU), copying result to 2760* WORK(IR) 2761* (CWorkspace: need 2*M*M+3*M, 2762* prefer 2*M*M+2*M+2*M*NB) 2763* (RWorkspace: need M) 2764* 2765 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 2766 $ RWORK( IE ), WORK( ITAUQ ), 2767 $ WORK( ITAUP ), WORK( IWORK ), 2768 $ LWORK-IWORK+1, IERR ) 2769 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, 2770 $ WORK( IR ), LDWRKR ) 2771* 2772* Generate right bidiagonalizing vectors in WORK(IU) 2773* (CWorkspace: need 2*M*M+3*M-1, 2774* prefer 2*M*M+2*M+(M-1)*NB) 2775* (RWorkspace: 0) 2776* 2777 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 2778 $ WORK( ITAUP ), WORK( IWORK ), 2779 $ LWORK-IWORK+1, IERR ) 2780* 2781* Generate left bidiagonalizing vectors in WORK(IR) 2782* (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 2783* (RWorkspace: 0) 2784* 2785 CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, 2786 $ WORK( ITAUQ ), WORK( IWORK ), 2787 $ LWORK-IWORK+1, IERR ) 2788 IRWORK = IE + M 2789* 2790* Perform bidiagonal QR iteration, computing left 2791* singular vectors of L in WORK(IR) and computing 2792* right singular vectors of L in WORK(IU) 2793* (CWorkspace: need 2*M*M) 2794* (RWorkspace: need BDSPAC) 2795* 2796 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 2797 $ WORK( IU ), LDWRKU, WORK( IR ), 2798 $ LDWRKR, CDUM, 1, RWORK( IRWORK ), 2799 $ INFO ) 2800* 2801* Multiply right singular vectors of L in WORK(IU) by 2802* Q in A, storing result in VT 2803* (CWorkspace: need M*M) 2804* (RWorkspace: 0) 2805* 2806 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 2807 $ LDWRKU, A, LDA, CZERO, VT, LDVT ) 2808* 2809* Copy left singular vectors of L to A 2810* (CWorkspace: need M*M) 2811* (RWorkspace: 0) 2812* 2813 CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A, 2814 $ LDA ) 2815* 2816 ELSE 2817* 2818* Insufficient workspace for a fast algorithm 2819* 2820 ITAU = 1 2821 IWORK = ITAU + M 2822* 2823* Compute A=L*Q, copying result to VT 2824* (CWorkspace: need 2*M, prefer M+M*NB) 2825* (RWorkspace: 0) 2826* 2827 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2828 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2829 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 2830* 2831* Generate Q in VT 2832* (CWorkspace: need 2*M, prefer M+M*NB) 2833* (RWorkspace: 0) 2834* 2835 CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 2836 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2837 IE = 1 2838 ITAUQ = ITAU 2839 ITAUP = ITAUQ + M 2840 IWORK = ITAUP + M 2841* 2842* Zero out above L in A 2843* 2844 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 2845 $ A( 1, 2 ), LDA ) 2846* 2847* Bidiagonalize L in A 2848* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 2849* (RWorkspace: need M) 2850* 2851 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 2852 $ WORK( ITAUQ ), WORK( ITAUP ), 2853 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2854* 2855* Multiply right vectors bidiagonalizing L by Q in VT 2856* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 2857* (RWorkspace: 0) 2858* 2859 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 2860 $ WORK( ITAUP ), VT, LDVT, 2861 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2862* 2863* Generate left bidiagonalizing vectors of L in A 2864* (CWorkspace: need 3*M, prefer 2*M+M*NB) 2865* (RWorkspace: 0) 2866* 2867 CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 2868 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2869 IRWORK = IE + M 2870* 2871* Perform bidiagonal QR iteration, computing left 2872* singular vectors of A in A and computing right 2873* singular vectors of A in VT 2874* (CWorkspace: 0) 2875* (RWorkspace: need BDSPAC) 2876* 2877 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 2878 $ LDVT, A, LDA, CDUM, 1, 2879 $ RWORK( IRWORK ), INFO ) 2880* 2881 END IF 2882* 2883 ELSE IF( WNTUAS ) THEN 2884* 2885* Path 6t(N much larger than M, JOBU='S' or 'A', 2886* JOBVT='S') 2887* M right singular vectors to be computed in VT and 2888* M left singular vectors to be computed in U 2889* 2890 IF( LWORK.GE.M*M+3*M ) THEN 2891* 2892* Sufficient workspace for a fast algorithm 2893* 2894 IU = 1 2895 IF( LWORK.GE.WRKBL+LDA*M ) THEN 2896* 2897* WORK(IU) is LDA by N 2898* 2899 LDWRKU = LDA 2900 ELSE 2901* 2902* WORK(IU) is LDA by M 2903* 2904 LDWRKU = M 2905 END IF 2906 ITAU = IU + LDWRKU*M 2907 IWORK = ITAU + M 2908* 2909* Compute A=L*Q 2910* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2911* (RWorkspace: 0) 2912* 2913 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2914 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2915* 2916* Copy L to WORK(IU), zeroing out above it 2917* 2918 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 2919 $ LDWRKU ) 2920 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 2921 $ WORK( IU+LDWRKU ), LDWRKU ) 2922* 2923* Generate Q in A 2924* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 2925* (RWorkspace: 0) 2926* 2927 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 2928 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2929 IE = 1 2930 ITAUQ = ITAU 2931 ITAUP = ITAUQ + M 2932 IWORK = ITAUP + M 2933* 2934* Bidiagonalize L in WORK(IU), copying result to U 2935* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 2936* (RWorkspace: need M) 2937* 2938 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 2939 $ RWORK( IE ), WORK( ITAUQ ), 2940 $ WORK( ITAUP ), WORK( IWORK ), 2941 $ LWORK-IWORK+1, IERR ) 2942 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U, 2943 $ LDU ) 2944* 2945* Generate right bidiagonalizing vectors in WORK(IU) 2946* (CWorkspace: need M*M+3*M-1, 2947* prefer M*M+2*M+(M-1)*NB) 2948* (RWorkspace: 0) 2949* 2950 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 2951 $ WORK( ITAUP ), WORK( IWORK ), 2952 $ LWORK-IWORK+1, IERR ) 2953* 2954* Generate left bidiagonalizing vectors in U 2955* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 2956* (RWorkspace: 0) 2957* 2958 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 2959 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2960 IRWORK = IE + M 2961* 2962* Perform bidiagonal QR iteration, computing left 2963* singular vectors of L in U and computing right 2964* singular vectors of L in WORK(IU) 2965* (CWorkspace: need M*M) 2966* (RWorkspace: need BDSPAC) 2967* 2968 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 2969 $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1, 2970 $ RWORK( IRWORK ), INFO ) 2971* 2972* Multiply right singular vectors of L in WORK(IU) by 2973* Q in A, storing result in VT 2974* (CWorkspace: need M*M) 2975* (RWorkspace: 0) 2976* 2977 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 2978 $ LDWRKU, A, LDA, CZERO, VT, LDVT ) 2979* 2980 ELSE 2981* 2982* Insufficient workspace for a fast algorithm 2983* 2984 ITAU = 1 2985 IWORK = ITAU + M 2986* 2987* Compute A=L*Q, copying result to VT 2988* (CWorkspace: need 2*M, prefer M+M*NB) 2989* (RWorkspace: 0) 2990* 2991 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 2992 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 2993 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 2994* 2995* Generate Q in VT 2996* (CWorkspace: need 2*M, prefer M+M*NB) 2997* (RWorkspace: 0) 2998* 2999 CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 3000 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3001* 3002* Copy L to U, zeroing out above it 3003* 3004 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 3005 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3006 $ U( 1, 2 ), LDU ) 3007 IE = 1 3008 ITAUQ = ITAU 3009 ITAUP = ITAUQ + M 3010 IWORK = ITAUP + M 3011* 3012* Bidiagonalize L in U 3013* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 3014* (RWorkspace: need M) 3015* 3016 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 3017 $ WORK( ITAUQ ), WORK( ITAUP ), 3018 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3019* 3020* Multiply right bidiagonalizing vectors in U by Q 3021* in VT 3022* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 3023* (RWorkspace: 0) 3024* 3025 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 3026 $ WORK( ITAUP ), VT, LDVT, 3027 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3028* 3029* Generate left bidiagonalizing vectors in U 3030* (CWorkspace: need 3*M, prefer 2*M+M*NB) 3031* (RWorkspace: 0) 3032* 3033 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 3034 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3035 IRWORK = IE + M 3036* 3037* Perform bidiagonal QR iteration, computing left 3038* singular vectors of A in U and computing right 3039* singular vectors of A in VT 3040* (CWorkspace: 0) 3041* (RWorkspace: need BDSPAC) 3042* 3043 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 3044 $ LDVT, U, LDU, CDUM, 1, 3045 $ RWORK( IRWORK ), INFO ) 3046* 3047 END IF 3048* 3049 END IF 3050* 3051 ELSE IF( WNTVA ) THEN 3052* 3053 IF( WNTUN ) THEN 3054* 3055* Path 7t(N much larger than M, JOBU='N', JOBVT='A') 3056* N right singular vectors to be computed in VT and 3057* no left singular vectors to be computed 3058* 3059 IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN 3060* 3061* Sufficient workspace for a fast algorithm 3062* 3063 IR = 1 3064 IF( LWORK.GE.WRKBL+LDA*M ) THEN 3065* 3066* WORK(IR) is LDA by M 3067* 3068 LDWRKR = LDA 3069 ELSE 3070* 3071* WORK(IR) is M by M 3072* 3073 LDWRKR = M 3074 END IF 3075 ITAU = IR + LDWRKR*M 3076 IWORK = ITAU + M 3077* 3078* Compute A=L*Q, copying result to VT 3079* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 3080* (RWorkspace: 0) 3081* 3082 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 3083 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3084 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3085* 3086* Copy L to WORK(IR), zeroing out above it 3087* 3088 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), 3089 $ LDWRKR ) 3090 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3091 $ WORK( IR+LDWRKR ), LDWRKR ) 3092* 3093* Generate Q in VT 3094* (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) 3095* (RWorkspace: 0) 3096* 3097 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 3098 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3099 IE = 1 3100 ITAUQ = ITAU 3101 ITAUP = ITAUQ + M 3102 IWORK = ITAUP + M 3103* 3104* Bidiagonalize L in WORK(IR) 3105* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 3106* (RWorkspace: need M) 3107* 3108 CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, 3109 $ RWORK( IE ), WORK( ITAUQ ), 3110 $ WORK( ITAUP ), WORK( IWORK ), 3111 $ LWORK-IWORK+1, IERR ) 3112* 3113* Generate right bidiagonalizing vectors in WORK(IR) 3114* (CWorkspace: need M*M+3*M-1, 3115* prefer M*M+2*M+(M-1)*NB) 3116* (RWorkspace: 0) 3117* 3118 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 3119 $ WORK( ITAUP ), WORK( IWORK ), 3120 $ LWORK-IWORK+1, IERR ) 3121 IRWORK = IE + M 3122* 3123* Perform bidiagonal QR iteration, computing right 3124* singular vectors of L in WORK(IR) 3125* (CWorkspace: need M*M) 3126* (RWorkspace: need BDSPAC) 3127* 3128 CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 3129 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 3130 $ RWORK( IRWORK ), INFO ) 3131* 3132* Multiply right singular vectors of L in WORK(IR) by 3133* Q in VT, storing result in A 3134* (CWorkspace: need M*M) 3135* (RWorkspace: 0) 3136* 3137 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ), 3138 $ LDWRKR, VT, LDVT, CZERO, A, LDA ) 3139* 3140* Copy right singular vectors of A from A to VT 3141* 3142 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 3143* 3144 ELSE 3145* 3146* Insufficient workspace for a fast algorithm 3147* 3148 ITAU = 1 3149 IWORK = ITAU + M 3150* 3151* Compute A=L*Q, copying result to VT 3152* (CWorkspace: need 2*M, prefer M+M*NB) 3153* (RWorkspace: 0) 3154* 3155 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 3156 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3157 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3158* 3159* Generate Q in VT 3160* (CWorkspace: need M+N, prefer M+N*NB) 3161* (RWorkspace: 0) 3162* 3163 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 3164 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3165 IE = 1 3166 ITAUQ = ITAU 3167 ITAUP = ITAUQ + M 3168 IWORK = ITAUP + M 3169* 3170* Zero out above L in A 3171* 3172 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3173 $ A( 1, 2 ), LDA ) 3174* 3175* Bidiagonalize L in A 3176* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 3177* (RWorkspace: need M) 3178* 3179 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 3180 $ WORK( ITAUQ ), WORK( ITAUP ), 3181 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3182* 3183* Multiply right bidiagonalizing vectors in A by Q 3184* in VT 3185* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 3186* (RWorkspace: 0) 3187* 3188 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 3189 $ WORK( ITAUP ), VT, LDVT, 3190 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3191 IRWORK = IE + M 3192* 3193* Perform bidiagonal QR iteration, computing right 3194* singular vectors of A in VT 3195* (CWorkspace: 0) 3196* (RWorkspace: need BDSPAC) 3197* 3198 CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT, 3199 $ LDVT, CDUM, 1, CDUM, 1, 3200 $ RWORK( IRWORK ), INFO ) 3201* 3202 END IF 3203* 3204 ELSE IF( WNTUO ) THEN 3205* 3206* Path 8t(N much larger than M, JOBU='O', JOBVT='A') 3207* N right singular vectors to be computed in VT and 3208* M left singular vectors to be overwritten on A 3209* 3210 IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN 3211* 3212* Sufficient workspace for a fast algorithm 3213* 3214 IU = 1 3215 IF( LWORK.GE.WRKBL+2*LDA*M ) THEN 3216* 3217* WORK(IU) is LDA by M and WORK(IR) is LDA by M 3218* 3219 LDWRKU = LDA 3220 IR = IU + LDWRKU*M 3221 LDWRKR = LDA 3222 ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN 3223* 3224* WORK(IU) is LDA by M and WORK(IR) is M by M 3225* 3226 LDWRKU = LDA 3227 IR = IU + LDWRKU*M 3228 LDWRKR = M 3229 ELSE 3230* 3231* WORK(IU) is M by M and WORK(IR) is M by M 3232* 3233 LDWRKU = M 3234 IR = IU + LDWRKU*M 3235 LDWRKR = M 3236 END IF 3237 ITAU = IR + LDWRKR*M 3238 IWORK = ITAU + M 3239* 3240* Compute A=L*Q, copying result to VT 3241* (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 3242* (RWorkspace: 0) 3243* 3244 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 3245 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3246 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3247* 3248* Generate Q in VT 3249* (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) 3250* (RWorkspace: 0) 3251* 3252 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 3253 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3254* 3255* Copy L to WORK(IU), zeroing out above it 3256* 3257 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 3258 $ LDWRKU ) 3259 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3260 $ WORK( IU+LDWRKU ), LDWRKU ) 3261 IE = 1 3262 ITAUQ = ITAU 3263 ITAUP = ITAUQ + M 3264 IWORK = ITAUP + M 3265* 3266* Bidiagonalize L in WORK(IU), copying result to 3267* WORK(IR) 3268* (CWorkspace: need 2*M*M+3*M, 3269* prefer 2*M*M+2*M+2*M*NB) 3270* (RWorkspace: need M) 3271* 3272 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 3273 $ RWORK( IE ), WORK( ITAUQ ), 3274 $ WORK( ITAUP ), WORK( IWORK ), 3275 $ LWORK-IWORK+1, IERR ) 3276 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, 3277 $ WORK( IR ), LDWRKR ) 3278* 3279* Generate right bidiagonalizing vectors in WORK(IU) 3280* (CWorkspace: need 2*M*M+3*M-1, 3281* prefer 2*M*M+2*M+(M-1)*NB) 3282* (RWorkspace: 0) 3283* 3284 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 3285 $ WORK( ITAUP ), WORK( IWORK ), 3286 $ LWORK-IWORK+1, IERR ) 3287* 3288* Generate left bidiagonalizing vectors in WORK(IR) 3289* (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 3290* (RWorkspace: 0) 3291* 3292 CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, 3293 $ WORK( ITAUQ ), WORK( IWORK ), 3294 $ LWORK-IWORK+1, IERR ) 3295 IRWORK = IE + M 3296* 3297* Perform bidiagonal QR iteration, computing left 3298* singular vectors of L in WORK(IR) and computing 3299* right singular vectors of L in WORK(IU) 3300* (CWorkspace: need 2*M*M) 3301* (RWorkspace: need BDSPAC) 3302* 3303 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 3304 $ WORK( IU ), LDWRKU, WORK( IR ), 3305 $ LDWRKR, CDUM, 1, RWORK( IRWORK ), 3306 $ INFO ) 3307* 3308* Multiply right singular vectors of L in WORK(IU) by 3309* Q in VT, storing result in A 3310* (CWorkspace: need M*M) 3311* (RWorkspace: 0) 3312* 3313 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 3314 $ LDWRKU, VT, LDVT, CZERO, A, LDA ) 3315* 3316* Copy right singular vectors of A from A to VT 3317* 3318 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 3319* 3320* Copy left singular vectors of A from WORK(IR) to A 3321* 3322 CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A, 3323 $ LDA ) 3324* 3325 ELSE 3326* 3327* Insufficient workspace for a fast algorithm 3328* 3329 ITAU = 1 3330 IWORK = ITAU + M 3331* 3332* Compute A=L*Q, copying result to VT 3333* (CWorkspace: need 2*M, prefer M+M*NB) 3334* (RWorkspace: 0) 3335* 3336 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 3337 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3338 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3339* 3340* Generate Q in VT 3341* (CWorkspace: need M+N, prefer M+N*NB) 3342* (RWorkspace: 0) 3343* 3344 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 3345 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3346 IE = 1 3347 ITAUQ = ITAU 3348 ITAUP = ITAUQ + M 3349 IWORK = ITAUP + M 3350* 3351* Zero out above L in A 3352* 3353 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3354 $ A( 1, 2 ), LDA ) 3355* 3356* Bidiagonalize L in A 3357* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 3358* (RWorkspace: need M) 3359* 3360 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 3361 $ WORK( ITAUQ ), WORK( ITAUP ), 3362 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3363* 3364* Multiply right bidiagonalizing vectors in A by Q 3365* in VT 3366* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 3367* (RWorkspace: 0) 3368* 3369 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 3370 $ WORK( ITAUP ), VT, LDVT, 3371 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3372* 3373* Generate left bidiagonalizing vectors in A 3374* (CWorkspace: need 3*M, prefer 2*M+M*NB) 3375* (RWorkspace: 0) 3376* 3377 CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 3378 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3379 IRWORK = IE + M 3380* 3381* Perform bidiagonal QR iteration, computing left 3382* singular vectors of A in A and computing right 3383* singular vectors of A in VT 3384* (CWorkspace: 0) 3385* (RWorkspace: need BDSPAC) 3386* 3387 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 3388 $ LDVT, A, LDA, CDUM, 1, 3389 $ RWORK( IRWORK ), INFO ) 3390* 3391 END IF 3392* 3393 ELSE IF( WNTUAS ) THEN 3394* 3395* Path 9t(N much larger than M, JOBU='S' or 'A', 3396* JOBVT='A') 3397* N right singular vectors to be computed in VT and 3398* M left singular vectors to be computed in U 3399* 3400 IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN 3401* 3402* Sufficient workspace for a fast algorithm 3403* 3404 IU = 1 3405 IF( LWORK.GE.WRKBL+LDA*M ) THEN 3406* 3407* WORK(IU) is LDA by M 3408* 3409 LDWRKU = LDA 3410 ELSE 3411* 3412* WORK(IU) is M by M 3413* 3414 LDWRKU = M 3415 END IF 3416 ITAU = IU + LDWRKU*M 3417 IWORK = ITAU + M 3418* 3419* Compute A=L*Q, copying result to VT 3420* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 3421* (RWorkspace: 0) 3422* 3423 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 3424 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3425 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3426* 3427* Generate Q in VT 3428* (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) 3429* (RWorkspace: 0) 3430* 3431 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 3432 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3433* 3434* Copy L to WORK(IU), zeroing out above it 3435* 3436 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 3437 $ LDWRKU ) 3438 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3439 $ WORK( IU+LDWRKU ), LDWRKU ) 3440 IE = 1 3441 ITAUQ = ITAU 3442 ITAUP = ITAUQ + M 3443 IWORK = ITAUP + M 3444* 3445* Bidiagonalize L in WORK(IU), copying result to U 3446* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 3447* (RWorkspace: need M) 3448* 3449 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 3450 $ RWORK( IE ), WORK( ITAUQ ), 3451 $ WORK( ITAUP ), WORK( IWORK ), 3452 $ LWORK-IWORK+1, IERR ) 3453 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U, 3454 $ LDU ) 3455* 3456* Generate right bidiagonalizing vectors in WORK(IU) 3457* (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) 3458* (RWorkspace: 0) 3459* 3460 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 3461 $ WORK( ITAUP ), WORK( IWORK ), 3462 $ LWORK-IWORK+1, IERR ) 3463* 3464* Generate left bidiagonalizing vectors in U 3465* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 3466* (RWorkspace: 0) 3467* 3468 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 3469 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3470 IRWORK = IE + M 3471* 3472* Perform bidiagonal QR iteration, computing left 3473* singular vectors of L in U and computing right 3474* singular vectors of L in WORK(IU) 3475* (CWorkspace: need M*M) 3476* (RWorkspace: need BDSPAC) 3477* 3478 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 3479 $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1, 3480 $ RWORK( IRWORK ), INFO ) 3481* 3482* Multiply right singular vectors of L in WORK(IU) by 3483* Q in VT, storing result in A 3484* (CWorkspace: need M*M) 3485* (RWorkspace: 0) 3486* 3487 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 3488 $ LDWRKU, VT, LDVT, CZERO, A, LDA ) 3489* 3490* Copy right singular vectors of A from A to VT 3491* 3492 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 3493* 3494 ELSE 3495* 3496* Insufficient workspace for a fast algorithm 3497* 3498 ITAU = 1 3499 IWORK = ITAU + M 3500* 3501* Compute A=L*Q, copying result to VT 3502* (CWorkspace: need 2*M, prefer M+M*NB) 3503* (RWorkspace: 0) 3504* 3505 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 3506 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3507 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3508* 3509* Generate Q in VT 3510* (CWorkspace: need M+N, prefer M+N*NB) 3511* (RWorkspace: 0) 3512* 3513 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 3514 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3515* 3516* Copy L to U, zeroing out above it 3517* 3518 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 3519 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 3520 $ U( 1, 2 ), LDU ) 3521 IE = 1 3522 ITAUQ = ITAU 3523 ITAUP = ITAUQ + M 3524 IWORK = ITAUP + M 3525* 3526* Bidiagonalize L in U 3527* (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 3528* (RWorkspace: need M) 3529* 3530 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 3531 $ WORK( ITAUQ ), WORK( ITAUP ), 3532 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3533* 3534* Multiply right bidiagonalizing vectors in U by Q 3535* in VT 3536* (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 3537* (RWorkspace: 0) 3538* 3539 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 3540 $ WORK( ITAUP ), VT, LDVT, 3541 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3542* 3543* Generate left bidiagonalizing vectors in U 3544* (CWorkspace: need 3*M, prefer 2*M+M*NB) 3545* (RWorkspace: 0) 3546* 3547 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 3548 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3549 IRWORK = IE + M 3550* 3551* Perform bidiagonal QR iteration, computing left 3552* singular vectors of A in U and computing right 3553* singular vectors of A in VT 3554* (CWorkspace: 0) 3555* (RWorkspace: need BDSPAC) 3556* 3557 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 3558 $ LDVT, U, LDU, CDUM, 1, 3559 $ RWORK( IRWORK ), INFO ) 3560* 3561 END IF 3562* 3563 END IF 3564* 3565 END IF 3566* 3567 ELSE 3568* 3569* N .LT. MNTHR 3570* 3571* Path 10t(N greater than M, but not much larger) 3572* Reduce to bidiagonal form without LQ decomposition 3573* 3574 IE = 1 3575 ITAUQ = 1 3576 ITAUP = ITAUQ + M 3577 IWORK = ITAUP + M 3578* 3579* Bidiagonalize A 3580* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 3581* (RWorkspace: M) 3582* 3583 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 3584 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 3585 $ IERR ) 3586 IF( WNTUAS ) THEN 3587* 3588* If left singular vectors desired in U, copy result to U 3589* and generate left bidiagonalizing vectors in U 3590* (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) 3591* (RWorkspace: 0) 3592* 3593 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 3594 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 3595 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3596 END IF 3597 IF( WNTVAS ) THEN 3598* 3599* If right singular vectors desired in VT, copy result to 3600* VT and generate right bidiagonalizing vectors in VT 3601* (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB) 3602* (RWorkspace: 0) 3603* 3604 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 3605 IF( WNTVA ) 3606 $ NRVT = N 3607 IF( WNTVS ) 3608 $ NRVT = M 3609 CALL ZUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ), 3610 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3611 END IF 3612 IF( WNTUO ) THEN 3613* 3614* If left singular vectors desired in A, generate left 3615* bidiagonalizing vectors in A 3616* (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) 3617* (RWorkspace: 0) 3618* 3619 CALL ZUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), 3620 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3621 END IF 3622 IF( WNTVO ) THEN 3623* 3624* If right singular vectors desired in A, generate right 3625* bidiagonalizing vectors in A 3626* (CWorkspace: need 3*M, prefer 2*M+M*NB) 3627* (RWorkspace: 0) 3628* 3629 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 3630 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 3631 END IF 3632 IRWORK = IE + M 3633 IF( WNTUAS .OR. WNTUO ) 3634 $ NRU = M 3635 IF( WNTUN ) 3636 $ NRU = 0 3637 IF( WNTVAS .OR. WNTVO ) 3638 $ NCVT = N 3639 IF( WNTVN ) 3640 $ NCVT = 0 3641 IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN 3642* 3643* Perform bidiagonal QR iteration, if desired, computing 3644* left singular vectors in U and computing right singular 3645* vectors in VT 3646* (CWorkspace: 0) 3647* (RWorkspace: need BDSPAC) 3648* 3649 CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT, 3650 $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ), 3651 $ INFO ) 3652 ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN 3653* 3654* Perform bidiagonal QR iteration, if desired, computing 3655* left singular vectors in U and computing right singular 3656* vectors in A 3657* (CWorkspace: 0) 3658* (RWorkspace: need BDSPAC) 3659* 3660 CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A, 3661 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 3662 $ INFO ) 3663 ELSE 3664* 3665* Perform bidiagonal QR iteration, if desired, computing 3666* left singular vectors in A and computing right singular 3667* vectors in VT 3668* (CWorkspace: 0) 3669* (RWorkspace: need BDSPAC) 3670* 3671 CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT, 3672 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 3673 $ INFO ) 3674 END IF 3675* 3676 END IF 3677* 3678 END IF 3679* 3680* Undo scaling if necessary 3681* 3682 IF( ISCL.EQ.1 ) THEN 3683 IF( ANRM.GT.BIGNUM ) 3684 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 3685 $ IERR ) 3686 IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM ) 3687 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, 3688 $ RWORK( IE ), MINMN, IERR ) 3689 IF( ANRM.LT.SMLNUM ) 3690 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 3691 $ IERR ) 3692 IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM ) 3693 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, 3694 $ RWORK( IE ), MINMN, IERR ) 3695 END IF 3696* 3697* Return optimal workspace in WORK(1) 3698* 3699 WORK( 1 ) = MAXWRK 3700* 3701 RETURN 3702* 3703* End of ZGESVD 3704* 3705 END 3706