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