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