1*> \brief <b> SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SGESVDQ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgesvdq.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgesvdq.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgesvdq.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, 22* S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, 23* WORK, LWORK, RWORK, LRWORK, INFO ) 24* 25* .. Scalar Arguments .. 26* IMPLICIT NONE 27* CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV 28* INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK, 29* INFO 30* .. 31* .. Array Arguments .. 32* REAL A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * ) 33* REAL S( * ), RWORK( * ) 34* INTEGER IWORK( * ) 35* .. 36* 37* 38*> \par Purpose: 39* ============= 40*> 41*> \verbatim 42*> 43*> SGESVDQ computes the singular value decomposition (SVD) of a real 44*> M-by-N matrix A, where M >= N. The SVD of A is written as 45*> [++] [xx] [x0] [xx] 46*> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx] 47*> [++] [xx] 48*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal 49*> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements 50*> of SIGMA are the singular values of A. The columns of U and V are the 51*> left and the right singular vectors of A, respectively. 52*> \endverbatim 53* 54* Arguments: 55* ========== 56* 57*> \param[in] JOBA 58*> \verbatim 59*> JOBA is CHARACTER*1 60*> Specifies the level of accuracy in the computed SVD 61*> = 'A' The requested accuracy corresponds to having the backward 62*> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F, 63*> where EPS = SLAMCH('Epsilon'). This authorises CGESVDQ to 64*> truncate the computed triangular factor in a rank revealing 65*> QR factorization whenever the truncated part is below the 66*> threshold of the order of EPS * ||A||_F. This is aggressive 67*> truncation level. 68*> = 'M' Similarly as with 'A', but the truncation is more gentle: it 69*> is allowed only when there is a drop on the diagonal of the 70*> triangular factor in the QR factorization. This is medium 71*> truncation level. 72*> = 'H' High accuracy requested. No numerical rank determination based 73*> on the rank revealing QR factorization is attempted. 74*> = 'E' Same as 'H', and in addition the condition number of column 75*> scaled A is estimated and returned in RWORK(1). 76*> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1) 77*> \endverbatim 78*> 79*> \param[in] JOBP 80*> \verbatim 81*> JOBP is CHARACTER*1 82*> = 'P' The rows of A are ordered in decreasing order with respect to 83*> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost 84*> of extra data movement. Recommended for numerical robustness. 85*> = 'N' No row pivoting. 86*> \endverbatim 87*> 88*> \param[in] JOBR 89*> \verbatim 90*> JOBR is CHARACTER*1 91*> = 'T' After the initial pivoted QR factorization, SGESVD is applied to 92*> the transposed R**T of the computed triangular factor R. This involves 93*> some extra data movement (matrix transpositions). Useful for 94*> experiments, research and development. 95*> = 'N' The triangular factor R is given as input to SGESVD. This may be 96*> preferred as it involves less data movement. 97*> \endverbatim 98*> 99*> \param[in] JOBU 100*> \verbatim 101*> JOBU is CHARACTER*1 102*> = 'A' All M left singular vectors are computed and returned in the 103*> matrix U. See the description of U. 104*> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned 105*> in the matrix U. See the description of U. 106*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular 107*> vectors are computed and returned in the matrix U. 108*> = 'F' The N left singular vectors are returned in factored form as the 109*> product of the Q factor from the initial QR factorization and the 110*> N left singular vectors of (R**T , 0)**T. If row pivoting is used, 111*> then the necessary information on the row pivoting is stored in 112*> IWORK(N+1:N+M-1). 113*> = 'N' The left singular vectors are not computed. 114*> \endverbatim 115*> 116*> \param[in] JOBV 117*> \verbatim 118*> JOBV is CHARACTER*1 119*> = 'A', 'V' All N right singular vectors are computed and returned in 120*> the matrix V. 121*> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular 122*> vectors are computed and returned in the matrix V. This option is 123*> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal. 124*> = 'N' The right singular vectors are not computed. 125*> \endverbatim 126*> 127*> \param[in] M 128*> \verbatim 129*> M is INTEGER 130*> The number of rows of the input matrix A. M >= 0. 131*> \endverbatim 132*> 133*> \param[in] N 134*> \verbatim 135*> N is INTEGER 136*> The number of columns of the input matrix A. M >= N >= 0. 137*> \endverbatim 138*> 139*> \param[in,out] A 140*> \verbatim 141*> A is REAL array of dimensions LDA x N 142*> On entry, the input matrix A. 143*> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains 144*> the Householder vectors as stored by SGEQP3. If JOBU = 'F', these Householder 145*> vectors together with WORK(1:N) can be used to restore the Q factors from 146*> the initial pivoted QR factorization of A. See the description of U. 147*> \endverbatim 148*> 149*> \param[in] LDA 150*> \verbatim 151*> LDA is INTEGER. 152*> The leading dimension of the array A. LDA >= max(1,M). 153*> \endverbatim 154*> 155*> \param[out] S 156*> \verbatim 157*> S is REAL array of dimension N. 158*> The singular values of A, ordered so that S(i) >= S(i+1). 159*> \endverbatim 160*> 161*> \param[out] U 162*> \verbatim 163*> U is REAL array, dimension 164*> LDU x M if JOBU = 'A'; see the description of LDU. In this case, 165*> on exit, U contains the M left singular vectors. 166*> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this 167*> case, U contains the leading N or the leading NUMRANK left singular vectors. 168*> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U 169*> contains N x N orthogonal matrix that can be used to form the left 170*> singular vectors. 171*> If JOBU = 'N', U is not referenced. 172*> \endverbatim 173*> 174*> \param[in] LDU 175*> \verbatim 176*> LDU is INTEGER. 177*> The leading dimension of the array U. 178*> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M). 179*> If JOBU = 'F', LDU >= max(1,N). 180*> Otherwise, LDU >= 1. 181*> \endverbatim 182*> 183*> \param[out] V 184*> \verbatim 185*> V is REAL array, dimension 186*> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' . 187*> If JOBV = 'A', or 'V', V contains the N-by-N orthogonal matrix V**T; 188*> If JOBV = 'R', V contains the first NUMRANK rows of V**T (the right 189*> singular vectors, stored rowwise, of the NUMRANK largest singular values). 190*> If JOBV = 'N' and JOBA = 'E', V is used as a workspace. 191*> If JOBV = 'N', and JOBA.NE.'E', V is not referenced. 192*> \endverbatim 193*> 194*> \param[in] LDV 195*> \verbatim 196*> LDV is INTEGER 197*> The leading dimension of the array V. 198*> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N). 199*> Otherwise, LDV >= 1. 200*> \endverbatim 201*> 202*> \param[out] NUMRANK 203*> \verbatim 204*> NUMRANK is INTEGER 205*> NUMRANK is the numerical rank first determined after the rank 206*> revealing QR factorization, following the strategy specified by the 207*> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK 208*> leading singular values and vectors are then requested in the call 209*> of SGESVD. The final value of NUMRANK might be further reduced if 210*> some singular values are computed as zeros. 211*> \endverbatim 212*> 213*> \param[out] IWORK 214*> \verbatim 215*> IWORK is INTEGER array, dimension (max(1, LIWORK)). 216*> On exit, IWORK(1:N) contains column pivoting permutation of the 217*> rank revealing QR factorization. 218*> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence 219*> of row swaps used in row pivoting. These can be used to restore the 220*> left singular vectors in the case JOBU = 'F'. 221*> 222*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, 223*> IWORK(1) returns the minimal LIWORK. 224*> \endverbatim 225*> 226*> \param[in] LIWORK 227*> \verbatim 228*> LIWORK is INTEGER 229*> The dimension of the array IWORK. 230*> LIWORK >= N + M - 1, if JOBP = 'P' and JOBA .NE. 'E'; 231*> LIWORK >= N if JOBP = 'N' and JOBA .NE. 'E'; 232*> LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E'; 233*> LIWORK >= N + N if JOBP = 'N' and JOBA = 'E'. 234*> 235*> If LIWORK = -1, then a workspace query is assumed; the routine 236*> only calculates and returns the optimal and minimal sizes 237*> for the WORK, IWORK, and RWORK arrays, and no error 238*> message related to LWORK is issued by XERBLA. 239*> \endverbatim 240*> 241*> \param[out] WORK 242*> \verbatim 243*> WORK is REAL array, dimension (max(2, LWORK)), used as a workspace. 244*> On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters 245*> needed to recover the Q factor from the QR factorization computed by 246*> SGEQP3. 247*> 248*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, 249*> WORK(1) returns the optimal LWORK, and 250*> WORK(2) returns the minimal LWORK. 251*> \endverbatim 252*> 253*> \param[in,out] LWORK 254*> \verbatim 255*> LWORK is INTEGER 256*> The dimension of the array WORK. It is determined as follows: 257*> Let LWQP3 = 3*N+1, LWCON = 3*N, and let 258*> LWORQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U' 259*> { MAX( M, 1 ), if JOBU = 'A' 260*> LWSVD = MAX( 5*N, 1 ) 261*> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ), 262*> LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 ) 263*> Then the minimal value of LWORK is: 264*> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed; 265*> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed, 266*> and a scaled condition estimate requested; 267*> 268*> = N + MAX( LWQP3, LWSVD, LWORQ ) if the singular values and the left 269*> singular vectors are requested; 270*> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular values and the left 271*> singular vectors are requested, and also 272*> a scaled condition estimate requested; 273*> 274*> = N + MAX( LWQP3, LWSVD ) if the singular values and the right 275*> singular vectors are requested; 276*> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right 277*> singular vectors are requested, and also 278*> a scaled condition etimate requested; 279*> 280*> = N + MAX( LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV = 'R'; 281*> independent of JOBR; 282*> = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested, 283*> JOBV = 'R' and, also a scaled condition 284*> estimate requested; independent of JOBR; 285*> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ), 286*> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the 287*> full SVD is requested with JOBV = 'A' or 'V', and 288*> JOBR ='N' 289*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ), 290*> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) ) 291*> if the full SVD is requested with JOBV = 'A' or 'V', and 292*> JOBR ='N', and also a scaled condition number estimate 293*> requested. 294*> = MAX( N + MAX( LWQP3, LWSVD, LWORQ ), 295*> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the 296*> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T' 297*> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ), 298*> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) 299*> if the full SVD is requested with JOBV = 'A' or 'V', and 300*> JOBR ='T', and also a scaled condition number estimate 301*> requested. 302*> Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ). 303*> 304*> If LWORK = -1, then a workspace query is assumed; the routine 305*> only calculates and returns the optimal and minimal sizes 306*> for the WORK, IWORK, and RWORK arrays, and no error 307*> message related to LWORK is issued by XERBLA. 308*> \endverbatim 309*> 310*> \param[out] RWORK 311*> \verbatim 312*> RWORK is REAL array, dimension (max(1, LRWORK)). 313*> On exit, 314*> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition 315*> number of column scaled A. If A = C * D where D is diagonal and C 316*> has unit columns in the Euclidean norm, then, assuming full column rank, 317*> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1). 318*> Otherwise, RWORK(1) = -1. 319*> 2. RWORK(2) contains the number of singular values computed as 320*> exact zeros in SGESVD applied to the upper triangular or trapezoidal 321*> R (from the initial QR factorization). In case of early exit (no call to 322*> SGESVD, such as in the case of zero matrix) RWORK(2) = -1. 323*> 324*> If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0, 325*> RWORK(1) returns the minimal LRWORK. 326*> \endverbatim 327*> 328*> \param[in] LRWORK 329*> \verbatim 330*> LRWORK is INTEGER. 331*> The dimension of the array RWORK. 332*> If JOBP ='P', then LRWORK >= MAX(2, M). 333*> Otherwise, LRWORK >= 2 334*> 335*> If LRWORK = -1, then a workspace query is assumed; the routine 336*> only calculates and returns the optimal and minimal sizes 337*> for the WORK, IWORK, and RWORK arrays, and no error 338*> message related to LWORK is issued by XERBLA. 339*> \endverbatim 340*> 341*> \param[out] INFO 342*> \verbatim 343*> INFO is INTEGER 344*> = 0: successful exit. 345*> < 0: if INFO = -i, the i-th argument had an illegal value. 346*> > 0: if SBDSQR did not converge, INFO specifies how many superdiagonals 347*> of an intermediate bidiagonal form B (computed in SGESVD) did not 348*> converge to zero. 349*> \endverbatim 350* 351*> \par Further Details: 352* ======================== 353*> 354*> \verbatim 355*> 356*> 1. The data movement (matrix transpose) is coded using simple nested 357*> DO-loops because BLAS and LAPACK do not provide corresponding subroutines. 358*> Those DO-loops are easily identified in this source code - by the CONTINUE 359*> statements labeled with 11**. In an optimized version of this code, the 360*> nested DO loops should be replaced with calls to an optimized subroutine. 361*> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause 362*> column norm overflow. This is the minial precaution and it is left to the 363*> SVD routine (CGESVD) to do its own preemptive scaling if potential over- 364*> or underflows are detected. To avoid repeated scanning of the array A, 365*> an optimal implementation would do all necessary scaling before calling 366*> CGESVD and the scaling in CGESVD can be switched off. 367*> 3. Other comments related to code optimization are given in comments in the 368*> code, enlosed in [[double brackets]]. 369*> \endverbatim 370* 371*> \par Bugs, examples and comments 372* =========================== 373* 374*> \verbatim 375*> Please report all bugs and send interesting examples and/or comments to 376*> drmac@math.hr. Thank you. 377*> \endverbatim 378* 379*> \par References 380* =============== 381* 382*> \verbatim 383*> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for 384*> Computing the SVD with High Accuracy. ACM Trans. Math. Softw. 385*> 44(1): 11:1-11:30 (2017) 386*> 387*> SIGMA library, xGESVDQ section updated February 2016. 388*> Developed and coded by Zlatko Drmac, Department of Mathematics 389*> University of Zagreb, Croatia, drmac@math.hr 390*> \endverbatim 391* 392* 393*> \par Contributors: 394* ================== 395*> 396*> \verbatim 397*> Developed and coded by Zlatko Drmac, Department of Mathematics 398*> University of Zagreb, Croatia, drmac@math.hr 399*> \endverbatim 400* 401* Authors: 402* ======== 403* 404*> \author Univ. of Tennessee 405*> \author Univ. of California Berkeley 406*> \author Univ. of Colorado Denver 407*> \author NAG Ltd. 408* 409*> \ingroup realGEsing 410* 411* ===================================================================== 412 SUBROUTINE SGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, 413 $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, 414 $ WORK, LWORK, RWORK, LRWORK, INFO ) 415* .. Scalar Arguments .. 416 IMPLICIT NONE 417 CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV 418 INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK, 419 $ INFO 420* .. 421* .. Array Arguments .. 422 REAL A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * ) 423 REAL S( * ), RWORK( * ) 424 INTEGER IWORK( * ) 425* 426* ===================================================================== 427* 428* .. Parameters .. 429 REAL ZERO, ONE 430 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 431* .. 432* .. Local Scalars .. 433 INTEGER IERR, IWOFF, NR, N1, OPTRATIO, p, q 434 INTEGER LWCON, LWQP3, LWRK_SGELQF, LWRK_SGESVD, LWRK_SGESVD2, 435 $ LWRK_SGEQP3, LWRK_SGEQRF, LWRK_SORMLQ, LWRK_SORMQR, 436 $ LWRK_SORMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWORQ, 437 $ LWORQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2, 438 $ IMINWRK, RMINWRK 439 LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV, 440 $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA, 441 $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR 442 REAL BIG, EPSLN, RTMP, SCONDA, SFMIN 443* .. 444* .. Local Arrays 445 REAL RDUMMY(1) 446* .. 447* .. External Subroutines (BLAS, LAPACK) 448 EXTERNAL SGELQF, SGEQP3, SGEQRF, SGESVD, SLACPY, SLAPMT, 449 $ SLASCL, SLASET, SLASWP, SSCAL, SPOCON, SORMLQ, 450 $ SORMQR, XERBLA 451* .. 452* .. External Functions (BLAS, LAPACK) 453 LOGICAL LSAME 454 INTEGER ISAMAX 455 REAL SLANGE, SNRM2, SLAMCH 456 EXTERNAL SLANGE, LSAME, ISAMAX, SNRM2, SLAMCH 457* .. 458* .. Intrinsic Functions .. 459 INTRINSIC ABS, MAX, MIN, REAL, SQRT 460* .. 461* .. Executable Statements .. 462* 463* Test the input arguments 464* 465 WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' ) 466 WNTUR = LSAME( JOBU, 'R' ) 467 WNTUA = LSAME( JOBU, 'A' ) 468 WNTUF = LSAME( JOBU, 'F' ) 469 LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA 470 LSVEC = LSVC0 .OR. WNTUF 471 DNTWU = LSAME( JOBU, 'N' ) 472* 473 WNTVR = LSAME( JOBV, 'R' ) 474 WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' ) 475 RSVEC = WNTVR .OR. WNTVA 476 DNTWV = LSAME( JOBV, 'N' ) 477* 478 ACCLA = LSAME( JOBA, 'A' ) 479 ACCLM = LSAME( JOBA, 'M' ) 480 CONDA = LSAME( JOBA, 'E' ) 481 ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA 482* 483 ROWPRM = LSAME( JOBP, 'P' ) 484 RTRANS = LSAME( JOBR, 'T' ) 485* 486 IF ( ROWPRM ) THEN 487 IF ( CONDA ) THEN 488 IMINWRK = MAX( 1, N + M - 1 + N ) 489 ELSE 490 IMINWRK = MAX( 1, N + M - 1 ) 491 END IF 492 RMINWRK = MAX( 2, M ) 493 ELSE 494 IF ( CONDA ) THEN 495 IMINWRK = MAX( 1, N + N ) 496 ELSE 497 IMINWRK = MAX( 1, N ) 498 END IF 499 RMINWRK = 2 500 END IF 501 LQUERY = (LIWORK .EQ. -1 .OR. LWORK .EQ. -1 .OR. LRWORK .EQ. -1) 502 INFO = 0 503 IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN 504 INFO = -1 505 ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN 506 INFO = -2 507 ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN 508 INFO = -3 509 ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN 510 INFO = -4 511 ELSE IF ( WNTUR .AND. WNTVA ) THEN 512 INFO = -5 513 ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN 514 INFO = -5 515 ELSE IF ( M.LT.0 ) THEN 516 INFO = -6 517 ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN 518 INFO = -7 519 ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN 520 INFO = -9 521 ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR. 522 $ ( WNTUF .AND. LDU.LT.N ) ) THEN 523 INFO = -12 524 ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR. 525 $ ( CONDA .AND. LDV.LT.N ) ) THEN 526 INFO = -14 527 ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN 528 INFO = -17 529 END IF 530* 531* 532 IF ( INFO .EQ. 0 ) THEN 533* .. compute the minimal and the optimal workspace lengths 534* [[The expressions for computing the minimal and the optimal 535* values of LWORK are written with a lot of redundancy and 536* can be simplified. However, this detailed form is easier for 537* maintenance and modifications of the code.]] 538* 539* .. minimal workspace length for SGEQP3 of an M x N matrix 540 LWQP3 = 3 * N + 1 541* .. minimal workspace length for SORMQR to build left singular vectors 542 IF ( WNTUS .OR. WNTUR ) THEN 543 LWORQ = MAX( N , 1 ) 544 ELSE IF ( WNTUA ) THEN 545 LWORQ = MAX( M , 1 ) 546 END IF 547* .. minimal workspace length for SPOCON of an N x N matrix 548 LWCON = 3 * N 549* .. SGESVD of an N x N matrix 550 LWSVD = MAX( 5 * N, 1 ) 551 IF ( LQUERY ) THEN 552 CALL SGEQP3( M, N, A, LDA, IWORK, RDUMMY, RDUMMY, -1, 553 $ IERR ) 554 LWRK_SGEQP3 = INT( RDUMMY(1) ) 555 IF ( WNTUS .OR. WNTUR ) THEN 556 CALL SORMQR( 'L', 'N', M, N, N, A, LDA, RDUMMY, U, 557 $ LDU, RDUMMY, -1, IERR ) 558 LWRK_SORMQR = INT( RDUMMY(1) ) 559 ELSE IF ( WNTUA ) THEN 560 CALL SORMQR( 'L', 'N', M, M, N, A, LDA, RDUMMY, U, 561 $ LDU, RDUMMY, -1, IERR ) 562 LWRK_SORMQR = INT( RDUMMY(1) ) 563 ELSE 564 LWRK_SORMQR = 0 565 END IF 566 END IF 567 MINWRK = 2 568 OPTWRK = 2 569 IF ( .NOT. (LSVEC .OR. RSVEC )) THEN 570* .. minimal and optimal sizes of the workspace if 571* only the singular values are requested 572 IF ( CONDA ) THEN 573 MINWRK = MAX( N+LWQP3, LWCON, LWSVD ) 574 ELSE 575 MINWRK = MAX( N+LWQP3, LWSVD ) 576 END IF 577 IF ( LQUERY ) THEN 578 CALL SGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU, 579 $ V, LDV, RDUMMY, -1, IERR ) 580 LWRK_SGESVD = INT( RDUMMY(1) ) 581 IF ( CONDA ) THEN 582 OPTWRK = MAX( N+LWRK_SGEQP3, N+LWCON, LWRK_SGESVD ) 583 ELSE 584 OPTWRK = MAX( N+LWRK_SGEQP3, LWRK_SGESVD ) 585 END IF 586 END IF 587 ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN 588* .. minimal and optimal sizes of the workspace if the 589* singular values and the left singular vectors are requested 590 IF ( CONDA ) THEN 591 MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) 592 ELSE 593 MINWRK = N + MAX( LWQP3, LWSVD, LWORQ ) 594 END IF 595 IF ( LQUERY ) THEN 596 IF ( RTRANS ) THEN 597 CALL SGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU, 598 $ V, LDV, RDUMMY, -1, IERR ) 599 ELSE 600 CALL SGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU, 601 $ V, LDV, RDUMMY, -1, IERR ) 602 END IF 603 LWRK_SGESVD = INT( RDUMMY(1) ) 604 IF ( CONDA ) THEN 605 OPTWRK = N + MAX( LWRK_SGEQP3, LWCON, LWRK_SGESVD, 606 $ LWRK_SORMQR ) 607 ELSE 608 OPTWRK = N + MAX( LWRK_SGEQP3, LWRK_SGESVD, 609 $ LWRK_SORMQR ) 610 END IF 611 END IF 612 ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN 613* .. minimal and optimal sizes of the workspace if the 614* singular values and the right singular vectors are requested 615 IF ( CONDA ) THEN 616 MINWRK = N + MAX( LWQP3, LWCON, LWSVD ) 617 ELSE 618 MINWRK = N + MAX( LWQP3, LWSVD ) 619 END IF 620 IF ( LQUERY ) THEN 621 IF ( RTRANS ) THEN 622 CALL SGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU, 623 $ V, LDV, RDUMMY, -1, IERR ) 624 ELSE 625 CALL SGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU, 626 $ V, LDV, RDUMMY, -1, IERR ) 627 END IF 628 LWRK_SGESVD = INT( RDUMMY(1) ) 629 IF ( CONDA ) THEN 630 OPTWRK = N + MAX( LWRK_SGEQP3, LWCON, LWRK_SGESVD ) 631 ELSE 632 OPTWRK = N + MAX( LWRK_SGEQP3, LWRK_SGESVD ) 633 END IF 634 END IF 635 ELSE 636* .. minimal and optimal sizes of the workspace if the 637* full SVD is requested 638 IF ( RTRANS ) THEN 639 MINWRK = MAX( LWQP3, LWSVD, LWORQ ) 640 IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON ) 641 MINWRK = MINWRK + N 642 IF ( WNTVA ) THEN 643* .. minimal workspace length for N x N/2 SGEQRF 644 LWQRF = MAX( N/2, 1 ) 645* .. minimal workspace length for N/2 x N/2 SGESVD 646 LWSVD2 = MAX( 5 * (N/2), 1 ) 647 LWORQ2 = MAX( N, 1 ) 648 MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, 649 $ N/2+LWORQ2, LWORQ ) 650 IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON ) 651 MINWRK2 = N + MINWRK2 652 MINWRK = MAX( MINWRK, MINWRK2 ) 653 END IF 654 ELSE 655 MINWRK = MAX( LWQP3, LWSVD, LWORQ ) 656 IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON ) 657 MINWRK = MINWRK + N 658 IF ( WNTVA ) THEN 659* .. minimal workspace length for N/2 x N SGELQF 660 LWLQF = MAX( N/2, 1 ) 661 LWSVD2 = MAX( 5 * (N/2), 1 ) 662 LWUNLQ = MAX( N , 1 ) 663 MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, 664 $ N/2+LWUNLQ, LWORQ ) 665 IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON ) 666 MINWRK2 = N + MINWRK2 667 MINWRK = MAX( MINWRK, MINWRK2 ) 668 END IF 669 END IF 670 IF ( LQUERY ) THEN 671 IF ( RTRANS ) THEN 672 CALL SGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU, 673 $ V, LDV, RDUMMY, -1, IERR ) 674 LWRK_SGESVD = INT( RDUMMY(1) ) 675 OPTWRK = MAX(LWRK_SGEQP3,LWRK_SGESVD,LWRK_SORMQR) 676 IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON ) 677 OPTWRK = N + OPTWRK 678 IF ( WNTVA ) THEN 679 CALL SGEQRF(N,N/2,U,LDU,RDUMMY,RDUMMY,-1,IERR) 680 LWRK_SGEQRF = INT( RDUMMY(1) ) 681 CALL SGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU, 682 $ V, LDV, RDUMMY, -1, IERR ) 683 LWRK_SGESVD2 = INT( RDUMMY(1) ) 684 CALL SORMQR( 'R', 'C', N, N, N/2, U, LDU, RDUMMY, 685 $ V, LDV, RDUMMY, -1, IERR ) 686 LWRK_SORMQR2 = INT( RDUMMY(1) ) 687 OPTWRK2 = MAX( LWRK_SGEQP3, N/2+LWRK_SGEQRF, 688 $ N/2+LWRK_SGESVD2, N/2+LWRK_SORMQR2 ) 689 IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON ) 690 OPTWRK2 = N + OPTWRK2 691 OPTWRK = MAX( OPTWRK, OPTWRK2 ) 692 END IF 693 ELSE 694 CALL SGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU, 695 $ V, LDV, RDUMMY, -1, IERR ) 696 LWRK_SGESVD = INT( RDUMMY(1) ) 697 OPTWRK = MAX(LWRK_SGEQP3,LWRK_SGESVD,LWRK_SORMQR) 698 IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON ) 699 OPTWRK = N + OPTWRK 700 IF ( WNTVA ) THEN 701 CALL SGELQF(N/2,N,U,LDU,RDUMMY,RDUMMY,-1,IERR) 702 LWRK_SGELQF = INT( RDUMMY(1) ) 703 CALL SGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU, 704 $ V, LDV, RDUMMY, -1, IERR ) 705 LWRK_SGESVD2 = INT( RDUMMY(1) ) 706 CALL SORMLQ( 'R', 'N', N, N, N/2, U, LDU, RDUMMY, 707 $ V, LDV, RDUMMY,-1,IERR ) 708 LWRK_SORMLQ = INT( RDUMMY(1) ) 709 OPTWRK2 = MAX( LWRK_SGEQP3, N/2+LWRK_SGELQF, 710 $ N/2+LWRK_SGESVD2, N/2+LWRK_SORMLQ ) 711 IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON ) 712 OPTWRK2 = N + OPTWRK2 713 OPTWRK = MAX( OPTWRK, OPTWRK2 ) 714 END IF 715 END IF 716 END IF 717 END IF 718* 719 MINWRK = MAX( 2, MINWRK ) 720 OPTWRK = MAX( 2, OPTWRK ) 721 IF ( LWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19 722* 723 END IF 724* 725 IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN 726 INFO = -21 727 END IF 728 IF( INFO.NE.0 ) THEN 729 CALL XERBLA( 'SGESVDQ', -INFO ) 730 RETURN 731 ELSE IF ( LQUERY ) THEN 732* 733* Return optimal workspace 734* 735 IWORK(1) = IMINWRK 736 WORK(1) = OPTWRK 737 WORK(2) = MINWRK 738 RWORK(1) = RMINWRK 739 RETURN 740 END IF 741* 742* Quick return if the matrix is void. 743* 744 IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN 745* .. all output is void. 746 RETURN 747 END IF 748* 749 BIG = SLAMCH('O') 750 ASCALED = .FALSE. 751 IWOFF = 1 752 IF ( ROWPRM ) THEN 753 IWOFF = M 754* .. reordering the rows in decreasing sequence in the 755* ell-infinity norm - this enhances numerical robustness in 756* the case of differently scaled rows. 757 DO 1904 p = 1, M 758* RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) ) 759* [[SLANGE will return NaN if an entry of the p-th row is Nan]] 760 RWORK(p) = SLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY ) 761* .. check for NaN's and Inf's 762 IF ( ( RWORK(p) .NE. RWORK(p) ) .OR. 763 $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN 764 INFO = -8 765 CALL XERBLA( 'SGESVDQ', -INFO ) 766 RETURN 767 END IF 768 1904 CONTINUE 769 DO 1952 p = 1, M - 1 770 q = ISAMAX( M-p+1, RWORK(p), 1 ) + p - 1 771 IWORK(N+p) = q 772 IF ( p .NE. q ) THEN 773 RTMP = RWORK(p) 774 RWORK(p) = RWORK(q) 775 RWORK(q) = RTMP 776 END IF 777 1952 CONTINUE 778* 779 IF ( RWORK(1) .EQ. ZERO ) THEN 780* Quick return: A is the M x N zero matrix. 781 NUMRANK = 0 782 CALL SLASET( 'G', N, 1, ZERO, ZERO, S, N ) 783 IF ( WNTUS ) CALL SLASET('G', M, N, ZERO, ONE, U, LDU) 784 IF ( WNTUA ) CALL SLASET('G', M, M, ZERO, ONE, U, LDU) 785 IF ( WNTVA ) CALL SLASET('G', N, N, ZERO, ONE, V, LDV) 786 IF ( WNTUF ) THEN 787 CALL SLASET( 'G', N, 1, ZERO, ZERO, WORK, N ) 788 CALL SLASET( 'G', M, N, ZERO, ONE, U, LDU ) 789 END IF 790 DO 5001 p = 1, N 791 IWORK(p) = p 792 5001 CONTINUE 793 IF ( ROWPRM ) THEN 794 DO 5002 p = N + 1, N + M - 1 795 IWORK(p) = p - N 796 5002 CONTINUE 797 END IF 798 IF ( CONDA ) RWORK(1) = -1 799 RWORK(2) = -1 800 RETURN 801 END IF 802* 803 IF ( RWORK(1) .GT. BIG / SQRT(REAL(M)) ) THEN 804* .. to prevent overflow in the QR factorization, scale the 805* matrix by 1/sqrt(M) if too large entry detected 806 CALL SLASCL('G',0,0,SQRT(REAL(M)),ONE, M,N, A,LDA, IERR) 807 ASCALED = .TRUE. 808 END IF 809 CALL SLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 ) 810 END IF 811* 812* .. At this stage, preemptive scaling is done only to avoid column 813* norms overflows during the QR factorization. The SVD procedure should 814* have its own scaling to save the singular values from overflows and 815* underflows. That depends on the SVD procedure. 816* 817 IF ( .NOT.ROWPRM ) THEN 818 RTMP = SLANGE( 'M', M, N, A, LDA, RDUMMY ) 819 IF ( ( RTMP .NE. RTMP ) .OR. 820 $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN 821 INFO = -8 822 CALL XERBLA( 'SGESVDQ', -INFO ) 823 RETURN 824 END IF 825 IF ( RTMP .GT. BIG / SQRT(REAL(M)) ) THEN 826* .. to prevent overflow in the QR factorization, scale the 827* matrix by 1/sqrt(M) if too large entry detected 828 CALL SLASCL('G',0,0, SQRT(REAL(M)),ONE, M,N, A,LDA, IERR) 829 ASCALED = .TRUE. 830 END IF 831 END IF 832* 833* .. QR factorization with column pivoting 834* 835* A * P = Q * [ R ] 836* [ 0 ] 837* 838 DO 1963 p = 1, N 839* .. all columns are free columns 840 IWORK(p) = 0 841 1963 CONTINUE 842 CALL SGEQP3( M, N, A, LDA, IWORK, WORK, WORK(N+1), LWORK-N, 843 $ IERR ) 844* 845* If the user requested accuracy level allows truncation in the 846* computed upper triangular factor, the matrix R is examined and, 847* if possible, replaced with its leading upper trapezoidal part. 848* 849 EPSLN = SLAMCH('E') 850 SFMIN = SLAMCH('S') 851* SMALL = SFMIN / EPSLN 852 NR = N 853* 854 IF ( ACCLA ) THEN 855* 856* Standard absolute error bound suffices. All sigma_i with 857* sigma_i < N*EPS*||A||_F are flushed to zero. This is an 858* aggressive enforcement of lower numerical rank by introducing a 859* backward error of the order of N*EPS*||A||_F. 860 NR = 1 861 RTMP = SQRT(REAL(N))*EPSLN 862 DO 3001 p = 2, N 863 IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002 864 NR = NR + 1 865 3001 CONTINUE 866 3002 CONTINUE 867* 868 ELSEIF ( ACCLM ) THEN 869* .. similarly as above, only slightly more gentle (less aggressive). 870* Sudden drop on the diagonal of R is used as the criterion for being 871* close-to-rank-deficient. The threshold is set to EPSLN=SLAMCH('E'). 872* [[This can be made more flexible by replacing this hard-coded value 873* with a user specified threshold.]] Also, the values that underflow 874* will be truncated. 875 NR = 1 876 DO 3401 p = 2, N 877 IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR. 878 $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402 879 NR = NR + 1 880 3401 CONTINUE 881 3402 CONTINUE 882* 883 ELSE 884* .. RRQR not authorized to determine numerical rank except in the 885* obvious case of zero pivots. 886* .. inspect R for exact zeros on the diagonal; 887* R(i,i)=0 => R(i:N,i:N)=0. 888 NR = 1 889 DO 3501 p = 2, N 890 IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502 891 NR = NR + 1 892 3501 CONTINUE 893 3502 CONTINUE 894* 895 IF ( CONDA ) THEN 896* Estimate the scaled condition number of A. Use the fact that it is 897* the same as the scaled condition number of R. 898* .. V is used as workspace 899 CALL SLACPY( 'U', N, N, A, LDA, V, LDV ) 900* Only the leading NR x NR submatrix of the triangular factor 901* is considered. Only if NR=N will this give a reliable error 902* bound. However, even for NR < N, this can be used on an 903* expert level and obtain useful information in the sense of 904* perturbation theory. 905 DO 3053 p = 1, NR 906 RTMP = SNRM2( p, V(1,p), 1 ) 907 CALL SSCAL( p, ONE/RTMP, V(1,p), 1 ) 908 3053 CONTINUE 909 IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN 910 CALL SPOCON( 'U', NR, V, LDV, ONE, RTMP, 911 $ WORK, IWORK(N+IWOFF), IERR ) 912 ELSE 913 CALL SPOCON( 'U', NR, V, LDV, ONE, RTMP, 914 $ WORK(N+1), IWORK(N+IWOFF), IERR ) 915 END IF 916 SCONDA = ONE / SQRT(RTMP) 917* For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1), 918* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA 919* See the reference [1] for more details. 920 END IF 921* 922 ENDIF 923* 924 IF ( WNTUR ) THEN 925 N1 = NR 926 ELSE IF ( WNTUS .OR. WNTUF) THEN 927 N1 = N 928 ELSE IF ( WNTUA ) THEN 929 N1 = M 930 END IF 931* 932 IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN 933*....................................................................... 934* .. only the singular values are requested 935*....................................................................... 936 IF ( RTRANS ) THEN 937* 938* .. compute the singular values of R**T = [A](1:NR,1:N)**T 939* .. set the lower triangle of [A] to [A](1:NR,1:N)**T and 940* the upper triangle of [A] to zero. 941 DO 1146 p = 1, MIN( N, NR ) 942 DO 1147 q = p + 1, N 943 A(q,p) = A(p,q) 944 IF ( q .LE. NR ) A(p,q) = ZERO 945 1147 CONTINUE 946 1146 CONTINUE 947* 948 CALL SGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU, 949 $ V, LDV, WORK, LWORK, INFO ) 950* 951 ELSE 952* 953* .. compute the singular values of R = [A](1:NR,1:N) 954* 955 IF ( NR .GT. 1 ) 956 $ CALL SLASET( 'L', NR-1,NR-1, ZERO,ZERO, A(2,1), LDA ) 957 CALL SGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU, 958 $ V, LDV, WORK, LWORK, INFO ) 959* 960 END IF 961* 962 ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN 963*....................................................................... 964* .. the singular values and the left singular vectors requested 965*......................................................................."""""""" 966 IF ( RTRANS ) THEN 967* .. apply SGESVD to R**T 968* .. copy R**T into [U] and overwrite [U] with the right singular 969* vectors of R 970 DO 1192 p = 1, NR 971 DO 1193 q = p, N 972 U(q,p) = A(p,q) 973 1193 CONTINUE 974 1192 CONTINUE 975 IF ( NR .GT. 1 ) 976 $ CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, U(1,2), LDU ) 977* .. the left singular vectors not computed, the NR right singular 978* vectors overwrite [U](1:NR,1:NR) as transposed. These 979* will be pre-multiplied by Q to build the left singular vectors of A. 980 CALL SGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU, 981 $ U, LDU, WORK(N+1), LWORK-N, INFO ) 982* 983 DO 1119 p = 1, NR 984 DO 1120 q = p + 1, NR 985 RTMP = U(q,p) 986 U(q,p) = U(p,q) 987 U(p,q) = RTMP 988 1120 CONTINUE 989 1119 CONTINUE 990* 991 ELSE 992* .. apply SGESVD to R 993* .. copy R into [U] and overwrite [U] with the left singular vectors 994 CALL SLACPY( 'U', NR, N, A, LDA, U, LDU ) 995 IF ( NR .GT. 1 ) 996 $ CALL SLASET( 'L', NR-1, NR-1, ZERO, ZERO, U(2,1), LDU ) 997* .. the right singular vectors not computed, the NR left singular 998* vectors overwrite [U](1:NR,1:NR) 999 CALL SGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU, 1000 $ V, LDV, WORK(N+1), LWORK-N, INFO ) 1001* .. now [U](1:NR,1:NR) contains the NR left singular vectors of 1002* R. These will be pre-multiplied by Q to build the left singular 1003* vectors of A. 1004 END IF 1005* 1006* .. assemble the left singular vector matrix U of dimensions 1007* (M x NR) or (M x N) or (M x M). 1008 IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) THEN 1009 CALL SLASET('A', M-NR, NR, ZERO, ZERO, U(NR+1,1), LDU) 1010 IF ( NR .LT. N1 ) THEN 1011 CALL SLASET( 'A',NR,N1-NR,ZERO,ZERO,U(1,NR+1), LDU ) 1012 CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE, 1013 $ U(NR+1,NR+1), LDU ) 1014 END IF 1015 END IF 1016* 1017* The Q matrix from the first QRF is built into the left singular 1018* vectors matrix U. 1019* 1020 IF ( .NOT.WNTUF ) 1021 $ CALL SORMQR( 'L', 'N', M, N1, N, A, LDA, WORK, U, 1022 $ LDU, WORK(N+1), LWORK-N, IERR ) 1023 IF ( ROWPRM .AND. .NOT.WNTUF ) 1024 $ CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 ) 1025* 1026 ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN 1027*....................................................................... 1028* .. the singular values and the right singular vectors requested 1029*....................................................................... 1030 IF ( RTRANS ) THEN 1031* .. apply SGESVD to R**T 1032* .. copy R**T into V and overwrite V with the left singular vectors 1033 DO 1165 p = 1, NR 1034 DO 1166 q = p, N 1035 V(q,p) = (A(p,q)) 1036 1166 CONTINUE 1037 1165 CONTINUE 1038 IF ( NR .GT. 1 ) 1039 $ CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV ) 1040* .. the left singular vectors of R**T overwrite V, the right singular 1041* vectors not computed 1042 IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN 1043 CALL SGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU, 1044 $ U, LDU, WORK(N+1), LWORK-N, INFO ) 1045* 1046 DO 1121 p = 1, NR 1047 DO 1122 q = p + 1, NR 1048 RTMP = V(q,p) 1049 V(q,p) = V(p,q) 1050 V(p,q) = RTMP 1051 1122 CONTINUE 1052 1121 CONTINUE 1053* 1054 IF ( NR .LT. N ) THEN 1055 DO 1103 p = 1, NR 1056 DO 1104 q = NR + 1, N 1057 V(p,q) = V(q,p) 1058 1104 CONTINUE 1059 1103 CONTINUE 1060 END IF 1061 CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK ) 1062 ELSE 1063* .. need all N right singular vectors and NR < N 1064* [!] This is simple implementation that augments [V](1:N,1:NR) 1065* by padding a zero block. In the case NR << N, a more efficient 1066* way is to first use the QR factorization. For more details 1067* how to implement this, see the " FULL SVD " branch. 1068 CALL SLASET('G', N, N-NR, ZERO, ZERO, V(1,NR+1), LDV) 1069 CALL SGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU, 1070 $ U, LDU, WORK(N+1), LWORK-N, INFO ) 1071* 1072 DO 1123 p = 1, N 1073 DO 1124 q = p + 1, N 1074 RTMP = V(q,p) 1075 V(q,p) = V(p,q) 1076 V(p,q) = RTMP 1077 1124 CONTINUE 1078 1123 CONTINUE 1079 CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK ) 1080 END IF 1081* 1082 ELSE 1083* .. aply SGESVD to R 1084* .. copy R into V and overwrite V with the right singular vectors 1085 CALL SLACPY( 'U', NR, N, A, LDA, V, LDV ) 1086 IF ( NR .GT. 1 ) 1087 $ CALL SLASET( 'L', NR-1, NR-1, ZERO, ZERO, V(2,1), LDV ) 1088* .. the right singular vectors overwrite V, the NR left singular 1089* vectors stored in U(1:NR,1:NR) 1090 IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN 1091 CALL SGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU, 1092 $ V, LDV, WORK(N+1), LWORK-N, INFO ) 1093 CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK ) 1094* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T 1095 ELSE 1096* .. need all N right singular vectors and NR < N 1097* [!] This is simple implementation that augments [V](1:NR,1:N) 1098* by padding a zero block. In the case NR << N, a more efficient 1099* way is to first use the LQ factorization. For more details 1100* how to implement this, see the " FULL SVD " branch. 1101 CALL SLASET('G', N-NR, N, ZERO,ZERO, V(NR+1,1), LDV) 1102 CALL SGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU, 1103 $ V, LDV, WORK(N+1), LWORK-N, INFO ) 1104 CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK ) 1105 END IF 1106* .. now [V] contains the transposed matrix of the right singular 1107* vectors of A. 1108 END IF 1109* 1110 ELSE 1111*....................................................................... 1112* .. FULL SVD requested 1113*....................................................................... 1114 IF ( RTRANS ) THEN 1115* 1116* .. apply SGESVD to R**T [[this option is left for R&D&T]] 1117* 1118 IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN 1119* .. copy R**T into [V] and overwrite [V] with the left singular 1120* vectors of R**T 1121 DO 1168 p = 1, NR 1122 DO 1169 q = p, N 1123 V(q,p) = A(p,q) 1124 1169 CONTINUE 1125 1168 CONTINUE 1126 IF ( NR .GT. 1 ) 1127 $ CALL SLASET( 'U', NR-1,NR-1, ZERO,ZERO, V(1,2), LDV ) 1128* 1129* .. the left singular vectors of R**T overwrite [V], the NR right 1130* singular vectors of R**T stored in [U](1:NR,1:NR) as transposed 1131 CALL SGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV, 1132 $ U, LDU, WORK(N+1), LWORK-N, INFO ) 1133* .. assemble V 1134 DO 1115 p = 1, NR 1135 DO 1116 q = p + 1, NR 1136 RTMP = V(q,p) 1137 V(q,p) = V(p,q) 1138 V(p,q) = RTMP 1139 1116 CONTINUE 1140 1115 CONTINUE 1141 IF ( NR .LT. N ) THEN 1142 DO 1101 p = 1, NR 1143 DO 1102 q = NR+1, N 1144 V(p,q) = V(q,p) 1145 1102 CONTINUE 1146 1101 CONTINUE 1147 END IF 1148 CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK ) 1149* 1150 DO 1117 p = 1, NR 1151 DO 1118 q = p + 1, NR 1152 RTMP = U(q,p) 1153 U(q,p) = U(p,q) 1154 U(p,q) = RTMP 1155 1118 CONTINUE 1156 1117 CONTINUE 1157* 1158 IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN 1159 CALL SLASET('A', M-NR,NR, ZERO,ZERO, U(NR+1,1), LDU) 1160 IF ( NR .LT. N1 ) THEN 1161 CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU) 1162 CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE, 1163 $ U(NR+1,NR+1), LDU ) 1164 END IF 1165 END IF 1166* 1167 ELSE 1168* .. need all N right singular vectors and NR < N 1169* .. copy R**T into [V] and overwrite [V] with the left singular 1170* vectors of R**T 1171* [[The optimal ratio N/NR for using QRF instead of padding 1172* with zeros. Here hard coded to 2; it must be at least 1173* two due to work space constraints.]] 1174* OPTRATIO = ILAENV(6, 'SGESVD', 'S' // 'O', NR,N,0,0) 1175* OPTRATIO = MAX( OPTRATIO, 2 ) 1176 OPTRATIO = 2 1177 IF ( OPTRATIO*NR .GT. N ) THEN 1178 DO 1198 p = 1, NR 1179 DO 1199 q = p, N 1180 V(q,p) = A(p,q) 1181 1199 CONTINUE 1182 1198 CONTINUE 1183 IF ( NR .GT. 1 ) 1184 $ CALL SLASET('U',NR-1,NR-1, ZERO,ZERO, V(1,2),LDV) 1185* 1186 CALL SLASET('A',N,N-NR,ZERO,ZERO,V(1,NR+1),LDV) 1187 CALL SGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV, 1188 $ U, LDU, WORK(N+1), LWORK-N, INFO ) 1189* 1190 DO 1113 p = 1, N 1191 DO 1114 q = p + 1, N 1192 RTMP = V(q,p) 1193 V(q,p) = V(p,q) 1194 V(p,q) = RTMP 1195 1114 CONTINUE 1196 1113 CONTINUE 1197 CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK ) 1198* .. assemble the left singular vector matrix U of dimensions 1199* (M x N1), i.e. (M x N) or (M x M). 1200* 1201 DO 1111 p = 1, N 1202 DO 1112 q = p + 1, N 1203 RTMP = U(q,p) 1204 U(q,p) = U(p,q) 1205 U(p,q) = RTMP 1206 1112 CONTINUE 1207 1111 CONTINUE 1208* 1209 IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN 1210 CALL SLASET('A',M-N,N,ZERO,ZERO,U(N+1,1),LDU) 1211 IF ( N .LT. N1 ) THEN 1212 CALL SLASET('A',N,N1-N,ZERO,ZERO,U(1,N+1),LDU) 1213 CALL SLASET('A',M-N,N1-N,ZERO,ONE, 1214 $ U(N+1,N+1), LDU ) 1215 END IF 1216 END IF 1217 ELSE 1218* .. copy R**T into [U] and overwrite [U] with the right 1219* singular vectors of R 1220 DO 1196 p = 1, NR 1221 DO 1197 q = p, N 1222 U(q,NR+p) = A(p,q) 1223 1197 CONTINUE 1224 1196 CONTINUE 1225 IF ( NR .GT. 1 ) 1226 $ CALL SLASET('U',NR-1,NR-1,ZERO,ZERO,U(1,NR+2),LDU) 1227 CALL SGEQRF( N, NR, U(1,NR+1), LDU, WORK(N+1), 1228 $ WORK(N+NR+1), LWORK-N-NR, IERR ) 1229 DO 1143 p = 1, NR 1230 DO 1144 q = 1, N 1231 V(q,p) = U(p,NR+q) 1232 1144 CONTINUE 1233 1143 CONTINUE 1234 CALL SLASET('U',NR-1,NR-1,ZERO,ZERO,V(1,2),LDV) 1235 CALL SGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU, 1236 $ V,LDV, WORK(N+NR+1),LWORK-N-NR, INFO ) 1237 CALL SLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV) 1238 CALL SLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV) 1239 CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV) 1240 CALL SORMQR('R','C', N, N, NR, U(1,NR+1), LDU, 1241 $ WORK(N+1),V,LDV,WORK(N+NR+1),LWORK-N-NR,IERR) 1242 CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK ) 1243* .. assemble the left singular vector matrix U of dimensions 1244* (M x NR) or (M x N) or (M x M). 1245 IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN 1246 CALL SLASET('A',M-NR,NR,ZERO,ZERO,U(NR+1,1),LDU) 1247 IF ( NR .LT. N1 ) THEN 1248 CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU) 1249 CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE, 1250 $ U(NR+1,NR+1),LDU) 1251 END IF 1252 END IF 1253 END IF 1254 END IF 1255* 1256 ELSE 1257* 1258* .. apply SGESVD to R [[this is the recommended option]] 1259* 1260 IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN 1261* .. copy R into [V] and overwrite V with the right singular vectors 1262 CALL SLACPY( 'U', NR, N, A, LDA, V, LDV ) 1263 IF ( NR .GT. 1 ) 1264 $ CALL SLASET( 'L', NR-1,NR-1, ZERO,ZERO, V(2,1), LDV ) 1265* .. the right singular vectors of R overwrite [V], the NR left 1266* singular vectors of R stored in [U](1:NR,1:NR) 1267 CALL SGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU, 1268 $ V, LDV, WORK(N+1), LWORK-N, INFO ) 1269 CALL SLAPMT( .FALSE., NR, N, V, LDV, IWORK ) 1270* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T 1271* .. assemble the left singular vector matrix U of dimensions 1272* (M x NR) or (M x N) or (M x M). 1273 IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN 1274 CALL SLASET('A', M-NR,NR, ZERO,ZERO, U(NR+1,1), LDU) 1275 IF ( NR .LT. N1 ) THEN 1276 CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU) 1277 CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE, 1278 $ U(NR+1,NR+1), LDU ) 1279 END IF 1280 END IF 1281* 1282 ELSE 1283* .. need all N right singular vectors and NR < N 1284* .. the requested number of the left singular vectors 1285* is then N1 (N or M) 1286* [[The optimal ratio N/NR for using LQ instead of padding 1287* with zeros. Here hard coded to 2; it must be at least 1288* two due to work space constraints.]] 1289* OPTRATIO = ILAENV(6, 'SGESVD', 'S' // 'O', NR,N,0,0) 1290* OPTRATIO = MAX( OPTRATIO, 2 ) 1291 OPTRATIO = 2 1292 IF ( OPTRATIO * NR .GT. N ) THEN 1293 CALL SLACPY( 'U', NR, N, A, LDA, V, LDV ) 1294 IF ( NR .GT. 1 ) 1295 $ CALL SLASET('L', NR-1,NR-1, ZERO,ZERO, V(2,1),LDV) 1296* .. the right singular vectors of R overwrite [V], the NR left 1297* singular vectors of R stored in [U](1:NR,1:NR) 1298 CALL SLASET('A', N-NR,N, ZERO,ZERO, V(NR+1,1),LDV) 1299 CALL SGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU, 1300 $ V, LDV, WORK(N+1), LWORK-N, INFO ) 1301 CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK ) 1302* .. now [V] contains the transposed matrix of the right 1303* singular vectors of A. The leading N left singular vectors 1304* are in [U](1:N,1:N) 1305* .. assemble the left singular vector matrix U of dimensions 1306* (M x N1), i.e. (M x N) or (M x M). 1307 IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN 1308 CALL SLASET('A',M-N,N,ZERO,ZERO,U(N+1,1),LDU) 1309 IF ( N .LT. N1 ) THEN 1310 CALL SLASET('A',N,N1-N,ZERO,ZERO,U(1,N+1),LDU) 1311 CALL SLASET( 'A',M-N,N1-N,ZERO,ONE, 1312 $ U(N+1,N+1), LDU ) 1313 END IF 1314 END IF 1315 ELSE 1316 CALL SLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU ) 1317 IF ( NR .GT. 1 ) 1318 $ CALL SLASET('L',NR-1,NR-1,ZERO,ZERO,U(NR+2,1),LDU) 1319 CALL SGELQF( NR, N, U(NR+1,1), LDU, WORK(N+1), 1320 $ WORK(N+NR+1), LWORK-N-NR, IERR ) 1321 CALL SLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV) 1322 IF ( NR .GT. 1 ) 1323 $ CALL SLASET('U',NR-1,NR-1,ZERO,ZERO,V(1,2),LDV) 1324 CALL SGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU, 1325 $ V, LDV, WORK(N+NR+1), LWORK-N-NR, INFO ) 1326 CALL SLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV) 1327 CALL SLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV) 1328 CALL SLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV) 1329 CALL SORMLQ('R','N',N,N,NR,U(NR+1,1),LDU,WORK(N+1), 1330 $ V, LDV, WORK(N+NR+1),LWORK-N-NR,IERR) 1331 CALL SLAPMT( .FALSE., N, N, V, LDV, IWORK ) 1332* .. assemble the left singular vector matrix U of dimensions 1333* (M x NR) or (M x N) or (M x M). 1334 IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN 1335 CALL SLASET('A',M-NR,NR,ZERO,ZERO,U(NR+1,1),LDU) 1336 IF ( NR .LT. N1 ) THEN 1337 CALL SLASET('A',NR,N1-NR,ZERO,ZERO,U(1,NR+1),LDU) 1338 CALL SLASET( 'A',M-NR,N1-NR,ZERO,ONE, 1339 $ U(NR+1,NR+1), LDU ) 1340 END IF 1341 END IF 1342 END IF 1343 END IF 1344* .. end of the "R**T or R" branch 1345 END IF 1346* 1347* The Q matrix from the first QRF is built into the left singular 1348* vectors matrix U. 1349* 1350 IF ( .NOT. WNTUF ) 1351 $ CALL SORMQR( 'L', 'N', M, N1, N, A, LDA, WORK, U, 1352 $ LDU, WORK(N+1), LWORK-N, IERR ) 1353 IF ( ROWPRM .AND. .NOT.WNTUF ) 1354 $ CALL SLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 ) 1355* 1356* ... end of the "full SVD" branch 1357 END IF 1358* 1359* Check whether some singular values are returned as zeros, e.g. 1360* due to underflow, and update the numerical rank. 1361 p = NR 1362 DO 4001 q = p, 1, -1 1363 IF ( S(q) .GT. ZERO ) GO TO 4002 1364 NR = NR - 1 1365 4001 CONTINUE 1366 4002 CONTINUE 1367* 1368* .. if numerical rank deficiency is detected, the truncated 1369* singular values are set to zero. 1370 IF ( NR .LT. N ) CALL SLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N ) 1371* .. undo scaling; this may cause overflow in the largest singular 1372* values. 1373 IF ( ASCALED ) 1374 $ CALL SLASCL( 'G',0,0, ONE,SQRT(REAL(M)), NR,1, S, N, IERR ) 1375 IF ( CONDA ) RWORK(1) = SCONDA 1376 RWORK(2) = p - NR 1377* .. p-NR is the number of singular values that are computed as 1378* exact zeros in SGESVD() applied to the (possibly truncated) 1379* full row rank triangular (trapezoidal) factor of A. 1380 NUMRANK = NR 1381* 1382 RETURN 1383* 1384* End of SGESVDQ 1385* 1386 END 1387