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