LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sgesvdq()

subroutine sgesvdq ( character joba,
character jobp,
character jobr,
character jobu,
character jobv,
integer m,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) s,
real, dimension( ldu, * ) u,
integer ldu,
real, dimension( ldv, * ) v,
integer ldv,
integer numrank,
integer, dimension( * ) iwork,
integer liwork,
real, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
integer lrwork,
integer info )

SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices

Download SGESVDQ + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SGESVDQ computes the singular value decomposition (SVD) of a real
!> M-by-N matrix A, where M >= N. The SVD of A is written as
!>                                    [++]   [xx]   [x0]   [xx]
!>              A = U * SIGMA * V^*,  [++] = [xx] * [ox] * [xx]
!>                                    [++]   [xx]
!> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
!> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
!> of SIGMA are the singular values of A. The columns of U and V are the
!> left and the right singular vectors of A, respectively.
!>
Parameters
[in]JOBA
!>  JOBA is CHARACTER*1
!>  Specifies the level of accuracy in the computed SVD
!>  = 'A' The requested accuracy corresponds to having the backward
!>        error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F,
!>        where EPS = SLAMCH('Epsilon'). This authorises CGESVDQ to
!>        truncate the computed triangular factor in a rank revealing
!>        QR factorization whenever the truncated part is below the
!>        threshold of the order of EPS * ||A||_F. This is aggressive
!>        truncation level.
!>  = 'M' Similarly as with 'A', but the truncation is more gentle: it
!>        is allowed only when there is a drop on the diagonal of the
!>        triangular factor in the QR factorization. This is medium
!>        truncation level.
!>  = 'H' High accuracy requested. No numerical rank determination based
!>        on the rank revealing QR factorization is attempted.
!>  = 'E' Same as 'H', and in addition the condition number of column
!>        scaled A is estimated and returned in  RWORK(1).
!>        N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
!>
[in]JOBP
!>  JOBP is CHARACTER*1
!>  = 'P' The rows of A are ordered in decreasing order with respect to
!>        ||A(i,:)||_\infty. This enhances numerical accuracy at the cost
!>        of extra data movement. Recommended for numerical robustness.
!>  = 'N' No row pivoting.
!>
[in]JOBR
!>          JOBR is CHARACTER*1
!>          = 'T' After the initial pivoted QR factorization, SGESVD is applied to
!>          the transposed R**T of the computed triangular factor R. This involves
!>          some extra data movement (matrix transpositions). Useful for
!>          experiments, research and development.
!>          = 'N' The triangular factor R is given as input to SGESVD. This may be
!>          preferred as it involves less data movement.
!>
[in]JOBU
!>          JOBU is CHARACTER*1
!>          = 'A' All M left singular vectors are computed and returned in the
!>          matrix U. See the description of U.
!>          = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned
!>          in the matrix U. See the description of U.
!>          = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular
!>          vectors are computed and returned in the matrix U.
!>          = 'F' The N left singular vectors are returned in factored form as the
!>          product of the Q factor from the initial QR factorization and the
!>          N left singular vectors of (R**T , 0)**T. If row pivoting is used,
!>          then the necessary information on the row pivoting is stored in
!>          IWORK(N+1:N+M-1).
!>          = 'N' The left singular vectors are not computed.
!>
[in]JOBV
!>          JOBV is CHARACTER*1
!>          = 'A', 'V' All N right singular vectors are computed and returned in
!>          the matrix V.
!>          = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular
!>          vectors are computed and returned in the matrix V. This option is
!>          allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal.
!>          = 'N' The right singular vectors are not computed.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the input matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the input matrix A.  M >= N >= 0.
!>
[in,out]A
!>          A is REAL array of dimensions LDA x N
!>          On entry, the input matrix A.
!>          On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains
!>          the Householder vectors as stored by SGEQP3. If JOBU = 'F', these Householder
!>          vectors together with WORK(1:N) can be used to restore the Q factors from
!>          the initial pivoted QR factorization of A. See the description of U.
!>
[in]LDA
!>          LDA is INTEGER.
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[out]S
!>          S is REAL array of dimension N.
!>          The singular values of A, ordered so that S(i) >= S(i+1).
!>
[out]U
!>          U is REAL array, dimension
!>          LDU x M if JOBU = 'A'; see the description of LDU. In this case,
!>          on exit, U contains the M left singular vectors.
!>          LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this
!>          case, U contains the leading N or the leading NUMRANK left singular vectors.
!>          LDU x N if JOBU = 'F' ; see the description of LDU. In this case U
!>          contains N x N orthogonal matrix that can be used to form the left
!>          singular vectors.
!>          If JOBU = 'N', U is not referenced.
!>
[in]LDU
!>          LDU is INTEGER.
!>          The leading dimension of the array U.
!>          If JOBU = 'A', 'S', 'U', 'R',  LDU >= max(1,M).
!>          If JOBU = 'F',                 LDU >= max(1,N).
!>          Otherwise,                     LDU >= 1.
!>
[out]V
!>          V is REAL array, dimension
!>          LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' .
!>          If JOBV = 'A', or 'V',  V contains the N-by-N orthogonal matrix  V**T;
!>          If JOBV = 'R', V contains the first NUMRANK rows of V**T (the right
!>          singular vectors, stored rowwise, of the NUMRANK largest singular values).
!>          If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
!>          If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V.
!>          If JOBV = 'A', 'V', 'R',  or JOBA = 'E', LDV >= max(1,N).
!>          Otherwise,                               LDV >= 1.
!>
[out]NUMRANK
!>          NUMRANK is INTEGER
!>          NUMRANK is the numerical rank first determined after the rank
!>          revealing QR factorization, following the strategy specified by the
!>          value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK
!>          leading singular values and vectors are then requested in the call
!>          of SGESVD. The final value of NUMRANK might be further reduced if
!>          some singular values are computed as zeros.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (max(1, LIWORK)).
!>          On exit, IWORK(1:N) contains column pivoting permutation of the
!>          rank revealing QR factorization.
!>          If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence
!>          of row swaps used in row pivoting. These can be used to restore the
!>          left singular vectors in the case JOBU = 'F'.
!>
!>          If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
!>          IWORK(1) returns the minimal LIWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.
!>          LIWORK >= N + M - 1,     if JOBP = 'P' and JOBA .NE. 'E';
!>          LIWORK >= N              if JOBP = 'N' and JOBA .NE. 'E';
!>          LIWORK >= N + M - 1 + N, if JOBP = 'P' and JOBA = 'E';
!>          LIWORK >= N + N          if JOBP = 'N' and JOBA = 'E'.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the routine
!>          only calculates and returns the optimal and minimal sizes
!>          for the WORK, IWORK, and RWORK arrays, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]WORK
!>          WORK is REAL array, dimension (max(2, LWORK)), used as a workspace.
!>          On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains parameters
!>          needed to recover the Q factor from the QR factorization computed by
!>          SGEQP3.
!>
!>          If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
!>          WORK(1) returns the optimal LWORK, and
!>          WORK(2) returns the minimal LWORK.
!>
[in,out]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. It is determined as follows:
!>          Let  LWQP3 = 3*N+1,  LWCON = 3*N, and let
!>          LWORQ = { MAX( N, 1 ),  if JOBU = 'R', 'S', or 'U'
!>                  { MAX( M, 1 ),  if JOBU = 'A'
!>          LWSVD = MAX( 5*N, 1 )
!>          LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ = MAX( N, 1 ),
!>          LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
!>          Then the minimal value of LWORK is:
!>          = MAX( N + LWQP3, LWSVD )        if only the singular values are needed;
!>          = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
!>                                   and a scaled condition estimate requested;
!>
!>          = N + MAX( LWQP3, LWSVD, LWORQ ) if the singular values and the left
!>                                   singular vectors are requested;
!>          = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular values and the left
!>                                   singular vectors are requested, and also
!>                                   a scaled condition estimate requested;
!>
!>          = N + MAX( LWQP3, LWSVD )        if the singular values and the right
!>                                   singular vectors are requested;
!>          = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
!>                                   singular vectors are requested, and also
!>                                   a scaled condition etimate requested;
!>
!>          = N + MAX( LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV = 'R';
!>                                   independent of JOBR;
!>          = N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is requested,
!>                                   JOBV = 'R' and, also a scaled condition
!>                                   estimate requested; independent of JOBR;
!>          = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
!>         N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) ) if the
!>                         full SVD is requested with JOBV = 'A' or 'V', and
!>                         JOBR ='N'
!>          = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
!>         N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ ) )
!>                         if the full SVD is requested with JOBV = 'A' or 'V', and
!>                         JOBR ='N', and also a scaled condition number estimate
!>                         requested.
!>          = MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
!>         N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) ) if the
!>                         full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
!>          = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
!>         N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
!>                         if the full SVD is requested with JOBV = 'A' or 'V', and
!>                         JOBR ='T', and also a scaled condition number estimate
!>                         requested.
!>          Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK ).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates and returns the optimal and minimal sizes
!>          for the WORK, IWORK, and RWORK arrays, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (max(1, LRWORK)).
!>          On exit,
!>          1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
!>          number of column scaled A. If A = C * D where D is diagonal and C
!>          has unit columns in the Euclidean norm, then, assuming full column rank,
!>          N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
!>          Otherwise, RWORK(1) = -1.
!>          2. RWORK(2) contains the number of singular values computed as
!>          exact zeros in SGESVD applied to the upper triangular or trapezoidal
!>          R (from the initial QR factorization). In case of early exit (no call to
!>          SGESVD, such as in the case of zero matrix) RWORK(2) = -1.
!>
!>          If LIWORK, LWORK, or LRWORK = -1, then on exit, if INFO = 0,
!>          RWORK(1) returns the minimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER.
!>          The dimension of the array RWORK.
!>          If JOBP ='P', then LRWORK >= MAX(2, M).
!>          Otherwise, LRWORK >= 2
!>
!>          If LRWORK = -1, then a workspace query is assumed; the routine
!>          only calculates and returns the optimal and minimal sizes
!>          for the WORK, IWORK, and RWORK arrays, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if SBDSQR did not converge, INFO specifies how many superdiagonals
!>          of an intermediate bidiagonal form B (computed in SGESVD) did not
!>          converge to zero.
!>
Further Details:
!>
!>   1. The data movement (matrix transpose) is coded using simple nested
!>   DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
!>   Those DO-loops are easily identified in this source code - by the CONTINUE
!>   statements labeled with 11**. In an optimized version of this code, the
!>   nested DO loops should be replaced with calls to an optimized subroutine.
!>   2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
!>   column norm overflow. This is the minial precaution and it is left to the
!>   SVD routine (CGESVD) to do its own preemptive scaling if potential over-
!>   or underflows are detected. To avoid repeated scanning of the array A,
!>   an optimal implementation would do all necessary scaling before calling
!>   CGESVD and the scaling in CGESVD can be switched off.
!>   3. Other comments related to code optimization are given in comments in the
!>   code, enclosed in [[double brackets]].
!>
Bugs, examples and comments
!>  Please report all bugs and send interesting examples and/or comments to
!>  drmac@math.hr. Thank you.
!>
References
!>  [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
!>      Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
!>      44(1): 11:1-11:30 (2017)
!>
!>  SIGMA library, xGESVDQ section updated February 2016.
!>  Developed and coded by Zlatko Drmac, Department of Mathematics
!>  University of Zagreb, Croatia, drmac@math.hr
!>
Contributors:
!> Developed and coded by Zlatko Drmac, Department of Mathematics
!>  University of Zagreb, Croatia, drmac@math.hr
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 410 of file sgesvdq.f.

413* .. Scalar Arguments ..
414 IMPLICIT NONE
415 CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
416 INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LWORK, LRWORK,
417 $ INFO
418* ..
419* .. Array Arguments ..
420 REAL A( LDA, * ), U( LDU, * ), V( LDV, * ), WORK( * )
421 REAL S( * ), RWORK( * )
422 INTEGER IWORK( * )
423*
424* =====================================================================
425*
426* .. Parameters ..
427 REAL ZERO, ONE
428 parameter( zero = 0.0e0, one = 1.0e0 )
429* ..
430* .. Local Scalars ..
431 INTEGER IERR, IWOFF, NR, N1, OPTRATIO, p, q
432 INTEGER LWCON, LWQP3, LWRK_SGELQF, LWRK_SGESVD, LWRK_SGESVD2,
433 $ LWRK_SGEQP3, LWRK_SGEQRF, LWRK_SORMLQ, LWRK_SORMQR,
434 $ LWRK_SORMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWORQ,
435 $ LWORQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2,
436 $ IMINWRK, RMINWRK
437 LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
438 $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
439 $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR
440 REAL BIG, EPSLN, RTMP, SCONDA, SFMIN
441* ..
442* .. Local Arrays
443 REAL RDUMMY(1)
444* ..
445* .. External Subroutines (BLAS, LAPACK)
446 EXTERNAL sgelqf, sgeqp3, sgeqrf, sgesvd, slacpy,
447 $ slapmt, slascl, slaset, slaswp, sscal,
448 $ spocon, sormlq, sormqr, xerbla
449* ..
450* .. External Functions (BLAS, LAPACK)
451 LOGICAL LSAME
452 INTEGER ISAMAX
453 REAL SLANGE, SNRM2, SLAMCH
454 EXTERNAL slange, lsame, isamax, snrm2, slamch
455* ..
456* .. Intrinsic Functions ..
457 INTRINSIC abs, max, min, real, sqrt
458* ..
459* .. Executable Statements ..
460*
461* Test the input arguments
462*
463 wntus = lsame( jobu, 'S' ) .OR. lsame( jobu, 'U' )
464 wntur = lsame( jobu, 'R' )
465 wntua = lsame( jobu, 'A' )
466 wntuf = lsame( jobu, 'F' )
467 lsvc0 = wntus .OR. wntur .OR. wntua
468 lsvec = lsvc0 .OR. wntuf
469 dntwu = lsame( jobu, 'N' )
470*
471 wntvr = lsame( jobv, 'R' )
472 wntva = lsame( jobv, 'A' ) .OR. lsame( jobv, 'V' )
473 rsvec = wntvr .OR. wntva
474 dntwv = lsame( jobv, 'N' )
475*
476 accla = lsame( joba, 'A' )
477 acclm = lsame( joba, 'M' )
478 conda = lsame( joba, 'E' )
479 acclh = lsame( joba, 'H' ) .OR. conda
480*
481 rowprm = lsame( jobp, 'P' )
482 rtrans = lsame( jobr, 'T' )
483*
484 IF ( rowprm ) THEN
485 IF ( conda ) THEN
486 iminwrk = max( 1, n + m - 1 + n )
487 ELSE
488 iminwrk = max( 1, n + m - 1 )
489 END IF
490 rminwrk = max( 2, m )
491 ELSE
492 IF ( conda ) THEN
493 iminwrk = max( 1, n + n )
494 ELSE
495 iminwrk = max( 1, n )
496 END IF
497 rminwrk = 2
498 END IF
499 lquery = (liwork .EQ. -1 .OR. lwork .EQ. -1 .OR. lrwork .EQ. -1)
500 info = 0
501 IF ( .NOT. ( accla .OR. acclm .OR. acclh ) ) THEN
502 info = -1
503 ELSE IF ( .NOT.( rowprm .OR. lsame( jobp, 'N' ) ) ) THEN
504 info = -2
505 ELSE IF ( .NOT.( rtrans .OR. lsame( jobr, 'N' ) ) ) THEN
506 info = -3
507 ELSE IF ( .NOT.( lsvec .OR. dntwu ) ) THEN
508 info = -4
509 ELSE IF ( wntur .AND. wntva ) THEN
510 info = -5
511 ELSE IF ( .NOT.( rsvec .OR. dntwv )) THEN
512 info = -5
513 ELSE IF ( m.LT.0 ) THEN
514 info = -6
515 ELSE IF ( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
516 info = -7
517 ELSE IF ( lda.LT.max( 1, m ) ) THEN
518 info = -9
519 ELSE IF ( ldu.LT.1 .OR. ( lsvc0 .AND. ldu.LT.m ) .OR.
520 $ ( wntuf .AND. ldu.LT.n ) ) THEN
521 info = -12
522 ELSE IF ( ldv.LT.1 .OR. ( rsvec .AND. ldv.LT.n ) .OR.
523 $ ( conda .AND. ldv.LT.n ) ) THEN
524 info = -14
525 ELSE IF ( liwork .LT. iminwrk .AND. .NOT. lquery ) THEN
526 info = -17
527 END IF
528*
529*
530 IF ( info .EQ. 0 ) THEN
531* .. compute the minimal and the optimal workspace lengths
532* [[The expressions for computing the minimal and the optimal
533* values of LWORK are written with a lot of redundancy and
534* can be simplified. However, this detailed form is easier for
535* maintenance and modifications of the code.]]
536*
537* .. minimal workspace length for SGEQP3 of an M x N matrix
538 lwqp3 = 3 * n + 1
539* .. minimal workspace length for SORMQR to build left singular vectors
540 IF ( wntus .OR. wntur ) THEN
541 lworq = max( n , 1 )
542 ELSE IF ( wntua ) THEN
543 lworq = max( m , 1 )
544 END IF
545* .. minimal workspace length for SPOCON of an N x N matrix
546 lwcon = 3 * n
547* .. SGESVD of an N x N matrix
548 lwsvd = max( 5 * n, 1 )
549 IF ( lquery ) THEN
550 CALL sgeqp3( m, n, a, lda, iwork, rdummy, rdummy, -1,
551 $ ierr )
552 lwrk_sgeqp3 = int( rdummy(1) )
553 IF ( wntus .OR. wntur ) THEN
554 CALL sormqr( 'L', 'N', m, n, n, a, lda, rdummy, u,
555 $ ldu, rdummy, -1, ierr )
556 lwrk_sormqr = int( rdummy(1) )
557 ELSE IF ( wntua ) THEN
558 CALL sormqr( 'L', 'N', m, m, n, a, lda, rdummy, u,
559 $ ldu, rdummy, -1, ierr )
560 lwrk_sormqr = int( rdummy(1) )
561 ELSE
562 lwrk_sormqr = 0
563 END IF
564 END IF
565 minwrk = 2
566 optwrk = 2
567 IF ( .NOT. (lsvec .OR. rsvec )) THEN
568* .. minimal and optimal sizes of the workspace if
569* only the singular values are requested
570 IF ( conda ) THEN
571 minwrk = max( n+lwqp3, lwcon, lwsvd )
572 ELSE
573 minwrk = max( n+lwqp3, lwsvd )
574 END IF
575 IF ( lquery ) THEN
576 CALL sgesvd( 'N', 'N', n, n, a, lda, s, u, ldu,
577 $ v, ldv, rdummy, -1, ierr )
578 lwrk_sgesvd = int( rdummy(1) )
579 IF ( conda ) THEN
580 optwrk = max( n+lwrk_sgeqp3, n+lwcon, lwrk_sgesvd )
581 ELSE
582 optwrk = max( n+lwrk_sgeqp3, lwrk_sgesvd )
583 END IF
584 END IF
585 ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
586* .. minimal and optimal sizes of the workspace if the
587* singular values and the left singular vectors are requested
588 IF ( conda ) THEN
589 minwrk = n + max( lwqp3, lwcon, lwsvd, lworq )
590 ELSE
591 minwrk = n + max( lwqp3, lwsvd, lworq )
592 END IF
593 IF ( lquery ) THEN
594 IF ( rtrans ) THEN
595 CALL sgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
596 $ v, ldv, rdummy, -1, ierr )
597 ELSE
598 CALL sgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
599 $ v, ldv, rdummy, -1, ierr )
600 END IF
601 lwrk_sgesvd = int( rdummy(1) )
602 IF ( conda ) THEN
603 optwrk = n + max( lwrk_sgeqp3, lwcon, lwrk_sgesvd,
604 $ lwrk_sormqr )
605 ELSE
606 optwrk = n + max( lwrk_sgeqp3, lwrk_sgesvd,
607 $ lwrk_sormqr )
608 END IF
609 END IF
610 ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
611* .. minimal and optimal sizes of the workspace if the
612* singular values and the right singular vectors are requested
613 IF ( conda ) THEN
614 minwrk = n + max( lwqp3, lwcon, lwsvd )
615 ELSE
616 minwrk = n + max( lwqp3, lwsvd )
617 END IF
618 IF ( lquery ) THEN
619 IF ( rtrans ) THEN
620 CALL sgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
621 $ v, ldv, rdummy, -1, ierr )
622 ELSE
623 CALL sgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
624 $ v, ldv, rdummy, -1, ierr )
625 END IF
626 lwrk_sgesvd = int( rdummy(1) )
627 IF ( conda ) THEN
628 optwrk = n + max( lwrk_sgeqp3, lwcon, lwrk_sgesvd )
629 ELSE
630 optwrk = n + max( lwrk_sgeqp3, lwrk_sgesvd )
631 END IF
632 END IF
633 ELSE
634* .. minimal and optimal sizes of the workspace if the
635* full SVD is requested
636 IF ( rtrans ) THEN
637 minwrk = max( lwqp3, lwsvd, lworq )
638 IF ( conda ) minwrk = max( minwrk, lwcon )
639 minwrk = minwrk + n
640 IF ( wntva ) THEN
641* .. minimal workspace length for N x N/2 SGEQRF
642 lwqrf = max( n/2, 1 )
643* .. minimal workspace length for N/2 x N/2 SGESVD
644 lwsvd2 = max( 5 * (n/2), 1 )
645 lworq2 = max( n, 1 )
646 minwrk2 = max( lwqp3, n/2+lwqrf, n/2+lwsvd2,
647 $ n/2+lworq2, lworq )
648 IF ( conda ) minwrk2 = max( minwrk2, lwcon )
649 minwrk2 = n + minwrk2
650 minwrk = max( minwrk, minwrk2 )
651 END IF
652 ELSE
653 minwrk = max( lwqp3, lwsvd, lworq )
654 IF ( conda ) minwrk = max( minwrk, lwcon )
655 minwrk = minwrk + n
656 IF ( wntva ) THEN
657* .. minimal workspace length for N/2 x N SGELQF
658 lwlqf = max( n/2, 1 )
659 lwsvd2 = max( 5 * (n/2), 1 )
660 lwunlq = max( n , 1 )
661 minwrk2 = max( lwqp3, n/2+lwlqf, n/2+lwsvd2,
662 $ n/2+lwunlq, lworq )
663 IF ( conda ) minwrk2 = max( minwrk2, lwcon )
664 minwrk2 = n + minwrk2
665 minwrk = max( minwrk, minwrk2 )
666 END IF
667 END IF
668 IF ( lquery ) THEN
669 IF ( rtrans ) THEN
670 CALL sgesvd( 'O', 'A', n, n, a, lda, s, u, ldu,
671 $ v, ldv, rdummy, -1, ierr )
672 lwrk_sgesvd = int( rdummy(1) )
673 optwrk = max(lwrk_sgeqp3,lwrk_sgesvd,lwrk_sormqr)
674 IF ( conda ) optwrk = max( optwrk, lwcon )
675 optwrk = n + optwrk
676 IF ( wntva ) THEN
677 CALL sgeqrf(n,n/2,u,ldu,rdummy,rdummy,-1,ierr)
678 lwrk_sgeqrf = int( rdummy(1) )
679 CALL sgesvd( 'S', 'O', n/2,n/2, v,ldv, s, u,
680 $ ldu,
681 $ v, ldv, rdummy, -1, ierr )
682 lwrk_sgesvd2 = int( rdummy(1) )
683 CALL sormqr( 'R', 'C', n, n, n/2, u, ldu,
684 $ 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,
704 $ ldu,
705 $ v, ldv, rdummy, -1, ierr )
706 lwrk_sgesvd2 = int( rdummy(1) )
707 CALL sormlq( 'R', 'N', n, n, n/2, u, ldu,
708 $ rdummy,
709 $ v, ldv, rdummy,-1,ierr )
710 lwrk_sormlq = int( rdummy(1) )
711 optwrk2 = max( lwrk_sgeqp3, n/2+lwrk_sgelqf,
712 $ n/2+lwrk_sgesvd2, n/2+lwrk_sormlq )
713 IF ( conda ) optwrk2 = max( optwrk2, lwcon )
714 optwrk2 = n + optwrk2
715 optwrk = max( optwrk, optwrk2 )
716 END IF
717 END IF
718 END IF
719 END IF
720*
721 minwrk = max( 2, minwrk )
722 optwrk = max( 2, optwrk )
723 IF ( lwork .LT. minwrk .AND. (.NOT.lquery) ) info = -19
724*
725 END IF
726*
727 IF (info .EQ. 0 .AND. lrwork .LT. rminwrk .AND. .NOT. lquery) THEN
728 info = -21
729 END IF
730 IF( info.NE.0 ) THEN
731 CALL xerbla( 'SGESVDQ', -info )
732 RETURN
733 ELSE IF ( lquery ) THEN
734*
735* Return optimal workspace
736*
737 iwork(1) = iminwrk
738 work(1) = real( optwrk )
739 work(2) = real( minwrk )
740 rwork(1) = real( rminwrk )
741 RETURN
742 END IF
743*
744* Quick return if the matrix is void.
745*
746 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) ) THEN
747* .. all output is void.
748 RETURN
749 END IF
750*
751 big = slamch('O')
752 ascaled = .false.
753 iwoff = 1
754 IF ( rowprm ) THEN
755 iwoff = m
756* .. reordering the rows in decreasing sequence in the
757* ell-infinity norm - this enhances numerical robustness in
758* the case of differently scaled rows.
759 DO 1904 p = 1, m
760* RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
761* [[SLANGE will return NaN if an entry of the p-th row is Nan]]
762 rwork(p) = slange( 'M', 1, n, a(p,1), lda, rdummy )
763* .. check for NaN's and Inf's
764 IF ( ( rwork(p) .NE. rwork(p) ) .OR.
765 $ ( (rwork(p)*zero) .NE. zero ) ) THEN
766 info = -8
767 CALL xerbla( 'SGESVDQ', -info )
768 RETURN
769 END IF
770 1904 CONTINUE
771 DO 1952 p = 1, m - 1
772 q = isamax( m-p+1, rwork(p), 1 ) + p - 1
773 iwork(n+p) = q
774 IF ( p .NE. q ) THEN
775 rtmp = rwork(p)
776 rwork(p) = rwork(q)
777 rwork(q) = rtmp
778 END IF
779 1952 CONTINUE
780*
781 IF ( rwork(1) .EQ. zero ) THEN
782* Quick return: A is the M x N zero matrix.
783 numrank = 0
784 CALL slaset( 'G', n, 1, zero, zero, s, n )
785 IF ( wntus ) CALL slaset('G', m, n, zero, one, u, ldu)
786 IF ( wntua ) CALL slaset('G', m, m, zero, one, u, ldu)
787 IF ( wntva ) CALL slaset('G', n, n, zero, one, v, ldv)
788 IF ( wntuf ) THEN
789 CALL slaset( 'G', n, 1, zero, zero, work, n )
790 CALL slaset( 'G', m, n, zero, one, u, ldu )
791 END IF
792 DO 5001 p = 1, n
793 iwork(p) = p
794 5001 CONTINUE
795 IF ( rowprm ) THEN
796 DO 5002 p = n + 1, n + m - 1
797 iwork(p) = p - n
798 5002 CONTINUE
799 END IF
800 IF ( conda ) rwork(1) = -1
801 rwork(2) = -1
802 RETURN
803 END IF
804*
805 IF ( rwork(1) .GT. big / sqrt(real(m)) ) THEN
806* .. to prevent overflow in the QR factorization, scale the
807* matrix by 1/sqrt(M) if too large entry detected
808 CALL slascl('G',0,0,sqrt(real(m)),one, m,n, a,lda,
809 $ ierr)
810 ascaled = .true.
811 END IF
812 CALL slaswp( n, a, lda, 1, m-1, iwork(n+1), 1 )
813 END IF
814*
815* .. At this stage, preemptive scaling is done only to avoid column
816* norms overflows during the QR factorization. The SVD procedure should
817* have its own scaling to save the singular values from overflows and
818* underflows. That depends on the SVD procedure.
819*
820 IF ( .NOT.rowprm ) THEN
821 rtmp = slange( 'M', m, n, a, lda, rdummy )
822 IF ( ( rtmp .NE. rtmp ) .OR.
823 $ ( (rtmp*zero) .NE. zero ) ) THEN
824 info = -8
825 CALL xerbla( 'SGESVDQ', -info )
826 RETURN
827 END IF
828 IF ( rtmp .GT. big / sqrt(real(m)) ) THEN
829* .. to prevent overflow in the QR factorization, scale the
830* matrix by 1/sqrt(M) if too large entry detected
831 CALL slascl('G',0,0, sqrt(real(m)),one, m,n, a,lda,
832 $ ierr)
833 ascaled = .true.
834 END IF
835 END IF
836*
837* .. QR factorization with column pivoting
838*
839* A * P = Q * [ R ]
840* [ 0 ]
841*
842 DO 1963 p = 1, n
843* .. all columns are free columns
844 iwork(p) = 0
845 1963 CONTINUE
846 CALL sgeqp3( m, n, a, lda, iwork, work, work(n+1), lwork-n,
847 $ ierr )
848*
849* If the user requested accuracy level allows truncation in the
850* computed upper triangular factor, the matrix R is examined and,
851* if possible, replaced with its leading upper trapezoidal part.
852*
853 epsln = slamch('E')
854 sfmin = slamch('S')
855* SMALL = SFMIN / EPSLN
856 nr = n
857*
858 IF ( accla ) THEN
859*
860* Standard absolute error bound suffices. All sigma_i with
861* sigma_i < N*EPS*||A||_F are flushed to zero. This is an
862* aggressive enforcement of lower numerical rank by introducing a
863* backward error of the order of N*EPS*||A||_F.
864 nr = 1
865 rtmp = sqrt(real(n))*epsln
866 DO 3001 p = 2, n
867 IF ( abs(a(p,p)) .LT. (rtmp*abs(a(1,1))) ) GO TO 3002
868 nr = nr + 1
869 3001 CONTINUE
870 3002 CONTINUE
871*
872 ELSEIF ( acclm ) THEN
873* .. similarly as above, only slightly more gentle (less aggressive).
874* Sudden drop on the diagonal of R is used as the criterion for being
875* close-to-rank-deficient. The threshold is set to EPSLN=SLAMCH('E').
876* [[This can be made more flexible by replacing this hard-coded value
877* with a user specified threshold.]] Also, the values that underflow
878* will be truncated.
879 nr = 1
880 DO 3401 p = 2, n
881 IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
882 $ ( abs(a(p,p)) .LT. sfmin ) ) GO TO 3402
883 nr = nr + 1
884 3401 CONTINUE
885 3402 CONTINUE
886*
887 ELSE
888* .. RRQR not authorized to determine numerical rank except in the
889* obvious case of zero pivots.
890* .. inspect R for exact zeros on the diagonal;
891* R(i,i)=0 => R(i:N,i:N)=0.
892 nr = 1
893 DO 3501 p = 2, n
894 IF ( abs(a(p,p)) .EQ. zero ) GO TO 3502
895 nr = nr + 1
896 3501 CONTINUE
897 3502 CONTINUE
898*
899 IF ( conda ) THEN
900* Estimate the scaled condition number of A. Use the fact that it is
901* the same as the scaled condition number of R.
902* .. V is used as workspace
903 CALL slacpy( 'U', n, n, a, lda, v, ldv )
904* Only the leading NR x NR submatrix of the triangular factor
905* is considered. Only if NR=N will this give a reliable error
906* bound. However, even for NR < N, this can be used on an
907* expert level and obtain useful information in the sense of
908* perturbation theory.
909 DO 3053 p = 1, nr
910 rtmp = snrm2( p, v(1,p), 1 )
911 CALL sscal( p, one/rtmp, v(1,p), 1 )
912 3053 CONTINUE
913 IF ( .NOT. ( lsvec .OR. rsvec ) ) THEN
914 CALL spocon( 'U', nr, v, ldv, one, rtmp,
915 $ work, iwork(n+iwoff), ierr )
916 ELSE
917 CALL spocon( 'U', nr, v, ldv, one, rtmp,
918 $ work(n+1), iwork(n+iwoff), ierr )
919 END IF
920 sconda = one / sqrt(rtmp)
921* For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
922* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
923* See the reference [1] for more details.
924 END IF
925*
926 ENDIF
927*
928 IF ( wntur ) THEN
929 n1 = nr
930 ELSE IF ( wntus .OR. wntuf) THEN
931 n1 = n
932 ELSE IF ( wntua ) THEN
933 n1 = m
934 END IF
935*
936 IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
937*.......................................................................
938* .. only the singular values are requested
939*.......................................................................
940 IF ( rtrans ) THEN
941*
942* .. compute the singular values of R**T = [A](1:NR,1:N)**T
943* .. set the lower triangle of [A] to [A](1:NR,1:N)**T and
944* the upper triangle of [A] to zero.
945 DO 1146 p = 1, min( n, nr )
946 DO 1147 q = p + 1, n
947 a(q,p) = a(p,q)
948 IF ( q .LE. nr ) a(p,q) = zero
949 1147 CONTINUE
950 1146 CONTINUE
951*
952 CALL sgesvd( 'N', 'N', n, nr, a, lda, s, u, ldu,
953 $ v, ldv, work, lwork, info )
954*
955 ELSE
956*
957* .. compute the singular values of R = [A](1:NR,1:N)
958*
959 IF ( nr .GT. 1 )
960 $ CALL slaset( 'L', nr-1,nr-1, zero,zero, a(2,1), lda )
961 CALL sgesvd( 'N', 'N', nr, n, a, lda, s, u, ldu,
962 $ v, ldv, work, lwork, info )
963*
964 END IF
965*
966 ELSE IF ( lsvec .AND. ( .NOT. rsvec) ) THEN
967*.......................................................................
968* .. the singular values and the left singular vectors requested
969*.......................................................................""""""""
970 IF ( rtrans ) THEN
971* .. apply SGESVD to R**T
972* .. copy R**T into [U] and overwrite [U] with the right singular
973* vectors of R
974 DO 1192 p = 1, nr
975 DO 1193 q = p, n
976 u(q,p) = a(p,q)
977 1193 CONTINUE
978 1192 CONTINUE
979 IF ( nr .GT. 1 )
980 $ CALL slaset( 'U', nr-1,nr-1, zero,zero, u(1,2), ldu )
981* .. the left singular vectors not computed, the NR right singular
982* vectors overwrite [U](1:NR,1:NR) as transposed. These
983* will be pre-multiplied by Q to build the left singular vectors of A.
984 CALL sgesvd( 'N', 'O', n, nr, u, ldu, s, u, ldu,
985 $ u, ldu, work(n+1), lwork-n, info )
986*
987 DO 1119 p = 1, nr
988 DO 1120 q = p + 1, nr
989 rtmp = u(q,p)
990 u(q,p) = u(p,q)
991 u(p,q) = rtmp
992 1120 CONTINUE
993 1119 CONTINUE
994*
995 ELSE
996* .. apply SGESVD to R
997* .. copy R into [U] and overwrite [U] with the left singular vectors
998 CALL slacpy( 'U', nr, n, a, lda, u, ldu )
999 IF ( nr .GT. 1 )
1000 $ CALL slaset( 'L', nr-1, nr-1, zero, zero, u(2,1),
1001 $ ldu )
1002* .. the right singular vectors not computed, the NR left singular
1003* vectors overwrite [U](1:NR,1:NR)
1004 CALL sgesvd( 'O', 'N', nr, n, u, ldu, s, u, ldu,
1005 $ v, ldv, work(n+1), lwork-n, info )
1006* .. now [U](1:NR,1:NR) contains the NR left singular vectors of
1007* R. These will be pre-multiplied by Q to build the left singular
1008* vectors of A.
1009 END IF
1010*
1011* .. assemble the left singular vector matrix U of dimensions
1012* (M x NR) or (M x N) or (M x M).
1013 IF ( ( nr .LT. m ) .AND. ( .NOT.wntuf ) ) THEN
1014 CALL slaset('A', m-nr, nr, zero, zero, u(nr+1,1), ldu)
1015 IF ( nr .LT. n1 ) THEN
1016 CALL slaset( 'A',nr,n1-nr,zero,zero,u(1,nr+1), ldu )
1017 CALL slaset( 'A',m-nr,n1-nr,zero,one,
1018 $ u(nr+1,nr+1), ldu )
1019 END IF
1020 END IF
1021*
1022* The Q matrix from the first QRF is built into the left singular
1023* vectors matrix U.
1024*
1025 IF ( .NOT.wntuf )
1026 $ CALL sormqr( 'L', 'N', m, n1, n, a, lda, work, u,
1027 $ ldu, work(n+1), lwork-n, ierr )
1028 IF ( rowprm .AND. .NOT.wntuf )
1029 $ CALL slaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1030*
1031 ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
1032*.......................................................................
1033* .. the singular values and the right singular vectors requested
1034*.......................................................................
1035 IF ( rtrans ) THEN
1036* .. apply SGESVD to R**T
1037* .. copy R**T into V and overwrite V with the left singular vectors
1038 DO 1165 p = 1, nr
1039 DO 1166 q = p, n
1040 v(q,p) = (a(p,q))
1041 1166 CONTINUE
1042 1165 CONTINUE
1043 IF ( nr .GT. 1 )
1044 $ CALL slaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1045* .. the left singular vectors of R**T overwrite V, the right singular
1046* vectors not computed
1047 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1048 CALL sgesvd( 'O', 'N', n, nr, v, ldv, s, u, ldu,
1049 $ u, ldu, work(n+1), lwork-n, info )
1050*
1051 DO 1121 p = 1, nr
1052 DO 1122 q = p + 1, nr
1053 rtmp = v(q,p)
1054 v(q,p) = v(p,q)
1055 v(p,q) = rtmp
1056 1122 CONTINUE
1057 1121 CONTINUE
1058*
1059 IF ( nr .LT. n ) THEN
1060 DO 1103 p = 1, nr
1061 DO 1104 q = nr + 1, n
1062 v(p,q) = v(q,p)
1063 1104 CONTINUE
1064 1103 CONTINUE
1065 END IF
1066 CALL slapmt( .false., nr, n, v, ldv, iwork )
1067 ELSE
1068* .. need all N right singular vectors and NR < N
1069* [!] This is simple implementation that augments [V](1:N,1:NR)
1070* by padding a zero block. In the case NR << N, a more efficient
1071* way is to first use the QR factorization. For more details
1072* how to implement this, see the " FULL SVD " branch.
1073 CALL slaset('G', n, n-nr, zero, zero, v(1,nr+1), ldv)
1074 CALL sgesvd( 'O', 'N', n, n, v, ldv, s, u, ldu,
1075 $ u, ldu, work(n+1), lwork-n, info )
1076*
1077 DO 1123 p = 1, n
1078 DO 1124 q = p + 1, n
1079 rtmp = v(q,p)
1080 v(q,p) = v(p,q)
1081 v(p,q) = rtmp
1082 1124 CONTINUE
1083 1123 CONTINUE
1084 CALL slapmt( .false., n, n, v, ldv, iwork )
1085 END IF
1086*
1087 ELSE
1088* .. aply SGESVD to R
1089* .. copy R into V and overwrite V with the right singular vectors
1090 CALL slacpy( 'U', nr, n, a, lda, v, ldv )
1091 IF ( nr .GT. 1 )
1092 $ CALL slaset( 'L', nr-1, nr-1, zero, zero, v(2,1),
1093 $ ldv )
1094* .. the right singular vectors overwrite V, the NR left singular
1095* vectors stored in U(1:NR,1:NR)
1096 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1097 CALL sgesvd( 'N', 'O', nr, n, v, ldv, s, u, ldu,
1098 $ v, ldv, work(n+1), lwork-n, info )
1099 CALL slapmt( .false., nr, n, v, ldv, iwork )
1100* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
1101 ELSE
1102* .. need all N right singular vectors and NR < N
1103* [!] This is simple implementation that augments [V](1:NR,1:N)
1104* by padding a zero block. In the case NR << N, a more efficient
1105* way is to first use the LQ factorization. For more details
1106* how to implement this, see the " FULL SVD " branch.
1107 CALL slaset('G', n-nr, n, zero,zero, v(nr+1,1), ldv)
1108 CALL sgesvd( 'N', 'O', n, n, v, ldv, s, u, ldu,
1109 $ v, ldv, work(n+1), lwork-n, info )
1110 CALL slapmt( .false., n, n, v, ldv, iwork )
1111 END IF
1112* .. now [V] contains the transposed matrix of the right singular
1113* vectors of A.
1114 END IF
1115*
1116 ELSE
1117*.......................................................................
1118* .. FULL SVD requested
1119*.......................................................................
1120 IF ( rtrans ) THEN
1121*
1122* .. apply SGESVD to R**T [[this option is left for R&D&T]]
1123*
1124 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1125* .. copy R**T into [V] and overwrite [V] with the left singular
1126* vectors of R**T
1127 DO 1168 p = 1, nr
1128 DO 1169 q = p, n
1129 v(q,p) = a(p,q)
1130 1169 CONTINUE
1131 1168 CONTINUE
1132 IF ( nr .GT. 1 )
1133 $ CALL slaset( 'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1134*
1135* .. the left singular vectors of R**T overwrite [V], the NR right
1136* singular vectors of R**T stored in [U](1:NR,1:NR) as transposed
1137 CALL sgesvd( 'O', 'A', n, nr, v, ldv, s, v, ldv,
1138 $ u, ldu, work(n+1), lwork-n, info )
1139* .. assemble V
1140 DO 1115 p = 1, nr
1141 DO 1116 q = p + 1, nr
1142 rtmp = v(q,p)
1143 v(q,p) = v(p,q)
1144 v(p,q) = rtmp
1145 1116 CONTINUE
1146 1115 CONTINUE
1147 IF ( nr .LT. n ) THEN
1148 DO 1101 p = 1, nr
1149 DO 1102 q = nr+1, n
1150 v(p,q) = v(q,p)
1151 1102 CONTINUE
1152 1101 CONTINUE
1153 END IF
1154 CALL slapmt( .false., nr, n, v, ldv, iwork )
1155*
1156 DO 1117 p = 1, nr
1157 DO 1118 q = p + 1, nr
1158 rtmp = u(q,p)
1159 u(q,p) = u(p,q)
1160 u(p,q) = rtmp
1161 1118 CONTINUE
1162 1117 CONTINUE
1163*
1164 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1165 CALL slaset('A', m-nr,nr, zero,zero, u(nr+1,1),
1166 $ ldu)
1167 IF ( nr .LT. n1 ) THEN
1168 CALL slaset('A',nr,n1-nr,zero,zero,u(1,nr+1),
1169 $ ldu)
1170 CALL slaset( 'A',m-nr,n1-nr,zero,one,
1171 $ u(nr+1,nr+1), ldu )
1172 END IF
1173 END IF
1174*
1175 ELSE
1176* .. need all N right singular vectors and NR < N
1177* .. copy R**T into [V] and overwrite [V] with the left singular
1178* vectors of R**T
1179* [[The optimal ratio N/NR for using QRF instead of padding
1180* with zeros. Here hard coded to 2; it must be at least
1181* two due to work space constraints.]]
1182* OPTRATIO = ILAENV(6, 'SGESVD', 'S' // 'O', NR,N,0,0)
1183* OPTRATIO = MAX( OPTRATIO, 2 )
1184 optratio = 2
1185 IF ( optratio*nr .GT. n ) THEN
1186 DO 1198 p = 1, nr
1187 DO 1199 q = p, n
1188 v(q,p) = a(p,q)
1189 1199 CONTINUE
1190 1198 CONTINUE
1191 IF ( nr .GT. 1 )
1192 $ CALL slaset('U',nr-1,nr-1, zero,zero, v(1,2),ldv)
1193*
1194 CALL slaset('A',n,n-nr,zero,zero,v(1,nr+1),ldv)
1195 CALL sgesvd( 'O', 'A', n, n, v, ldv, s, v, ldv,
1196 $ u, ldu, work(n+1), lwork-n, info )
1197*
1198 DO 1113 p = 1, n
1199 DO 1114 q = p + 1, n
1200 rtmp = v(q,p)
1201 v(q,p) = v(p,q)
1202 v(p,q) = rtmp
1203 1114 CONTINUE
1204 1113 CONTINUE
1205 CALL slapmt( .false., n, n, v, ldv, iwork )
1206* .. assemble the left singular vector matrix U of dimensions
1207* (M x N1), i.e. (M x N) or (M x M).
1208*
1209 DO 1111 p = 1, n
1210 DO 1112 q = p + 1, n
1211 rtmp = u(q,p)
1212 u(q,p) = u(p,q)
1213 u(p,q) = rtmp
1214 1112 CONTINUE
1215 1111 CONTINUE
1216*
1217 IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1218 CALL slaset('A',m-n,n,zero,zero,u(n+1,1),ldu)
1219 IF ( n .LT. n1 ) THEN
1220 CALL slaset('A',n,n1-n,zero,zero,u(1,n+1),
1221 $ ldu)
1222 CALL slaset('A',m-n,n1-n,zero,one,
1223 $ u(n+1,n+1), ldu )
1224 END IF
1225 END IF
1226 ELSE
1227* .. copy R**T into [U] and overwrite [U] with the right
1228* singular vectors of R
1229 DO 1196 p = 1, nr
1230 DO 1197 q = p, n
1231 u(q,nr+p) = a(p,q)
1232 1197 CONTINUE
1233 1196 CONTINUE
1234 IF ( nr .GT. 1 )
1235 $ CALL slaset('U',nr-1,nr-1,zero,zero,u(1,nr+2),ldu)
1236 CALL sgeqrf( n, nr, u(1,nr+1), ldu, work(n+1),
1237 $ work(n+nr+1), lwork-n-nr, ierr )
1238 DO 1143 p = 1, nr
1239 DO 1144 q = 1, n
1240 v(q,p) = u(p,nr+q)
1241 1144 CONTINUE
1242 1143 CONTINUE
1243 CALL slaset('U',nr-1,nr-1,zero,zero,v(1,2),ldv)
1244 CALL sgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1245 $ v,ldv, work(n+nr+1),lwork-n-nr, info )
1246 CALL slaset('A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1247 CALL slaset('A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1248 CALL slaset('A',n-nr,n-nr,zero,one,v(nr+1,nr+1),
1249 $ ldv)
1250 CALL sormqr('R','C', n, n, nr, u(1,nr+1), ldu,
1251 $ work(n+1),v,ldv,work(n+nr+1),lwork-n-nr,ierr)
1252 CALL slapmt( .false., n, n, v, ldv, iwork )
1253* .. assemble the left singular vector matrix U of dimensions
1254* (M x NR) or (M x N) or (M x M).
1255 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1256 CALL slaset('A',m-nr,nr,zero,zero,u(nr+1,1),ldu)
1257 IF ( nr .LT. n1 ) THEN
1258 CALL slaset('A',nr,n1-nr,zero,zero,u(1,nr+1),
1259 $ ldu)
1260 CALL slaset( 'A',m-nr,n1-nr,zero,one,
1261 $ u(nr+1,nr+1),ldu)
1262 END IF
1263 END IF
1264 END IF
1265 END IF
1266*
1267 ELSE
1268*
1269* .. apply SGESVD to R [[this is the recommended option]]
1270*
1271 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1272* .. copy R into [V] and overwrite V with the right singular vectors
1273 CALL slacpy( 'U', nr, n, a, lda, v, ldv )
1274 IF ( nr .GT. 1 )
1275 $ CALL slaset( 'L', nr-1,nr-1, zero,zero, v(2,1), ldv )
1276* .. the right singular vectors of R overwrite [V], the NR left
1277* singular vectors of R stored in [U](1:NR,1:NR)
1278 CALL sgesvd( 'S', 'O', nr, n, v, ldv, s, u, ldu,
1279 $ v, ldv, work(n+1), lwork-n, info )
1280 CALL slapmt( .false., nr, n, v, ldv, iwork )
1281* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**T
1282* .. assemble the left singular vector matrix U of dimensions
1283* (M x NR) or (M x N) or (M x M).
1284 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1285 CALL slaset('A', m-nr,nr, zero,zero, u(nr+1,1),
1286 $ ldu)
1287 IF ( nr .LT. n1 ) THEN
1288 CALL slaset('A',nr,n1-nr,zero,zero,u(1,nr+1),
1289 $ ldu)
1290 CALL slaset( 'A',m-nr,n1-nr,zero,one,
1291 $ u(nr+1,nr+1), ldu )
1292 END IF
1293 END IF
1294*
1295 ELSE
1296* .. need all N right singular vectors and NR < N
1297* .. the requested number of the left singular vectors
1298* is then N1 (N or M)
1299* [[The optimal ratio N/NR for using LQ instead of padding
1300* with zeros. Here hard coded to 2; it must be at least
1301* two due to work space constraints.]]
1302* OPTRATIO = ILAENV(6, 'SGESVD', 'S' // 'O', NR,N,0,0)
1303* OPTRATIO = MAX( OPTRATIO, 2 )
1304 optratio = 2
1305 IF ( optratio * nr .GT. n ) THEN
1306 CALL slacpy( 'U', nr, n, a, lda, v, ldv )
1307 IF ( nr .GT. 1 )
1308 $ CALL slaset('L', nr-1,nr-1, zero,zero, v(2,1),ldv)
1309* .. the right singular vectors of R overwrite [V], the NR left
1310* singular vectors of R stored in [U](1:NR,1:NR)
1311 CALL slaset('A', n-nr,n, zero,zero, v(nr+1,1),ldv)
1312 CALL sgesvd( 'S', 'O', n, n, v, ldv, s, u, ldu,
1313 $ v, ldv, work(n+1), lwork-n, info )
1314 CALL slapmt( .false., n, n, v, ldv, iwork )
1315* .. now [V] contains the transposed matrix of the right
1316* singular vectors of A. The leading N left singular vectors
1317* are in [U](1:N,1:N)
1318* .. assemble the left singular vector matrix U of dimensions
1319* (M x N1), i.e. (M x N) or (M x M).
1320 IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1321 CALL slaset('A',m-n,n,zero,zero,u(n+1,1),ldu)
1322 IF ( n .LT. n1 ) THEN
1323 CALL slaset('A',n,n1-n,zero,zero,u(1,n+1),
1324 $ ldu)
1325 CALL slaset( 'A',m-n,n1-n,zero,one,
1326 $ u(n+1,n+1), ldu )
1327 END IF
1328 END IF
1329 ELSE
1330 CALL slacpy( 'U', nr, n, a, lda, u(nr+1,1), ldu )
1331 IF ( nr .GT. 1 )
1332 $ CALL slaset('L',nr-1,nr-1,zero,zero,u(nr+2,1),ldu)
1333 CALL sgelqf( nr, n, u(nr+1,1), ldu, work(n+1),
1334 $ work(n+nr+1), lwork-n-nr, ierr )
1335 CALL slacpy('L',nr,nr,u(nr+1,1),ldu,v,ldv)
1336 IF ( nr .GT. 1 )
1337 $ CALL slaset('U',nr-1,nr-1,zero,zero,v(1,2),ldv)
1338 CALL sgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1339 $ v, ldv, work(n+nr+1), lwork-n-nr, info )
1340 CALL slaset('A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1341 CALL slaset('A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1342 CALL slaset('A',n-nr,n-nr,zero,one,v(nr+1,nr+1),
1343 $ ldv)
1344 CALL sormlq('R','N',n,n,nr,u(nr+1,1),ldu,work(n+1),
1345 $ v, ldv, work(n+nr+1),lwork-n-nr,ierr)
1346 CALL slapmt( .false., n, n, v, ldv, iwork )
1347* .. assemble the left singular vector matrix U of dimensions
1348* (M x NR) or (M x N) or (M x M).
1349 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1350 CALL slaset('A',m-nr,nr,zero,zero,u(nr+1,1),ldu)
1351 IF ( nr .LT. n1 ) THEN
1352 CALL slaset('A',nr,n1-nr,zero,zero,u(1,nr+1),
1353 $ ldu)
1354 CALL slaset( 'A',m-nr,n1-nr,zero,one,
1355 $ u(nr+1,nr+1), ldu )
1356 END IF
1357 END IF
1358 END IF
1359 END IF
1360* .. end of the "R**T or R" branch
1361 END IF
1362*
1363* The Q matrix from the first QRF is built into the left singular
1364* vectors matrix U.
1365*
1366 IF ( .NOT. wntuf )
1367 $ CALL sormqr( 'L', 'N', m, n1, n, a, lda, work, u,
1368 $ ldu, work(n+1), lwork-n, ierr )
1369 IF ( rowprm .AND. .NOT.wntuf )
1370 $ CALL slaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1371*
1372* ... end of the "full SVD" branch
1373 END IF
1374*
1375* Check whether some singular values are returned as zeros, e.g.
1376* due to underflow, and update the numerical rank.
1377 p = nr
1378 DO 4001 q = p, 1, -1
1379 IF ( s(q) .GT. zero ) GO TO 4002
1380 nr = nr - 1
1381 4001 CONTINUE
1382 4002 CONTINUE
1383*
1384* .. if numerical rank deficiency is detected, the truncated
1385* singular values are set to zero.
1386 IF ( nr .LT. n ) CALL slaset( 'G', n-nr,1, zero,zero, s(nr+1),
1387 $ n )
1388* .. undo scaling; this may cause overflow in the largest singular
1389* values.
1390 IF ( ascaled )
1391 $ CALL slascl( 'G',0,0, one,sqrt(real(m)), nr,1, s, n, ierr )
1392 IF ( conda ) rwork(1) = sconda
1393 rwork(2) = real( p - nr )
1394* .. p-NR is the number of singular values that are computed as
1395* exact zeros in SGESVD() applied to the (possibly truncated)
1396* full row rank triangular (trapezoidal) factor of A.
1397 numrank = nr
1398*
1399 RETURN
1400*
1401* End of SGESVDQ
1402*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgelqf
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:142
sgeqp3
subroutine sgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
SGEQP3
Definition sgeqp3.f:149
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:144
sgesvd
subroutine sgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
SGESVD computes the singular value decomposition (SVD) for GE matrices
Definition sgesvd.f:210
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:112
slapmt
subroutine slapmt(forwrd, m, n, x, ldx, k)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition slapmt.f:102
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
slaswp
subroutine slaswp(n, a, lda, k1, k2, ipiv, incx)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition slaswp.f:113
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
spocon
subroutine spocon(uplo, n, a, lda, anorm, rcond, work, iwork, info)
SPOCON
Definition spocon.f:119
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sormlq
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
Definition sormlq.f:166
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:166
Here is the call graph for this function:
Here is the caller graph for this function: