1*> \brief \b DGESDD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DGESDD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesdd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesdd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesdd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, 22* WORK, LWORK, 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 A( LDA, * ), S( * ), U( LDU, * ), 31* $ VT( LDVT, * ), WORK( * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> DGESDD computes the singular value decomposition (SVD) of a real 41*> M-by-N matrix A, optionally computing the left and right singular 42*> vectors. If singular vectors are desired, it uses a 43*> divide-and-conquer algorithm. 44*> 45*> The SVD is written 46*> 47*> A = U * SIGMA * transpose(V) 48*> 49*> where SIGMA is an M-by-N matrix which is zero except for its 50*> min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and 51*> V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA 52*> are the singular values of A; they are real and non-negative, and 53*> are returned in descending order. The first min(m,n) columns of 54*> U and V are the left and right singular vectors of A. 55*> 56*> Note that the routine returns VT = V**T, not V. 57*> 58*> The divide and conquer algorithm makes very mild assumptions about 59*> floating point arithmetic. It will work on machines with a guard 60*> digit in add/subtract, or on those binary machines without guard 61*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 62*> Cray-2. It could conceivably fail on hexadecimal or decimal machines 63*> without guard digits, but we know of none. 64*> \endverbatim 65* 66* Arguments: 67* ========== 68* 69*> \param[in] JOBZ 70*> \verbatim 71*> JOBZ is CHARACTER*1 72*> Specifies options for computing all or part of the matrix U: 73*> = 'A': all M columns of U and all N rows of V**T are 74*> returned in the arrays U and VT; 75*> = 'S': the first min(M,N) columns of U and the first 76*> min(M,N) rows of V**T are returned in the arrays U 77*> and VT; 78*> = 'O': If M >= N, the first N columns of U are overwritten 79*> on the array A and all rows of V**T are returned in 80*> the array VT; 81*> otherwise, all columns of U are returned in the 82*> array U and the first M rows of V**T are overwritten 83*> in the array A; 84*> = 'N': no columns of U or rows of V**T are computed. 85*> \endverbatim 86*> 87*> \param[in] M 88*> \verbatim 89*> M is INTEGER 90*> The number of rows of the input matrix A. M >= 0. 91*> \endverbatim 92*> 93*> \param[in] N 94*> \verbatim 95*> N is INTEGER 96*> The number of columns of the input matrix A. N >= 0. 97*> \endverbatim 98*> 99*> \param[in,out] A 100*> \verbatim 101*> A is DOUBLE PRECISION array, dimension (LDA,N) 102*> On entry, the M-by-N matrix A. 103*> On exit, 104*> if JOBZ = 'O', A is overwritten with the first N columns 105*> of U (the left singular vectors, stored 106*> columnwise) if M >= N; 107*> A is overwritten with the first M rows 108*> of V**T (the right singular vectors, stored 109*> rowwise) otherwise. 110*> if JOBZ .ne. 'O', the contents of A are destroyed. 111*> \endverbatim 112*> 113*> \param[in] LDA 114*> \verbatim 115*> LDA is INTEGER 116*> The leading dimension of the array A. LDA >= max(1,M). 117*> \endverbatim 118*> 119*> \param[out] S 120*> \verbatim 121*> S is DOUBLE PRECISION array, dimension (min(M,N)) 122*> The singular values of A, sorted so that S(i) >= S(i+1). 123*> \endverbatim 124*> 125*> \param[out] U 126*> \verbatim 127*> U is DOUBLE PRECISION array, dimension (LDU,UCOL) 128*> UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; 129*> UCOL = min(M,N) if JOBZ = 'S'. 130*> If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M 131*> orthogonal matrix U; 132*> if JOBZ = 'S', U contains the first min(M,N) columns of U 133*> (the left singular vectors, stored columnwise); 134*> if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. 135*> \endverbatim 136*> 137*> \param[in] LDU 138*> \verbatim 139*> LDU is INTEGER 140*> The leading dimension of the array U. LDU >= 1; if 141*> JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. 142*> \endverbatim 143*> 144*> \param[out] VT 145*> \verbatim 146*> VT is DOUBLE PRECISION array, dimension (LDVT,N) 147*> If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the 148*> N-by-N orthogonal matrix V**T; 149*> if JOBZ = 'S', VT contains the first min(M,N) rows of 150*> V**T (the right singular vectors, stored rowwise); 151*> if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. 152*> \endverbatim 153*> 154*> \param[in] LDVT 155*> \verbatim 156*> LDVT is INTEGER 157*> The leading dimension of the array VT. LDVT >= 1; 158*> if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; 159*> if JOBZ = 'S', LDVT >= min(M,N). 160*> \endverbatim 161*> 162*> \param[out] WORK 163*> \verbatim 164*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 165*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK; 166*> \endverbatim 167*> 168*> \param[in] LWORK 169*> \verbatim 170*> LWORK is INTEGER 171*> The dimension of the array WORK. LWORK >= 1. 172*> If LWORK = -1, a workspace query is assumed. The optimal 173*> size for the WORK array is calculated and stored in WORK(1), 174*> and no other work except argument checking is performed. 175*> 176*> Let mx = max(M,N) and mn = min(M,N). 177*> If JOBZ = 'N', LWORK >= 3*mn + max( mx, 7*mn ). 178*> If JOBZ = 'O', LWORK >= 3*mn + max( mx, 5*mn*mn + 4*mn ). 179*> If JOBZ = 'S', LWORK >= 4*mn*mn + 7*mn. 180*> If JOBZ = 'A', LWORK >= 4*mn*mn + 6*mn + mx. 181*> These are not tight minimums in all cases; see comments inside code. 182*> For good performance, LWORK should generally be larger; 183*> a query is recommended. 184*> \endverbatim 185*> 186*> \param[out] IWORK 187*> \verbatim 188*> IWORK is INTEGER array, dimension (8*min(M,N)) 189*> \endverbatim 190*> 191*> \param[out] INFO 192*> \verbatim 193*> INFO is INTEGER 194*> < 0: if INFO = -i, the i-th argument had an illegal value. 195*> = -4: if A had a NAN entry. 196*> > 0: DBDSDC did not converge, updating process failed. 197*> = 0: successful exit. 198*> \endverbatim 199* 200* Authors: 201* ======== 202* 203*> \author Univ. of Tennessee 204*> \author Univ. of California Berkeley 205*> \author Univ. of Colorado Denver 206*> \author NAG Ltd. 207* 208*> \ingroup doubleGEsing 209* 210*> \par Contributors: 211* ================== 212*> 213*> Ming Gu and Huan Ren, Computer Science Division, University of 214*> California at Berkeley, USA 215*> 216* ===================================================================== 217 SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, 218 $ WORK, LWORK, IWORK, INFO ) 219 implicit none 220* 221* -- LAPACK driver routine -- 222* -- LAPACK is a software package provided by Univ. of Tennessee, -- 223* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 224* 225* .. Scalar Arguments .. 226 CHARACTER JOBZ 227 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 228* .. 229* .. Array Arguments .. 230 INTEGER IWORK( * ) 231 DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), 232 $ VT( LDVT, * ), WORK( * ) 233* .. 234* 235* ===================================================================== 236* 237* .. Parameters .. 238 DOUBLE PRECISION ZERO, ONE 239 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 240* .. 241* .. Local Scalars .. 242 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS 243 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL, 244 $ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, 245 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, 246 $ MNTHR, NWORK, WRKBL 247 INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM, 248 $ LWORK_DGEBRD_NN, LWORK_DGELQF_MN, 249 $ LWORK_DGEQRF_MN, 250 $ LWORK_DORGBR_P_MM, LWORK_DORGBR_Q_NN, 251 $ LWORK_DORGLQ_MN, LWORK_DORGLQ_NN, 252 $ LWORK_DORGQR_MM, LWORK_DORGQR_MN, 253 $ LWORK_DORMBR_PRT_MM, LWORK_DORMBR_QLN_MM, 254 $ LWORK_DORMBR_PRT_MN, LWORK_DORMBR_QLN_MN, 255 $ LWORK_DORMBR_PRT_NN, LWORK_DORMBR_QLN_NN 256 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM 257* .. 258* .. Local Arrays .. 259 INTEGER IDUM( 1 ) 260 DOUBLE PRECISION DUM( 1 ) 261* .. 262* .. External Subroutines .. 263 EXTERNAL DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY, 264 $ DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR, 265 $ XERBLA 266* .. 267* .. External Functions .. 268 LOGICAL LSAME, DISNAN 269 DOUBLE PRECISION DLAMCH, DLANGE 270 EXTERNAL DLAMCH, DLANGE, LSAME, DISNAN 271* .. 272* .. Intrinsic Functions .. 273 INTRINSIC INT, MAX, MIN, SQRT 274* .. 275* .. Executable Statements .. 276* 277* Test the input arguments 278* 279 INFO = 0 280 MINMN = MIN( M, N ) 281 WNTQA = LSAME( JOBZ, 'A' ) 282 WNTQS = LSAME( JOBZ, 'S' ) 283 WNTQAS = WNTQA .OR. WNTQS 284 WNTQO = LSAME( JOBZ, 'O' ) 285 WNTQN = LSAME( JOBZ, 'N' ) 286 LQUERY = ( LWORK.EQ.-1 ) 287* 288 IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN 289 INFO = -1 290 ELSE IF( M.LT.0 ) THEN 291 INFO = -2 292 ELSE IF( N.LT.0 ) THEN 293 INFO = -3 294 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 295 INFO = -5 296 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR. 297 $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN 298 INFO = -8 299 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR. 300 $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR. 301 $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN 302 INFO = -10 303 END IF 304* 305* Compute workspace 306* Note: Comments in the code beginning "Workspace:" describe the 307* minimal amount of workspace allocated at that point in the code, 308* as well as the preferred amount for good performance. 309* NB refers to the optimal block size for the immediately 310* following subroutine, as returned by ILAENV. 311* 312 IF( INFO.EQ.0 ) THEN 313 MINWRK = 1 314 MAXWRK = 1 315 BDSPAC = 0 316 MNTHR = INT( MINMN*11.0D0 / 6.0D0 ) 317 IF( M.GE.N .AND. MINMN.GT.0 ) THEN 318* 319* Compute space needed for DBDSDC 320* 321 IF( WNTQN ) THEN 322* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6) 323* keep 7*N for backwards compatibility. 324 BDSPAC = 7*N 325 ELSE 326 BDSPAC = 3*N*N + 4*N 327 END IF 328* 329* Compute space preferred for each routine 330 CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1), 331 $ DUM(1), DUM(1), -1, IERR ) 332 LWORK_DGEBRD_MN = INT( DUM(1) ) 333* 334 CALL DGEBRD( N, N, DUM(1), N, DUM(1), DUM(1), DUM(1), 335 $ DUM(1), DUM(1), -1, IERR ) 336 LWORK_DGEBRD_NN = INT( DUM(1) ) 337* 338 CALL DGEQRF( M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR ) 339 LWORK_DGEQRF_MN = INT( DUM(1) ) 340* 341 CALL DORGBR( 'Q', N, N, N, DUM(1), N, DUM(1), DUM(1), -1, 342 $ IERR ) 343 LWORK_DORGBR_Q_NN = INT( DUM(1) ) 344* 345 CALL DORGQR( M, M, N, DUM(1), M, DUM(1), DUM(1), -1, IERR ) 346 LWORK_DORGQR_MM = INT( DUM(1) ) 347* 348 CALL DORGQR( M, N, N, DUM(1), M, DUM(1), DUM(1), -1, IERR ) 349 LWORK_DORGQR_MN = INT( DUM(1) ) 350* 351 CALL DORMBR( 'P', 'R', 'T', N, N, N, DUM(1), N, 352 $ DUM(1), DUM(1), N, DUM(1), -1, IERR ) 353 LWORK_DORMBR_PRT_NN = INT( DUM(1) ) 354* 355 CALL DORMBR( 'Q', 'L', 'N', N, N, N, DUM(1), N, 356 $ DUM(1), DUM(1), N, DUM(1), -1, IERR ) 357 LWORK_DORMBR_QLN_NN = INT( DUM(1) ) 358* 359 CALL DORMBR( 'Q', 'L', 'N', M, N, N, DUM(1), M, 360 $ DUM(1), DUM(1), M, DUM(1), -1, IERR ) 361 LWORK_DORMBR_QLN_MN = INT( DUM(1) ) 362* 363 CALL DORMBR( 'Q', 'L', 'N', M, M, N, DUM(1), M, 364 $ DUM(1), DUM(1), M, DUM(1), -1, IERR ) 365 LWORK_DORMBR_QLN_MM = INT( DUM(1) ) 366* 367 IF( M.GE.MNTHR ) THEN 368 IF( WNTQN ) THEN 369* 370* Path 1 (M >> N, JOBZ='N') 371* 372 WRKBL = N + LWORK_DGEQRF_MN 373 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN ) 374 MAXWRK = MAX( WRKBL, BDSPAC + N ) 375 MINWRK = BDSPAC + N 376 ELSE IF( WNTQO ) THEN 377* 378* Path 2 (M >> N, JOBZ='O') 379* 380 WRKBL = N + LWORK_DGEQRF_MN 381 WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MN ) 382 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN ) 383 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN ) 384 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN ) 385 WRKBL = MAX( WRKBL, 3*N + BDSPAC ) 386 MAXWRK = WRKBL + 2*N*N 387 MINWRK = BDSPAC + 2*N*N + 3*N 388 ELSE IF( WNTQS ) THEN 389* 390* Path 3 (M >> N, JOBZ='S') 391* 392 WRKBL = N + LWORK_DGEQRF_MN 393 WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MN ) 394 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN ) 395 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN ) 396 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN ) 397 WRKBL = MAX( WRKBL, 3*N + BDSPAC ) 398 MAXWRK = WRKBL + N*N 399 MINWRK = BDSPAC + N*N + 3*N 400 ELSE IF( WNTQA ) THEN 401* 402* Path 4 (M >> N, JOBZ='A') 403* 404 WRKBL = N + LWORK_DGEQRF_MN 405 WRKBL = MAX( WRKBL, N + LWORK_DORGQR_MM ) 406 WRKBL = MAX( WRKBL, 3*N + LWORK_DGEBRD_NN ) 407 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_NN ) 408 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN ) 409 WRKBL = MAX( WRKBL, 3*N + BDSPAC ) 410 MAXWRK = WRKBL + N*N 411 MINWRK = N*N + MAX( 3*N + BDSPAC, N + M ) 412 END IF 413 ELSE 414* 415* Path 5 (M >= N, but not much larger) 416* 417 WRKBL = 3*N + LWORK_DGEBRD_MN 418 IF( WNTQN ) THEN 419* Path 5n (M >= N, jobz='N') 420 MAXWRK = MAX( WRKBL, 3*N + BDSPAC ) 421 MINWRK = 3*N + MAX( M, BDSPAC ) 422 ELSE IF( WNTQO ) THEN 423* Path 5o (M >= N, jobz='O') 424 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN ) 425 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN ) 426 WRKBL = MAX( WRKBL, 3*N + BDSPAC ) 427 MAXWRK = WRKBL + M*N 428 MINWRK = 3*N + MAX( M, N*N + BDSPAC ) 429 ELSE IF( WNTQS ) THEN 430* Path 5s (M >= N, jobz='S') 431 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MN ) 432 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN ) 433 MAXWRK = MAX( WRKBL, 3*N + BDSPAC ) 434 MINWRK = 3*N + MAX( M, BDSPAC ) 435 ELSE IF( WNTQA ) THEN 436* Path 5a (M >= N, jobz='A') 437 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_QLN_MM ) 438 WRKBL = MAX( WRKBL, 3*N + LWORK_DORMBR_PRT_NN ) 439 MAXWRK = MAX( WRKBL, 3*N + BDSPAC ) 440 MINWRK = 3*N + MAX( M, BDSPAC ) 441 END IF 442 END IF 443 ELSE IF( MINMN.GT.0 ) THEN 444* 445* Compute space needed for DBDSDC 446* 447 IF( WNTQN ) THEN 448* dbdsdc needs only 4*N (or 6*N for uplo=L for LAPACK <= 3.6) 449* keep 7*N for backwards compatibility. 450 BDSPAC = 7*M 451 ELSE 452 BDSPAC = 3*M*M + 4*M 453 END IF 454* 455* Compute space preferred for each routine 456 CALL DGEBRD( M, N, DUM(1), M, DUM(1), DUM(1), DUM(1), 457 $ DUM(1), DUM(1), -1, IERR ) 458 LWORK_DGEBRD_MN = INT( DUM(1) ) 459* 460 CALL DGEBRD( M, M, A, M, S, DUM(1), DUM(1), 461 $ DUM(1), DUM(1), -1, IERR ) 462 LWORK_DGEBRD_MM = INT( DUM(1) ) 463* 464 CALL DGELQF( M, N, A, M, DUM(1), DUM(1), -1, IERR ) 465 LWORK_DGELQF_MN = INT( DUM(1) ) 466* 467 CALL DORGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR ) 468 LWORK_DORGLQ_NN = INT( DUM(1) ) 469* 470 CALL DORGLQ( M, N, M, A, M, DUM(1), DUM(1), -1, IERR ) 471 LWORK_DORGLQ_MN = INT( DUM(1) ) 472* 473 CALL DORGBR( 'P', M, M, M, A, N, DUM(1), DUM(1), -1, IERR ) 474 LWORK_DORGBR_P_MM = INT( DUM(1) ) 475* 476 CALL DORMBR( 'P', 'R', 'T', M, M, M, DUM(1), M, 477 $ DUM(1), DUM(1), M, DUM(1), -1, IERR ) 478 LWORK_DORMBR_PRT_MM = INT( DUM(1) ) 479* 480 CALL DORMBR( 'P', 'R', 'T', M, N, M, DUM(1), M, 481 $ DUM(1), DUM(1), M, DUM(1), -1, IERR ) 482 LWORK_DORMBR_PRT_MN = INT( DUM(1) ) 483* 484 CALL DORMBR( 'P', 'R', 'T', N, N, M, DUM(1), N, 485 $ DUM(1), DUM(1), N, DUM(1), -1, IERR ) 486 LWORK_DORMBR_PRT_NN = INT( DUM(1) ) 487* 488 CALL DORMBR( 'Q', 'L', 'N', M, M, M, DUM(1), M, 489 $ DUM(1), DUM(1), M, DUM(1), -1, IERR ) 490 LWORK_DORMBR_QLN_MM = INT( DUM(1) ) 491* 492 IF( N.GE.MNTHR ) THEN 493 IF( WNTQN ) THEN 494* 495* Path 1t (N >> M, JOBZ='N') 496* 497 WRKBL = M + LWORK_DGELQF_MN 498 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM ) 499 MAXWRK = MAX( WRKBL, BDSPAC + M ) 500 MINWRK = BDSPAC + M 501 ELSE IF( WNTQO ) THEN 502* 503* Path 2t (N >> M, JOBZ='O') 504* 505 WRKBL = M + LWORK_DGELQF_MN 506 WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_MN ) 507 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM ) 508 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM ) 509 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM ) 510 WRKBL = MAX( WRKBL, 3*M + BDSPAC ) 511 MAXWRK = WRKBL + 2*M*M 512 MINWRK = BDSPAC + 2*M*M + 3*M 513 ELSE IF( WNTQS ) THEN 514* 515* Path 3t (N >> M, JOBZ='S') 516* 517 WRKBL = M + LWORK_DGELQF_MN 518 WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_MN ) 519 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM ) 520 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM ) 521 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM ) 522 WRKBL = MAX( WRKBL, 3*M + BDSPAC ) 523 MAXWRK = WRKBL + M*M 524 MINWRK = BDSPAC + M*M + 3*M 525 ELSE IF( WNTQA ) THEN 526* 527* Path 4t (N >> M, JOBZ='A') 528* 529 WRKBL = M + LWORK_DGELQF_MN 530 WRKBL = MAX( WRKBL, M + LWORK_DORGLQ_NN ) 531 WRKBL = MAX( WRKBL, 3*M + LWORK_DGEBRD_MM ) 532 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM ) 533 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MM ) 534 WRKBL = MAX( WRKBL, 3*M + BDSPAC ) 535 MAXWRK = WRKBL + M*M 536 MINWRK = M*M + MAX( 3*M + BDSPAC, M + N ) 537 END IF 538 ELSE 539* 540* Path 5t (N > M, but not much larger) 541* 542 WRKBL = 3*M + LWORK_DGEBRD_MN 543 IF( WNTQN ) THEN 544* Path 5tn (N > M, jobz='N') 545 MAXWRK = MAX( WRKBL, 3*M + BDSPAC ) 546 MINWRK = 3*M + MAX( N, BDSPAC ) 547 ELSE IF( WNTQO ) THEN 548* Path 5to (N > M, jobz='O') 549 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM ) 550 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN ) 551 WRKBL = MAX( WRKBL, 3*M + BDSPAC ) 552 MAXWRK = WRKBL + M*N 553 MINWRK = 3*M + MAX( N, M*M + BDSPAC ) 554 ELSE IF( WNTQS ) THEN 555* Path 5ts (N > M, jobz='S') 556 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM ) 557 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_MN ) 558 MAXWRK = MAX( WRKBL, 3*M + BDSPAC ) 559 MINWRK = 3*M + MAX( N, BDSPAC ) 560 ELSE IF( WNTQA ) THEN 561* Path 5ta (N > M, jobz='A') 562 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_QLN_MM ) 563 WRKBL = MAX( WRKBL, 3*M + LWORK_DORMBR_PRT_NN ) 564 MAXWRK = MAX( WRKBL, 3*M + BDSPAC ) 565 MINWRK = 3*M + MAX( N, BDSPAC ) 566 END IF 567 END IF 568 END IF 569 570 MAXWRK = MAX( MAXWRK, MINWRK ) 571 WORK( 1 ) = MAXWRK 572* 573 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 574 INFO = -12 575 END IF 576 END IF 577* 578 IF( INFO.NE.0 ) THEN 579 CALL XERBLA( 'DGESDD', -INFO ) 580 RETURN 581 ELSE IF( LQUERY ) THEN 582 RETURN 583 END IF 584* 585* Quick return if possible 586* 587 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 588 RETURN 589 END IF 590* 591* Get machine constants 592* 593 EPS = DLAMCH( 'P' ) 594 SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS 595 BIGNUM = ONE / SMLNUM 596* 597* Scale A if max element outside range [SMLNUM,BIGNUM] 598* 599 ANRM = DLANGE( 'M', M, N, A, LDA, DUM ) 600 IF( DISNAN( ANRM ) ) THEN 601 INFO = -4 602 RETURN 603 END IF 604 ISCL = 0 605 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 606 ISCL = 1 607 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 608 ELSE IF( ANRM.GT.BIGNUM ) THEN 609 ISCL = 1 610 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 611 END IF 612* 613 IF( M.GE.N ) THEN 614* 615* A has at least as many rows as columns. If A has sufficiently 616* more rows than columns, first reduce using the QR 617* decomposition (if sufficient workspace available) 618* 619 IF( M.GE.MNTHR ) THEN 620* 621 IF( WNTQN ) THEN 622* 623* Path 1 (M >> N, JOBZ='N') 624* No singular vectors to be computed 625* 626 ITAU = 1 627 NWORK = ITAU + N 628* 629* Compute A=Q*R 630* Workspace: need N [tau] + N [work] 631* Workspace: prefer N [tau] + N*NB [work] 632* 633 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 634 $ LWORK - NWORK + 1, IERR ) 635* 636* Zero out below R 637* 638 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 639 IE = 1 640 ITAUQ = IE + N 641 ITAUP = ITAUQ + N 642 NWORK = ITAUP + N 643* 644* Bidiagonalize R in A 645* Workspace: need 3*N [e, tauq, taup] + N [work] 646* Workspace: prefer 3*N [e, tauq, taup] + 2*N*NB [work] 647* 648 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 649 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 650 $ IERR ) 651 NWORK = IE + N 652* 653* Perform bidiagonal SVD, computing singular values only 654* Workspace: need N [e] + BDSPAC 655* 656 CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1, 657 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 658* 659 ELSE IF( WNTQO ) THEN 660* 661* Path 2 (M >> N, JOBZ = 'O') 662* N left singular vectors to be overwritten on A and 663* N right singular vectors to be computed in VT 664* 665 IR = 1 666* 667* WORK(IR) is LDWRKR by N 668* 669 IF( LWORK .GE. LDA*N + N*N + 3*N + BDSPAC ) THEN 670 LDWRKR = LDA 671 ELSE 672 LDWRKR = ( LWORK - N*N - 3*N - BDSPAC ) / N 673 END IF 674 ITAU = IR + LDWRKR*N 675 NWORK = ITAU + N 676* 677* Compute A=Q*R 678* Workspace: need N*N [R] + N [tau] + N [work] 679* Workspace: prefer N*N [R] + N [tau] + N*NB [work] 680* 681 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 682 $ LWORK - NWORK + 1, IERR ) 683* 684* Copy R to WORK(IR), zeroing out below it 685* 686 CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 687 CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1), 688 $ LDWRKR ) 689* 690* Generate Q in A 691* Workspace: need N*N [R] + N [tau] + N [work] 692* Workspace: prefer N*N [R] + N [tau] + N*NB [work] 693* 694 CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), 695 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 696 IE = ITAU 697 ITAUQ = IE + N 698 ITAUP = ITAUQ + N 699 NWORK = ITAUP + N 700* 701* Bidiagonalize R in WORK(IR) 702* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work] 703* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work] 704* 705 CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), 706 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 707 $ LWORK - NWORK + 1, IERR ) 708* 709* WORK(IU) is N by N 710* 711 IU = NWORK 712 NWORK = IU + N*N 713* 714* Perform bidiagonal SVD, computing left singular vectors 715* of bidiagonal matrix in WORK(IU) and computing right 716* singular vectors of bidiagonal matrix in VT 717* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + BDSPAC 718* 719 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N, 720 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 721 $ INFO ) 722* 723* Overwrite WORK(IU) by left singular vectors of R 724* and VT by right singular vectors of R 725* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N [work] 726* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*N [U] + N*NB [work] 727* 728 CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 729 $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ), 730 $ LWORK - NWORK + 1, IERR ) 731 CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR, 732 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 733 $ LWORK - NWORK + 1, IERR ) 734* 735* Multiply Q in A by left singular vectors of R in 736* WORK(IU), storing result in WORK(IR) and copying to A 737* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N*N [U] 738* Workspace: prefer M*N [R] + 3*N [e, tauq, taup] + N*N [U] 739* 740 DO 10 I = 1, M, LDWRKR 741 CHUNK = MIN( M - I + 1, LDWRKR ) 742 CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ), 743 $ LDA, WORK( IU ), N, ZERO, WORK( IR ), 744 $ LDWRKR ) 745 CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 746 $ A( I, 1 ), LDA ) 747 10 CONTINUE 748* 749 ELSE IF( WNTQS ) THEN 750* 751* Path 3 (M >> N, JOBZ='S') 752* N left singular vectors to be computed in U and 753* N right singular vectors to be computed in VT 754* 755 IR = 1 756* 757* WORK(IR) is N by N 758* 759 LDWRKR = N 760 ITAU = IR + LDWRKR*N 761 NWORK = ITAU + N 762* 763* Compute A=Q*R 764* Workspace: need N*N [R] + N [tau] + N [work] 765* Workspace: prefer N*N [R] + N [tau] + N*NB [work] 766* 767 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 768 $ LWORK - NWORK + 1, IERR ) 769* 770* Copy R to WORK(IR), zeroing out below it 771* 772 CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 773 CALL DLASET( 'L', N - 1, N - 1, ZERO, ZERO, WORK(IR+1), 774 $ LDWRKR ) 775* 776* Generate Q in A 777* Workspace: need N*N [R] + N [tau] + N [work] 778* Workspace: prefer N*N [R] + N [tau] + N*NB [work] 779* 780 CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ), 781 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 782 IE = ITAU 783 ITAUQ = IE + N 784 ITAUP = ITAUQ + N 785 NWORK = ITAUP + N 786* 787* Bidiagonalize R in WORK(IR) 788* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work] 789* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + 2*N*NB [work] 790* 791 CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ), 792 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 793 $ LWORK - NWORK + 1, IERR ) 794* 795* Perform bidiagonal SVD, computing left singular vectors 796* of bidiagoal matrix in U and computing right singular 797* vectors of bidiagonal matrix in VT 798* Workspace: need N*N [R] + 3*N [e, tauq, taup] + BDSPAC 799* 800 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, 801 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 802 $ INFO ) 803* 804* Overwrite U by left singular vectors of R and VT 805* by right singular vectors of R 806* Workspace: need N*N [R] + 3*N [e, tauq, taup] + N [work] 807* Workspace: prefer N*N [R] + 3*N [e, tauq, taup] + N*NB [work] 808* 809 CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 810 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 811 $ LWORK - NWORK + 1, IERR ) 812* 813 CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR, 814 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 815 $ LWORK - NWORK + 1, IERR ) 816* 817* Multiply Q in A by left singular vectors of R in 818* WORK(IR), storing result in U 819* Workspace: need N*N [R] 820* 821 CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) 822 CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ), 823 $ LDWRKR, ZERO, U, LDU ) 824* 825 ELSE IF( WNTQA ) THEN 826* 827* Path 4 (M >> N, JOBZ='A') 828* M left singular vectors to be computed in U and 829* N right singular vectors to be computed in VT 830* 831 IU = 1 832* 833* WORK(IU) is N by N 834* 835 LDWRKU = N 836 ITAU = IU + LDWRKU*N 837 NWORK = ITAU + N 838* 839* Compute A=Q*R, copying result to U 840* Workspace: need N*N [U] + N [tau] + N [work] 841* Workspace: prefer N*N [U] + N [tau] + N*NB [work] 842* 843 CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 844 $ LWORK - NWORK + 1, IERR ) 845 CALL DLACPY( 'L', M, N, A, LDA, U, LDU ) 846* 847* Generate Q in U 848* Workspace: need N*N [U] + N [tau] + M [work] 849* Workspace: prefer N*N [U] + N [tau] + M*NB [work] 850 CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ), 851 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 852* 853* Produce R in A, zeroing out other entries 854* 855 CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA ) 856 IE = ITAU 857 ITAUQ = IE + N 858 ITAUP = ITAUQ + N 859 NWORK = ITAUP + N 860* 861* Bidiagonalize R in A 862* Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work] 863* Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + 2*N*NB [work] 864* 865 CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 866 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 867 $ IERR ) 868* 869* Perform bidiagonal SVD, computing left singular vectors 870* of bidiagonal matrix in WORK(IU) and computing right 871* singular vectors of bidiagonal matrix in VT 872* Workspace: need N*N [U] + 3*N [e, tauq, taup] + BDSPAC 873* 874 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N, 875 $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 876 $ INFO ) 877* 878* Overwrite WORK(IU) by left singular vectors of R and VT 879* by right singular vectors of R 880* Workspace: need N*N [U] + 3*N [e, tauq, taup] + N [work] 881* Workspace: prefer N*N [U] + 3*N [e, tauq, taup] + N*NB [work] 882* 883 CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA, 884 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 885 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 886 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, 887 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 888 $ LWORK - NWORK + 1, IERR ) 889* 890* Multiply Q in U by left singular vectors of R in 891* WORK(IU), storing result in A 892* Workspace: need N*N [U] 893* 894 CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ), 895 $ LDWRKU, ZERO, A, LDA ) 896* 897* Copy left singular vectors of A from A to U 898* 899 CALL DLACPY( 'F', M, N, A, LDA, U, LDU ) 900* 901 END IF 902* 903 ELSE 904* 905* M .LT. MNTHR 906* 907* Path 5 (M >= N, but not much larger) 908* Reduce to bidiagonal form without QR decomposition 909* 910 IE = 1 911 ITAUQ = IE + N 912 ITAUP = ITAUQ + N 913 NWORK = ITAUP + N 914* 915* Bidiagonalize A 916* Workspace: need 3*N [e, tauq, taup] + M [work] 917* Workspace: prefer 3*N [e, tauq, taup] + (M+N)*NB [work] 918* 919 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 920 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 921 $ IERR ) 922 IF( WNTQN ) THEN 923* 924* Path 5n (M >= N, JOBZ='N') 925* Perform bidiagonal SVD, only computing singular values 926* Workspace: need 3*N [e, tauq, taup] + BDSPAC 927* 928 CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1, 929 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 930 ELSE IF( WNTQO ) THEN 931* Path 5o (M >= N, JOBZ='O') 932 IU = NWORK 933 IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN 934* 935* WORK( IU ) is M by N 936* 937 LDWRKU = M 938 NWORK = IU + LDWRKU*N 939 CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ), 940 $ LDWRKU ) 941* IR is unused; silence compile warnings 942 IR = -1 943 ELSE 944* 945* WORK( IU ) is N by N 946* 947 LDWRKU = N 948 NWORK = IU + LDWRKU*N 949* 950* WORK(IR) is LDWRKR by N 951* 952 IR = NWORK 953 LDWRKR = ( LWORK - N*N - 3*N ) / N 954 END IF 955 NWORK = IU + LDWRKU*N 956* 957* Perform bidiagonal SVD, computing left singular vectors 958* of bidiagonal matrix in WORK(IU) and computing right 959* singular vectors of bidiagonal matrix in VT 960* Workspace: need 3*N [e, tauq, taup] + N*N [U] + BDSPAC 961* 962 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), 963 $ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ), 964 $ IWORK, INFO ) 965* 966* Overwrite VT by right singular vectors of A 967* Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work] 968* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work] 969* 970 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, 971 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 972 $ LWORK - NWORK + 1, IERR ) 973* 974 IF( LWORK .GE. M*N + 3*N + BDSPAC ) THEN 975* 976* Path 5o-fast 977* Overwrite WORK(IU) by left singular vectors of A 978* Workspace: need 3*N [e, tauq, taup] + M*N [U] + N [work] 979* Workspace: prefer 3*N [e, tauq, taup] + M*N [U] + N*NB [work] 980* 981 CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 982 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 983 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 984* 985* Copy left singular vectors of A from WORK(IU) to A 986* 987 CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) 988 ELSE 989* 990* Path 5o-slow 991* Generate Q in A 992* Workspace: need 3*N [e, tauq, taup] + N*N [U] + N [work] 993* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + N*NB [work] 994* 995 CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 996 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 997* 998* Multiply Q in A by left singular vectors of 999* bidiagonal matrix in WORK(IU), storing result in 1000* WORK(IR) and copying to A 1001* Workspace: need 3*N [e, tauq, taup] + N*N [U] + NB*N [R] 1002* Workspace: prefer 3*N [e, tauq, taup] + N*N [U] + M*N [R] 1003* 1004 DO 20 I = 1, M, LDWRKR 1005 CHUNK = MIN( M - I + 1, LDWRKR ) 1006 CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ), 1007 $ LDA, WORK( IU ), LDWRKU, ZERO, 1008 $ WORK( IR ), LDWRKR ) 1009 CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 1010 $ A( I, 1 ), LDA ) 1011 20 CONTINUE 1012 END IF 1013* 1014 ELSE IF( WNTQS ) THEN 1015* 1016* Path 5s (M >= N, JOBZ='S') 1017* Perform bidiagonal SVD, computing left singular vectors 1018* of bidiagonal matrix in U and computing right singular 1019* vectors of bidiagonal matrix in VT 1020* Workspace: need 3*N [e, tauq, taup] + BDSPAC 1021* 1022 CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU ) 1023 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, 1024 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 1025 $ INFO ) 1026* 1027* Overwrite U by left singular vectors of A and VT 1028* by right singular vectors of A 1029* Workspace: need 3*N [e, tauq, taup] + N [work] 1030* Workspace: prefer 3*N [e, tauq, taup] + N*NB [work] 1031* 1032 CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 1033 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1034 $ LWORK - NWORK + 1, IERR ) 1035 CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA, 1036 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1037 $ LWORK - NWORK + 1, IERR ) 1038 ELSE IF( WNTQA ) THEN 1039* 1040* Path 5a (M >= N, JOBZ='A') 1041* Perform bidiagonal SVD, computing left singular vectors 1042* of bidiagonal matrix in U and computing right singular 1043* vectors of bidiagonal matrix in VT 1044* Workspace: need 3*N [e, tauq, taup] + BDSPAC 1045* 1046 CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU ) 1047 CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT, 1048 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 1049 $ INFO ) 1050* 1051* Set the right corner of U to identity matrix 1052* 1053 IF( M.GT.N ) THEN 1054 CALL DLASET( 'F', M - N, M - N, ZERO, ONE, U(N+1,N+1), 1055 $ LDU ) 1056 END IF 1057* 1058* Overwrite U by left singular vectors of A and VT 1059* by right singular vectors of A 1060* Workspace: need 3*N [e, tauq, taup] + M [work] 1061* Workspace: prefer 3*N [e, tauq, taup] + M*NB [work] 1062* 1063 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1064 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1065 $ LWORK - NWORK + 1, IERR ) 1066 CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA, 1067 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1068 $ LWORK - NWORK + 1, IERR ) 1069 END IF 1070* 1071 END IF 1072* 1073 ELSE 1074* 1075* A has more columns than rows. If A has sufficiently more 1076* columns than rows, first reduce using the LQ decomposition (if 1077* sufficient workspace available) 1078* 1079 IF( N.GE.MNTHR ) THEN 1080* 1081 IF( WNTQN ) THEN 1082* 1083* Path 1t (N >> M, JOBZ='N') 1084* No singular vectors to be computed 1085* 1086 ITAU = 1 1087 NWORK = ITAU + M 1088* 1089* Compute A=L*Q 1090* Workspace: need M [tau] + M [work] 1091* Workspace: prefer M [tau] + M*NB [work] 1092* 1093 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1094 $ LWORK - NWORK + 1, IERR ) 1095* 1096* Zero out above L 1097* 1098 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA ) 1099 IE = 1 1100 ITAUQ = IE + M 1101 ITAUP = ITAUQ + M 1102 NWORK = ITAUP + M 1103* 1104* Bidiagonalize L in A 1105* Workspace: need 3*M [e, tauq, taup] + M [work] 1106* Workspace: prefer 3*M [e, tauq, taup] + 2*M*NB [work] 1107* 1108 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 1109 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1110 $ IERR ) 1111 NWORK = IE + M 1112* 1113* Perform bidiagonal SVD, computing singular values only 1114* Workspace: need M [e] + BDSPAC 1115* 1116 CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1, 1117 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 1118* 1119 ELSE IF( WNTQO ) THEN 1120* 1121* Path 2t (N >> M, JOBZ='O') 1122* M right singular vectors to be overwritten on A and 1123* M left singular vectors to be computed in U 1124* 1125 IVT = 1 1126* 1127* WORK(IVT) is M by M 1128* WORK(IL) is M by M; it is later resized to M by chunk for gemm 1129* 1130 IL = IVT + M*M 1131 IF( LWORK .GE. M*N + M*M + 3*M + BDSPAC ) THEN 1132 LDWRKL = M 1133 CHUNK = N 1134 ELSE 1135 LDWRKL = M 1136 CHUNK = ( LWORK - M*M ) / M 1137 END IF 1138 ITAU = IL + LDWRKL*M 1139 NWORK = ITAU + M 1140* 1141* Compute A=L*Q 1142* Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work] 1143* Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] 1144* 1145 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1146 $ LWORK - NWORK + 1, IERR ) 1147* 1148* Copy L to WORK(IL), zeroing about above it 1149* 1150 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 1151 CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO, 1152 $ WORK( IL + LDWRKL ), LDWRKL ) 1153* 1154* Generate Q in A 1155* Workspace: need M*M [VT] + M*M [L] + M [tau] + M [work] 1156* Workspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work] 1157* 1158 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), 1159 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1160 IE = ITAU 1161 ITAUQ = IE + M 1162 ITAUP = ITAUQ + M 1163 NWORK = ITAUP + M 1164* 1165* Bidiagonalize L in WORK(IL) 1166* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work] 1167* Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work] 1168* 1169 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ), 1170 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 1171 $ LWORK - NWORK + 1, IERR ) 1172* 1173* Perform bidiagonal SVD, computing left singular vectors 1174* of bidiagonal matrix in U, and computing right singular 1175* vectors of bidiagonal matrix in WORK(IVT) 1176* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + BDSPAC 1177* 1178 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, 1179 $ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ), 1180 $ IWORK, INFO ) 1181* 1182* Overwrite U by left singular vectors of L and WORK(IVT) 1183* by right singular vectors of L 1184* Workspace: need M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M [work] 1185* Workspace: prefer M*M [VT] + M*M [L] + 3*M [e, tauq, taup] + M*NB [work] 1186* 1187 CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 1188 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1189 $ LWORK - NWORK + 1, IERR ) 1190 CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL, 1191 $ WORK( ITAUP ), WORK( IVT ), M, 1192 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1193* 1194* Multiply right singular vectors of L in WORK(IVT) by Q 1195* in A, storing result in WORK(IL) and copying to A 1196* Workspace: need M*M [VT] + M*M [L] 1197* Workspace: prefer M*M [VT] + M*N [L] 1198* At this point, L is resized as M by chunk. 1199* 1200 DO 30 I = 1, N, CHUNK 1201 BLK = MIN( N - I + 1, CHUNK ) 1202 CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M, 1203 $ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL ) 1204 CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL, 1205 $ A( 1, I ), LDA ) 1206 30 CONTINUE 1207* 1208 ELSE IF( WNTQS ) THEN 1209* 1210* Path 3t (N >> M, JOBZ='S') 1211* M right singular vectors to be computed in VT and 1212* M left singular vectors to be computed in U 1213* 1214 IL = 1 1215* 1216* WORK(IL) is M by M 1217* 1218 LDWRKL = M 1219 ITAU = IL + LDWRKL*M 1220 NWORK = ITAU + M 1221* 1222* Compute A=L*Q 1223* Workspace: need M*M [L] + M [tau] + M [work] 1224* Workspace: prefer M*M [L] + M [tau] + M*NB [work] 1225* 1226 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1227 $ LWORK - NWORK + 1, IERR ) 1228* 1229* Copy L to WORK(IL), zeroing out above it 1230* 1231 CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 1232 CALL DLASET( 'U', M - 1, M - 1, ZERO, ZERO, 1233 $ WORK( IL + LDWRKL ), LDWRKL ) 1234* 1235* Generate Q in A 1236* Workspace: need M*M [L] + M [tau] + M [work] 1237* Workspace: prefer M*M [L] + M [tau] + M*NB [work] 1238* 1239 CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ), 1240 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1241 IE = ITAU 1242 ITAUQ = IE + M 1243 ITAUP = ITAUQ + M 1244 NWORK = ITAUP + M 1245* 1246* Bidiagonalize L in WORK(IU). 1247* Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work] 1248* Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + 2*M*NB [work] 1249* 1250 CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ), 1251 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 1252 $ LWORK - NWORK + 1, IERR ) 1253* 1254* Perform bidiagonal SVD, computing left singular vectors 1255* of bidiagonal matrix in U and computing right singular 1256* vectors of bidiagonal matrix in VT 1257* Workspace: need M*M [L] + 3*M [e, tauq, taup] + BDSPAC 1258* 1259 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT, 1260 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 1261 $ INFO ) 1262* 1263* Overwrite U by left singular vectors of L and VT 1264* by right singular vectors of L 1265* Workspace: need M*M [L] + 3*M [e, tauq, taup] + M [work] 1266* Workspace: prefer M*M [L] + 3*M [e, tauq, taup] + M*NB [work] 1267* 1268 CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 1269 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1270 $ LWORK - NWORK + 1, IERR ) 1271 CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL, 1272 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1273 $ LWORK - NWORK + 1, IERR ) 1274* 1275* Multiply right singular vectors of L in WORK(IL) by 1276* Q in A, storing result in VT 1277* Workspace: need M*M [L] 1278* 1279 CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) 1280 CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL, 1281 $ A, LDA, ZERO, VT, LDVT ) 1282* 1283 ELSE IF( WNTQA ) THEN 1284* 1285* Path 4t (N >> M, JOBZ='A') 1286* N right singular vectors to be computed in VT and 1287* M left singular vectors to be computed in U 1288* 1289 IVT = 1 1290* 1291* WORK(IVT) is M by M 1292* 1293 LDWKVT = M 1294 ITAU = IVT + LDWKVT*M 1295 NWORK = ITAU + M 1296* 1297* Compute A=L*Q, copying result to VT 1298* Workspace: need M*M [VT] + M [tau] + M [work] 1299* Workspace: prefer M*M [VT] + M [tau] + M*NB [work] 1300* 1301 CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 1302 $ LWORK - NWORK + 1, IERR ) 1303 CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT ) 1304* 1305* Generate Q in VT 1306* Workspace: need M*M [VT] + M [tau] + N [work] 1307* Workspace: prefer M*M [VT] + M [tau] + N*NB [work] 1308* 1309 CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 1310 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1311* 1312* Produce L in A, zeroing out other entries 1313* 1314 CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA ) 1315 IE = ITAU 1316 ITAUQ = IE + M 1317 ITAUP = ITAUQ + M 1318 NWORK = ITAUP + M 1319* 1320* Bidiagonalize L in A 1321* Workspace: need M*M [VT] + 3*M [e, tauq, taup] + M [work] 1322* Workspace: prefer M*M [VT] + 3*M [e, tauq, taup] + 2*M*NB [work] 1323* 1324 CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 1325 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1326 $ IERR ) 1327* 1328* Perform bidiagonal SVD, computing left singular vectors 1329* of bidiagonal matrix in U and computing right singular 1330* vectors of bidiagonal matrix in WORK(IVT) 1331* Workspace: need M*M [VT] + 3*M [e, tauq, taup] + BDSPAC 1332* 1333 CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, 1334 $ WORK( IVT ), LDWKVT, DUM, IDUM, 1335 $ WORK( NWORK ), IWORK, INFO ) 1336* 1337* Overwrite U by left singular vectors of L and WORK(IVT) 1338* by right singular vectors of L 1339* Workspace: need M*M [VT] + 3*M [e, tauq, taup]+ M [work] 1340* Workspace: prefer M*M [VT] + 3*M [e, tauq, taup]+ M*NB [work] 1341* 1342 CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA, 1343 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1344 $ LWORK - NWORK + 1, IERR ) 1345 CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA, 1346 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 1347 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1348* 1349* Multiply right singular vectors of L in WORK(IVT) by 1350* Q in VT, storing result in A 1351* Workspace: need M*M [VT] 1352* 1353 CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT, 1354 $ VT, LDVT, ZERO, A, LDA ) 1355* 1356* Copy right singular vectors of A from A to VT 1357* 1358 CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT ) 1359* 1360 END IF 1361* 1362 ELSE 1363* 1364* N .LT. MNTHR 1365* 1366* Path 5t (N > M, but not much larger) 1367* Reduce to bidiagonal form without LQ decomposition 1368* 1369 IE = 1 1370 ITAUQ = IE + M 1371 ITAUP = ITAUQ + M 1372 NWORK = ITAUP + M 1373* 1374* Bidiagonalize A 1375* Workspace: need 3*M [e, tauq, taup] + N [work] 1376* Workspace: prefer 3*M [e, tauq, taup] + (M+N)*NB [work] 1377* 1378 CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ), 1379 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 1380 $ IERR ) 1381 IF( WNTQN ) THEN 1382* 1383* Path 5tn (N > M, JOBZ='N') 1384* Perform bidiagonal SVD, only computing singular values 1385* Workspace: need 3*M [e, tauq, taup] + BDSPAC 1386* 1387 CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1, 1388 $ DUM, IDUM, WORK( NWORK ), IWORK, INFO ) 1389 ELSE IF( WNTQO ) THEN 1390* Path 5to (N > M, JOBZ='O') 1391 LDWKVT = M 1392 IVT = NWORK 1393 IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN 1394* 1395* WORK( IVT ) is M by N 1396* 1397 CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ), 1398 $ LDWKVT ) 1399 NWORK = IVT + LDWKVT*N 1400* IL is unused; silence compile warnings 1401 IL = -1 1402 ELSE 1403* 1404* WORK( IVT ) is M by M 1405* 1406 NWORK = IVT + LDWKVT*M 1407 IL = NWORK 1408* 1409* WORK(IL) is M by CHUNK 1410* 1411 CHUNK = ( LWORK - M*M - 3*M ) / M 1412 END IF 1413* 1414* Perform bidiagonal SVD, computing left singular vectors 1415* of bidiagonal matrix in U and computing right singular 1416* vectors of bidiagonal matrix in WORK(IVT) 1417* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + BDSPAC 1418* 1419 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, 1420 $ WORK( IVT ), LDWKVT, DUM, IDUM, 1421 $ WORK( NWORK ), IWORK, INFO ) 1422* 1423* Overwrite U by left singular vectors of A 1424* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work] 1425* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work] 1426* 1427 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1428 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1429 $ LWORK - NWORK + 1, IERR ) 1430* 1431 IF( LWORK .GE. M*N + 3*M + BDSPAC ) THEN 1432* 1433* Path 5to-fast 1434* Overwrite WORK(IVT) by left singular vectors of A 1435* Workspace: need 3*M [e, tauq, taup] + M*N [VT] + M [work] 1436* Workspace: prefer 3*M [e, tauq, taup] + M*N [VT] + M*NB [work] 1437* 1438 CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA, 1439 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 1440 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1441* 1442* Copy right singular vectors of A from WORK(IVT) to A 1443* 1444 CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) 1445 ELSE 1446* 1447* Path 5to-slow 1448* Generate P**T in A 1449* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M [work] 1450* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*NB [work] 1451* 1452 CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 1453 $ WORK( NWORK ), LWORK - NWORK + 1, IERR ) 1454* 1455* Multiply Q in A by right singular vectors of 1456* bidiagonal matrix in WORK(IVT), storing result in 1457* WORK(IL) and copying to A 1458* Workspace: need 3*M [e, tauq, taup] + M*M [VT] + M*NB [L] 1459* Workspace: prefer 3*M [e, tauq, taup] + M*M [VT] + M*N [L] 1460* 1461 DO 40 I = 1, N, CHUNK 1462 BLK = MIN( N - I + 1, CHUNK ) 1463 CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), 1464 $ LDWKVT, A( 1, I ), LDA, ZERO, 1465 $ WORK( IL ), M ) 1466 CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ), 1467 $ LDA ) 1468 40 CONTINUE 1469 END IF 1470 ELSE IF( WNTQS ) THEN 1471* 1472* Path 5ts (N > M, JOBZ='S') 1473* Perform bidiagonal SVD, computing left singular vectors 1474* of bidiagonal matrix in U and computing right singular 1475* vectors of bidiagonal matrix in VT 1476* Workspace: need 3*M [e, tauq, taup] + BDSPAC 1477* 1478 CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT ) 1479 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT, 1480 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 1481 $ INFO ) 1482* 1483* Overwrite U by left singular vectors of A and VT 1484* by right singular vectors of A 1485* Workspace: need 3*M [e, tauq, taup] + M [work] 1486* Workspace: prefer 3*M [e, tauq, taup] + M*NB [work] 1487* 1488 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1489 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1490 $ LWORK - NWORK + 1, IERR ) 1491 CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA, 1492 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1493 $ LWORK - NWORK + 1, IERR ) 1494 ELSE IF( WNTQA ) THEN 1495* 1496* Path 5ta (N > M, JOBZ='A') 1497* Perform bidiagonal SVD, computing left singular vectors 1498* of bidiagonal matrix in U and computing right singular 1499* vectors of bidiagonal matrix in VT 1500* Workspace: need 3*M [e, tauq, taup] + BDSPAC 1501* 1502 CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT ) 1503 CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT, 1504 $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK, 1505 $ INFO ) 1506* 1507* Set the right corner of VT to identity matrix 1508* 1509 IF( N.GT.M ) THEN 1510 CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT(M+1,M+1), 1511 $ LDVT ) 1512 END IF 1513* 1514* Overwrite U by left singular vectors of A and VT 1515* by right singular vectors of A 1516* Workspace: need 3*M [e, tauq, taup] + N [work] 1517* Workspace: prefer 3*M [e, tauq, taup] + N*NB [work] 1518* 1519 CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 1520 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 1521 $ LWORK - NWORK + 1, IERR ) 1522 CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA, 1523 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 1524 $ LWORK - NWORK + 1, IERR ) 1525 END IF 1526* 1527 END IF 1528* 1529 END IF 1530* 1531* Undo scaling if necessary 1532* 1533 IF( ISCL.EQ.1 ) THEN 1534 IF( ANRM.GT.BIGNUM ) 1535 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 1536 $ IERR ) 1537 IF( ANRM.LT.SMLNUM ) 1538 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 1539 $ IERR ) 1540 END IF 1541* 1542* Return optimal workspace in WORK(1) 1543* 1544 WORK( 1 ) = MAXWRK 1545* 1546 RETURN 1547* 1548* End of DGESDD 1549* 1550 END 1551