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

◆ dgesvdq()

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

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

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

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