LAPACK 3.11.0
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, enlosed 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 412 of file sgesvdq.f.

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,
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*
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:143
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:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:114
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
Definition: sgeqp3.f:151
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:146
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
Definition: sgelqf.f:143
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:211
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:104
subroutine slaswp(N, A, LDA, K1, K2, IPIV, INCX)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: slaswp.f:115
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
Definition: sormlq.f:168
subroutine spocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SPOCON
Definition: spocon.f:121
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: