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

◆ cgesvdq()

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

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

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

Purpose:
!>
!> CGESVDQ computes the singular value decomposition (SVD) of a complex
!> 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 unitary 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, CGESVD is applied to
!>          the adjoint R**H 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 CGESVD. 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**H , 0)**H. 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 COMPLEX 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 CGEQP3. If JOBU = 'F', these Householder
!>          vectors together with CWORK(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 COMPLEX 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 unitary 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 COMPLEX 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 unitary matrix  V**H;
!>          If JOBV = 'R', V contains the first NUMRANK rows of V**H (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 CGESVD. 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, LCWORK, 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';
!>          LIWORK >= N           if JOBP = 'N'.
!>
!>          If LIWORK = -1, then a workspace query is assumed; the routine
!>          only calculates and returns the optimal and minimal sizes
!>          for the CWORK, IWORK, and RWORK arrays, and no error
!>          message related to LCWORK is issued by XERBLA.
!>
[out]CWORK
!>          CWORK is COMPLEX array, dimension (max(2, LCWORK)), used as a workspace.
!>          On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters
!>          needed to recover the Q factor from the QR factorization computed by
!>          CGEQP3.
!>
!>          If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
!>          CWORK(1) returns the optimal LCWORK, and
!>          CWORK(2) returns the minimal LCWORK.
!>
[in,out]LCWORK
!>          LCWORK is INTEGER
!>          The dimension of the array CWORK. It is determined as follows:
!>          Let  LWQP3 = N+1,  LWCON = 2*N, and let
!>          LWUNQ = { MAX( N, 1 ),  if JOBU = 'R', 'S', or 'U'
!>                  { MAX( M, 1 ),  if JOBU = 'A'
!>          LWSVD = MAX( 3*N, 1 )
!>          LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ),
!>          LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )
!>          Then the minimal value of LCWORK 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, LWUNQ ) if the singular values and the left
!>                                   singular vectors are requested;
!>          = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) 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, LWUNQ ) if the full SVD is requested with JOBV = 'R';
!>                                   independent of JOBR;
!>          = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is requested,
!>                                   JOBV = 'R' and, also a scaled condition
!>                                   estimate requested; independent of JOBR;
!>          = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
!>         N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) ) if the
!>                         full SVD is requested with JOBV = 'A' or 'V', and
!>                         JOBR ='N'
!>          = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
!>         N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ ) )
!>                         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, LWUNQ ),
!>         N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) if the
!>                         full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
!>          = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
!>         N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) )
!>                         if the full SVD is requested with JOBV = 'A', 'V' and
!>                         JOBR ='T', and also a scaled condition number estimate
!>                         requested.
!>          Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ).
!>
!>          If LCWORK = -1, then a workspace query is assumed; the routine
!>          only calculates and returns the optimal and minimal sizes
!>          for the CWORK, IWORK, and RWORK arrays, and no error
!>          message related to LCWORK 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 CGESVD applied to the upper triangular or trapezoidal
!>          R (from the initial QR factorization). In case of early exit (no call to
!>          CGESVD, such as in the case of zero matrix) RWORK(2) = -1.
!>
!>          If LIWORK, LCWORK, 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, 5*N);
!>          Otherwise, LRWORK >= MAX(2, 5*N).
!>
!>          If LRWORK = -1, then a workspace query is assumed; the routine
!>          only calculates and returns the optimal and minimal sizes
!>          for the CWORK, IWORK, and RWORK arrays, and no error
!>          message related to LCWORK 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 CBDSQR did not converge, INFO specifies how many superdiagonals
!>          of an intermediate bidiagonal form B (computed in CGESVD) 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 408 of file cgesvdq.f.

411* .. Scalar Arguments ..
412 IMPLICIT NONE
413 CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
414 INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK,
415 $ INFO
416* ..
417* .. Array Arguments ..
418 COMPLEX A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
419 REAL S( * ), RWORK( * )
420 INTEGER IWORK( * )
421*
422* =====================================================================
423*
424* .. Parameters ..
425 REAL ZERO, ONE
426 parameter( zero = 0.0e0, one = 1.0e0 )
427 COMPLEX CZERO, CONE
428 parameter( czero = ( 0.0e0, 0.0e0 ), cone = ( 1.0e0, 0.0e0 ) )
429* ..
430* .. Local Scalars ..
431 INTEGER IERR, NR, N1, OPTRATIO, p, q
432 INTEGER LWCON, LWQP3, LWRK_CGELQF, LWRK_CGESVD, LWRK_CGESVD2,
433 $ LWRK_CGEQP3, LWRK_CGEQRF, LWRK_CUNMLQ, LWRK_CUNMQR,
434 $ LWRK_CUNMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWUNQ,
435 $ LWUNQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2,
436 $ IMINWRK, RMINWRK
437 LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
438 $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
439 $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR
440 REAL BIG, EPSLN, RTMP, SCONDA, SFMIN
441 COMPLEX CTMP
442* ..
443* .. Local Arrays
444 COMPLEX CDUMMY(1)
445 REAL RDUMMY(1)
446* ..
447* .. External Subroutines (BLAS, LAPACK)
448 EXTERNAL cgelqf, cgeqp3, cgeqrf, cgesvd, clacpy,
449 $ clapmt, clascl, claset, claswp, csscal,
450 $ slaset, slascl, cpocon, cunmlq, cunmqr,
451 $ xerbla
452* ..
453* .. External Functions (BLAS, LAPACK)
454 LOGICAL LSAME
455 INTEGER ISAMAX
456 REAL CLANGE, SCNRM2, SLAMCH
457 EXTERNAL clange, lsame, isamax, scnrm2, slamch
458* ..
459* .. Intrinsic Functions ..
460 INTRINSIC abs, conjg, max, min, real, sqrt
461* ..
462* .. Executable Statements ..
463*
464* Test the input arguments
465*
466 wntus = lsame( jobu, 'S' ) .OR. lsame( jobu, 'U' )
467 wntur = lsame( jobu, 'R' )
468 wntua = lsame( jobu, 'A' )
469 wntuf = lsame( jobu, 'F' )
470 lsvc0 = wntus .OR. wntur .OR. wntua
471 lsvec = lsvc0 .OR. wntuf
472 dntwu = lsame( jobu, 'N' )
473*
474 wntvr = lsame( jobv, 'R' )
475 wntva = lsame( jobv, 'A' ) .OR. lsame( jobv, 'V' )
476 rsvec = wntvr .OR. wntva
477 dntwv = lsame( jobv, 'N' )
478*
479 accla = lsame( joba, 'A' )
480 acclm = lsame( joba, 'M' )
481 conda = lsame( joba, 'E' )
482 acclh = lsame( joba, 'H' ) .OR. conda
483*
484 rowprm = lsame( jobp, 'P' )
485 rtrans = lsame( jobr, 'T' )
486*
487 IF ( rowprm ) THEN
488 iminwrk = max( 1, n + m - 1 )
489 rminwrk = max( 2, m, 5*n )
490 ELSE
491 iminwrk = max( 1, n )
492 rminwrk = max( 2, 5*n )
493 END IF
494 lquery = (liwork .EQ. -1 .OR. lcwork .EQ. -1 .OR. lrwork .EQ. -1)
495 info = 0
496 IF ( .NOT. ( accla .OR. acclm .OR. acclh ) ) THEN
497 info = -1
498 ELSE IF ( .NOT.( rowprm .OR. lsame( jobp, 'N' ) ) ) THEN
499 info = -2
500 ELSE IF ( .NOT.( rtrans .OR. lsame( jobr, 'N' ) ) ) THEN
501 info = -3
502 ELSE IF ( .NOT.( lsvec .OR. dntwu ) ) THEN
503 info = -4
504 ELSE IF ( wntur .AND. wntva ) THEN
505 info = -5
506 ELSE IF ( .NOT.( rsvec .OR. dntwv )) THEN
507 info = -5
508 ELSE IF ( m.LT.0 ) THEN
509 info = -6
510 ELSE IF ( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
511 info = -7
512 ELSE IF ( lda.LT.max( 1, m ) ) THEN
513 info = -9
514 ELSE IF ( ldu.LT.1 .OR. ( lsvc0 .AND. ldu.LT.m ) .OR.
515 $ ( wntuf .AND. ldu.LT.n ) ) THEN
516 info = -12
517 ELSE IF ( ldv.LT.1 .OR. ( rsvec .AND. ldv.LT.n ) .OR.
518 $ ( conda .AND. ldv.LT.n ) ) THEN
519 info = -14
520 ELSE IF ( liwork .LT. iminwrk .AND. .NOT. lquery ) THEN
521 info = -17
522 END IF
523*
524*
525 IF ( info .EQ. 0 ) THEN
526*
527* Compute workspace
528* .. compute the minimal and the optimal workspace lengths
529* [[The expressions for computing the minimal and the optimal
530* values of LCWORK are written with a lot of redundancy and
531* can be simplified. However, this detailed form is easier for
532* maintenance and modifications of the code.]]
533*
534* .. minimal workspace length for CGEQP3 of an M x N matrix
535 lwqp3 = n+1
536* .. minimal workspace length for CUNMQR to build left singular vectors
537 IF ( wntus .OR. wntur ) THEN
538 lwunq = max( n , 1 )
539 ELSE IF ( wntua ) THEN
540 lwunq = max( m , 1 )
541 END IF
542* .. minimal workspace length for CPOCON of an N x N matrix
543 lwcon = 2 * n
544* .. CGESVD of an N x N matrix
545 lwsvd = max( 3 * n, 1 )
546 IF ( lquery ) THEN
547 CALL cgeqp3( m, n, a, lda, iwork, cdummy, cdummy, -1,
548 $ rdummy, ierr )
549 lwrk_cgeqp3 = int( cdummy(1) )
550 IF ( wntus .OR. wntur ) THEN
551 CALL cunmqr( 'L', 'N', m, n, n, a, lda, cdummy, u,
552 $ ldu, cdummy, -1, ierr )
553 lwrk_cunmqr = int( cdummy(1) )
554 ELSE IF ( wntua ) THEN
555 CALL cunmqr( 'L', 'N', m, m, n, a, lda, cdummy, u,
556 $ ldu, cdummy, -1, ierr )
557 lwrk_cunmqr = int( cdummy(1) )
558 ELSE
559 lwrk_cunmqr = 0
560 END IF
561 END IF
562 minwrk = 2
563 optwrk = 2
564 IF ( .NOT. (lsvec .OR. rsvec )) THEN
565* .. minimal and optimal sizes of the complex workspace if
566* only the singular values are requested
567 IF ( conda ) THEN
568 minwrk = max( n+lwqp3, lwcon, lwsvd )
569 ELSE
570 minwrk = max( n+lwqp3, lwsvd )
571 END IF
572 IF ( lquery ) THEN
573 CALL cgesvd( 'N', 'N', n, n, a, lda, s, u, ldu,
574 $ v, ldv, cdummy, -1, rdummy, ierr )
575 lwrk_cgesvd = int( cdummy(1) )
576 IF ( conda ) THEN
577 optwrk = max( n+lwrk_cgeqp3, n+lwcon, lwrk_cgesvd )
578 ELSE
579 optwrk = max( n+lwrk_cgeqp3, lwrk_cgesvd )
580 END IF
581 END IF
582 ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
583* .. minimal and optimal sizes of the complex workspace if the
584* singular values and the left singular vectors are requested
585 IF ( conda ) THEN
586 minwrk = n + max( lwqp3, lwcon, lwsvd, lwunq )
587 ELSE
588 minwrk = n + max( lwqp3, lwsvd, lwunq )
589 END IF
590 IF ( lquery ) THEN
591 IF ( rtrans ) THEN
592 CALL cgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
593 $ v, ldv, cdummy, -1, rdummy, ierr )
594 ELSE
595 CALL cgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
596 $ v, ldv, cdummy, -1, rdummy, ierr )
597 END IF
598 lwrk_cgesvd = int( cdummy(1) )
599 IF ( conda ) THEN
600 optwrk = n + max( lwrk_cgeqp3, lwcon, lwrk_cgesvd,
601 $ lwrk_cunmqr )
602 ELSE
603 optwrk = n + max( lwrk_cgeqp3, lwrk_cgesvd,
604 $ lwrk_cunmqr )
605 END IF
606 END IF
607 ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
608* .. minimal and optimal sizes of the complex workspace if the
609* singular values and the right singular vectors are requested
610 IF ( conda ) THEN
611 minwrk = n + max( lwqp3, lwcon, lwsvd )
612 ELSE
613 minwrk = n + max( lwqp3, lwsvd )
614 END IF
615 IF ( lquery ) THEN
616 IF ( rtrans ) THEN
617 CALL cgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
618 $ v, ldv, cdummy, -1, rdummy, ierr )
619 ELSE
620 CALL cgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
621 $ v, ldv, cdummy, -1, rdummy, ierr )
622 END IF
623 lwrk_cgesvd = int( cdummy(1) )
624 IF ( conda ) THEN
625 optwrk = n + max( lwrk_cgeqp3, lwcon, lwrk_cgesvd )
626 ELSE
627 optwrk = n + max( lwrk_cgeqp3, lwrk_cgesvd )
628 END IF
629 END IF
630 ELSE
631* .. minimal and optimal sizes of the complex workspace if the
632* full SVD is requested
633 IF ( rtrans ) THEN
634 minwrk = max( lwqp3, lwsvd, lwunq )
635 IF ( conda ) minwrk = max( minwrk, lwcon )
636 minwrk = minwrk + n
637 IF ( wntva ) THEN
638* .. minimal workspace length for N x N/2 CGEQRF
639 lwqrf = max( n/2, 1 )
640* .. minimal workspace length for N/2 x N/2 CGESVD
641 lwsvd2 = max( 3 * (n/2), 1 )
642 lwunq2 = max( n, 1 )
643 minwrk2 = max( lwqp3, n/2+lwqrf, n/2+lwsvd2,
644 $ n/2+lwunq2, lwunq )
645 IF ( conda ) minwrk2 = max( minwrk2, lwcon )
646 minwrk2 = n + minwrk2
647 minwrk = max( minwrk, minwrk2 )
648 END IF
649 ELSE
650 minwrk = max( lwqp3, lwsvd, lwunq )
651 IF ( conda ) minwrk = max( minwrk, lwcon )
652 minwrk = minwrk + n
653 IF ( wntva ) THEN
654* .. minimal workspace length for N/2 x N CGELQF
655 lwlqf = max( n/2, 1 )
656 lwsvd2 = max( 3 * (n/2), 1 )
657 lwunlq = max( n , 1 )
658 minwrk2 = max( lwqp3, n/2+lwlqf, n/2+lwsvd2,
659 $ n/2+lwunlq, lwunq )
660 IF ( conda ) minwrk2 = max( minwrk2, lwcon )
661 minwrk2 = n + minwrk2
662 minwrk = max( minwrk, minwrk2 )
663 END IF
664 END IF
665 IF ( lquery ) THEN
666 IF ( rtrans ) THEN
667 CALL cgesvd( 'O', 'A', n, n, a, lda, s, u, ldu,
668 $ v, ldv, cdummy, -1, rdummy, ierr )
669 lwrk_cgesvd = int( cdummy(1) )
670 optwrk = max(lwrk_cgeqp3,lwrk_cgesvd,lwrk_cunmqr)
671 IF ( conda ) optwrk = max( optwrk, lwcon )
672 optwrk = n + optwrk
673 IF ( wntva ) THEN
674 CALL cgeqrf(n,n/2,u,ldu,cdummy,cdummy,-1,ierr)
675 lwrk_cgeqrf = int( cdummy(1) )
676 CALL cgesvd( 'S', 'O', n/2,n/2, v,ldv, s, u,
677 $ ldu,
678 $ v, ldv, cdummy, -1, rdummy, ierr )
679 lwrk_cgesvd2 = int( cdummy(1) )
680 CALL cunmqr( 'R', 'C', n, n, n/2, u, ldu,
681 $ cdummy,
682 $ v, ldv, cdummy, -1, ierr )
683 lwrk_cunmqr2 = int( cdummy(1) )
684 optwrk2 = max( lwrk_cgeqp3, n/2+lwrk_cgeqrf,
685 $ n/2+lwrk_cgesvd2, n/2+lwrk_cunmqr2 )
686 IF ( conda ) optwrk2 = max( optwrk2, lwcon )
687 optwrk2 = n + optwrk2
688 optwrk = max( optwrk, optwrk2 )
689 END IF
690 ELSE
691 CALL cgesvd( 'S', 'O', n, n, a, lda, s, u, ldu,
692 $ v, ldv, cdummy, -1, rdummy, ierr )
693 lwrk_cgesvd = int( cdummy(1) )
694 optwrk = max(lwrk_cgeqp3,lwrk_cgesvd,lwrk_cunmqr)
695 IF ( conda ) optwrk = max( optwrk, lwcon )
696 optwrk = n + optwrk
697 IF ( wntva ) THEN
698 CALL cgelqf(n/2,n,u,ldu,cdummy,cdummy,-1,ierr)
699 lwrk_cgelqf = int( cdummy(1) )
700 CALL cgesvd( 'S','O', n/2,n/2, v, ldv, s, u,
701 $ ldu,
702 $ v, ldv, cdummy, -1, rdummy, ierr )
703 lwrk_cgesvd2 = int( cdummy(1) )
704 CALL cunmlq( 'R', 'N', n, n, n/2, u, ldu,
705 $ cdummy,
706 $ v, ldv, cdummy,-1,ierr )
707 lwrk_cunmlq = int( cdummy(1) )
708 optwrk2 = max( lwrk_cgeqp3, n/2+lwrk_cgelqf,
709 $ n/2+lwrk_cgesvd2, n/2+lwrk_cunmlq )
710 IF ( conda ) optwrk2 = max( optwrk2, lwcon )
711 optwrk2 = n + optwrk2
712 optwrk = max( optwrk, optwrk2 )
713 END IF
714 END IF
715 END IF
716 END IF
717*
718 minwrk = max( 2, minwrk )
719 optwrk = max( 2, optwrk )
720 IF ( lcwork .LT. minwrk .AND. (.NOT.lquery) ) info = -19
721*
722 END IF
723*
724 IF (info .EQ. 0 .AND. lrwork .LT. rminwrk .AND. .NOT. lquery) THEN
725 info = -21
726 END IF
727 IF( info.NE.0 ) THEN
728 CALL xerbla( 'CGESVDQ', -info )
729 RETURN
730 ELSE IF ( lquery ) THEN
731*
732* Return optimal workspace
733*
734 iwork(1) = iminwrk
735 cwork(1) = cmplx( optwrk )
736 cwork(2) = cmplx( minwrk )
737 rwork(1) = real( rminwrk )
738 RETURN
739 END IF
740*
741* Quick return if the matrix is void.
742*
743 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) ) THEN
744* .. all output is void.
745 RETURN
746 END IF
747*
748 big = slamch('O')
749 ascaled = .false.
750 IF ( rowprm ) THEN
751* .. reordering the rows in decreasing sequence in the
752* ell-infinity norm - this enhances numerical robustness in
753* the case of differently scaled rows.
754 DO 1904 p = 1, m
755* RWORK(p) = ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
756* [[CLANGE will return NaN if an entry of the p-th row is Nan]]
757 rwork(p) = clange( 'M', 1, n, a(p,1), lda, rdummy )
758* .. check for NaN's and Inf's
759 IF ( ( rwork(p) .NE. rwork(p) ) .OR.
760 $ ( (rwork(p)*zero) .NE. zero ) ) THEN
761 info = - 8
762 CALL xerbla( 'CGESVDQ', -info )
763 RETURN
764 END IF
765 1904 CONTINUE
766 DO 1952 p = 1, m - 1
767 q = isamax( m-p+1, rwork(p), 1 ) + p - 1
768 iwork(n+p) = q
769 IF ( p .NE. q ) THEN
770 rtmp = rwork(p)
771 rwork(p) = rwork(q)
772 rwork(q) = rtmp
773 END IF
774 1952 CONTINUE
775*
776 IF ( rwork(1) .EQ. zero ) THEN
777* Quick return: A is the M x N zero matrix.
778 numrank = 0
779 CALL slaset( 'G', n, 1, zero, zero, s, n )
780 IF ( wntus ) CALL claset('G', m, n, czero, cone, u,
781 $ ldu)
782 IF ( wntua ) CALL claset('G', m, m, czero, cone, u,
783 $ ldu)
784 IF ( wntva ) CALL claset('G', n, n, czero, cone, v,
785 $ ldv)
786 IF ( wntuf ) THEN
787 CALL claset( 'G', n, 1, czero, czero, cwork, n )
788 CALL claset( 'G', m, n, czero, cone, 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 clascl('G',0,0,sqrt(real(m)),one, m,n, a,lda,
807 $ ierr)
808 ascaled = .true.
809 END IF
810 CALL claswp( 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 = clange( 'M', m, n, a, lda, rwork )
820 IF ( ( rtmp .NE. rtmp ) .OR.
821 $ ( (rtmp*zero) .NE. zero ) ) THEN
822 info = - 8
823 CALL xerbla( 'CGESVDQ', -info )
824 RETURN
825 END IF
826 IF ( rtmp .GT. big / sqrt(real(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 clascl('G',0,0, sqrt(real(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 cgeqp3( m, n, a, lda, iwork, cwork, cwork(n+1), lcwork-n,
845 $ rwork, 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 = slamch('E')
852 sfmin = slamch('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(real(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=SLAMCH('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 clacpy( '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 = scnrm2( p, v(1,p), 1 )
909 CALL csscal( p, one/rtmp, v(1,p), 1 )
910 3053 CONTINUE
911 IF ( .NOT. ( lsvec .OR. rsvec ) ) THEN
912 CALL cpocon( 'U', nr, v, ldv, one, rtmp,
913 $ cwork, rwork, ierr )
914 ELSE
915 CALL cpocon( 'U', nr, v, ldv, one, rtmp,
916 $ cwork(n+1), rwork, 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**H = [A](1:NR,1:N)**H
941* .. set the lower triangle of [A] to [A](1:NR,1:N)**H and
942* the upper triangle of [A] to zero.
943 DO 1146 p = 1, min( n, nr )
944 a(p,p) = conjg(a(p,p))
945 DO 1147 q = p + 1, n
946 a(q,p) = conjg(a(p,q))
947 IF ( q .LE. nr ) a(p,q) = czero
948 1147 CONTINUE
949 1146 CONTINUE
950*
951 CALL cgesvd( 'N', 'N', n, nr, a, lda, s, u, ldu,
952 $ v, ldv, cwork, lcwork, rwork, info )
953*
954 ELSE
955*
956* .. compute the singular values of R = [A](1:NR,1:N)
957*
958 IF ( nr .GT. 1 )
959 $ CALL claset( 'L', nr-1,nr-1, czero,czero, a(2,1),
960 $ lda )
961 CALL cgesvd( 'N', 'N', nr, n, a, lda, s, u, ldu,
962 $ v, ldv, cwork, lcwork, rwork, info )
963*
964 END IF
965*
966 ELSE IF ( lsvec .AND. ( .NOT. rsvec) ) THEN
967*.......................................................................
968* .. the singular values and the left singular vectors requested
969*.......................................................................""""""""
970 IF ( rtrans ) THEN
971* .. apply CGESVD to R**H
972* .. copy R**H into [U] and overwrite [U] with the right singular
973* vectors of R
974 DO 1192 p = 1, nr
975 DO 1193 q = p, n
976 u(q,p) = conjg(a(p,q))
977 1193 CONTINUE
978 1192 CONTINUE
979 IF ( nr .GT. 1 )
980 $ CALL claset( 'U', nr-1,nr-1, czero,czero, u(1,2),
981 $ ldu )
982* .. the left singular vectors not computed, the NR right singular
983* vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These
984* will be pre-multiplied by Q to build the left singular vectors of A.
985 CALL cgesvd( 'N', 'O', n, nr, u, ldu, s, u, ldu,
986 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
987*
988 DO 1119 p = 1, nr
989 u(p,p) = conjg(u(p,p))
990 DO 1120 q = p + 1, nr
991 ctmp = conjg(u(q,p))
992 u(q,p) = conjg(u(p,q))
993 u(p,q) = ctmp
994 1120 CONTINUE
995 1119 CONTINUE
996*
997 ELSE
998* .. apply CGESVD to R
999* .. copy R into [U] and overwrite [U] with the left singular vectors
1000 CALL clacpy( 'U', nr, n, a, lda, u, ldu )
1001 IF ( nr .GT. 1 )
1002 $ CALL claset( 'L', nr-1, nr-1, czero, czero, u(2,1),
1003 $ ldu )
1004* .. the right singular vectors not computed, the NR left singular
1005* vectors overwrite [U](1:NR,1:NR)
1006 CALL cgesvd( 'O', 'N', nr, n, u, ldu, s, u, ldu,
1007 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1008* .. now [U](1:NR,1:NR) contains the NR left singular vectors of
1009* R. These will be pre-multiplied by Q to build the left singular
1010* vectors of A.
1011 END IF
1012*
1013* .. assemble the left singular vector matrix U of dimensions
1014* (M x NR) or (M x N) or (M x M).
1015 IF ( ( nr .LT. m ) .AND. ( .NOT.wntuf ) ) THEN
1016 CALL claset('A', m-nr, nr, czero, czero, u(nr+1,1), ldu)
1017 IF ( nr .LT. n1 ) THEN
1018 CALL claset( 'A',nr,n1-nr,czero,czero,u(1,nr+1),
1019 $ ldu )
1020 CALL claset( 'A',m-nr,n1-nr,czero,cone,
1021 $ u(nr+1,nr+1), ldu )
1022 END IF
1023 END IF
1024*
1025* The Q matrix from the first QRF is built into the left singular
1026* vectors matrix U.
1027*
1028 IF ( .NOT.wntuf )
1029 $ CALL cunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1030 $ ldu, cwork(n+1), lcwork-n, ierr )
1031 IF ( rowprm .AND. .NOT.wntuf )
1032 $ CALL claswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1033*
1034 ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
1035*.......................................................................
1036* .. the singular values and the right singular vectors requested
1037*.......................................................................
1038 IF ( rtrans ) THEN
1039* .. apply CGESVD to R**H
1040* .. copy R**H into V and overwrite V with the left singular vectors
1041 DO 1165 p = 1, nr
1042 DO 1166 q = p, n
1043 v(q,p) = conjg(a(p,q))
1044 1166 CONTINUE
1045 1165 CONTINUE
1046 IF ( nr .GT. 1 )
1047 $ CALL claset( 'U', nr-1,nr-1, czero,czero, v(1,2),
1048 $ ldv )
1049* .. the left singular vectors of R**H overwrite V, the right singular
1050* vectors not computed
1051 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1052 CALL cgesvd( 'O', 'N', n, nr, v, ldv, s, u, ldu,
1053 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1054*
1055 DO 1121 p = 1, nr
1056 v(p,p) = conjg(v(p,p))
1057 DO 1122 q = p + 1, nr
1058 ctmp = conjg(v(q,p))
1059 v(q,p) = conjg(v(p,q))
1060 v(p,q) = ctmp
1061 1122 CONTINUE
1062 1121 CONTINUE
1063*
1064 IF ( nr .LT. n ) THEN
1065 DO 1103 p = 1, nr
1066 DO 1104 q = nr + 1, n
1067 v(p,q) = conjg(v(q,p))
1068 1104 CONTINUE
1069 1103 CONTINUE
1070 END IF
1071 CALL clapmt( .false., nr, n, v, ldv, iwork )
1072 ELSE
1073* .. need all N right singular vectors and NR < N
1074* [!] This is simple implementation that augments [V](1:N,1:NR)
1075* by padding a zero block. In the case NR << N, a more efficient
1076* way is to first use the QR factorization. For more details
1077* how to implement this, see the " FULL SVD " branch.
1078 CALL claset('G', n, n-nr, czero, czero, v(1,nr+1),
1079 $ ldv)
1080 CALL cgesvd( 'O', 'N', n, n, v, ldv, s, u, ldu,
1081 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1082*
1083 DO 1123 p = 1, n
1084 v(p,p) = conjg(v(p,p))
1085 DO 1124 q = p + 1, n
1086 ctmp = conjg(v(q,p))
1087 v(q,p) = conjg(v(p,q))
1088 v(p,q) = ctmp
1089 1124 CONTINUE
1090 1123 CONTINUE
1091 CALL clapmt( .false., n, n, v, ldv, iwork )
1092 END IF
1093*
1094 ELSE
1095* .. aply CGESVD to R
1096* .. copy R into V and overwrite V with the right singular vectors
1097 CALL clacpy( 'U', nr, n, a, lda, v, ldv )
1098 IF ( nr .GT. 1 )
1099 $ CALL claset( 'L', nr-1, nr-1, czero, czero, v(2,1),
1100 $ ldv )
1101* .. the right singular vectors overwrite V, the NR left singular
1102* vectors stored in U(1:NR,1:NR)
1103 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1104 CALL cgesvd( 'N', 'O', nr, n, v, ldv, s, u, ldu,
1105 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1106 CALL clapmt( .false., nr, n, v, ldv, iwork )
1107* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1108 ELSE
1109* .. need all N right singular vectors and NR < N
1110* [!] This is simple implementation that augments [V](1:NR,1:N)
1111* by padding a zero block. In the case NR << N, a more efficient
1112* way is to first use the LQ factorization. For more details
1113* how to implement this, see the " FULL SVD " branch.
1114 CALL claset('G', n-nr, n, czero,czero, v(nr+1,1),
1115 $ ldv)
1116 CALL cgesvd( 'N', 'O', n, n, v, ldv, s, u, ldu,
1117 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1118 CALL clapmt( .false., n, n, v, ldv, iwork )
1119 END IF
1120* .. now [V] contains the adjoint of the matrix of the right singular
1121* vectors of A.
1122 END IF
1123*
1124 ELSE
1125*.......................................................................
1126* .. FULL SVD requested
1127*.......................................................................
1128 IF ( rtrans ) THEN
1129*
1130* .. apply CGESVD to R**H [[this option is left for R&D&T]]
1131*
1132 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1133* .. copy R**H into [V] and overwrite [V] with the left singular
1134* vectors of R**H
1135 DO 1168 p = 1, nr
1136 DO 1169 q = p, n
1137 v(q,p) = conjg(a(p,q))
1138 1169 CONTINUE
1139 1168 CONTINUE
1140 IF ( nr .GT. 1 )
1141 $ CALL claset( 'U', nr-1,nr-1, czero,czero, v(1,2),
1142 $ ldv )
1143*
1144* .. the left singular vectors of R**H overwrite [V], the NR right
1145* singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate
1146* transposed
1147 CALL cgesvd( 'O', 'A', n, nr, v, ldv, s, v, ldv,
1148 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1149* .. assemble V
1150 DO 1115 p = 1, nr
1151 v(p,p) = conjg(v(p,p))
1152 DO 1116 q = p + 1, nr
1153 ctmp = conjg(v(q,p))
1154 v(q,p) = conjg(v(p,q))
1155 v(p,q) = ctmp
1156 1116 CONTINUE
1157 1115 CONTINUE
1158 IF ( nr .LT. n ) THEN
1159 DO 1101 p = 1, nr
1160 DO 1102 q = nr+1, n
1161 v(p,q) = conjg(v(q,p))
1162 1102 CONTINUE
1163 1101 CONTINUE
1164 END IF
1165 CALL clapmt( .false., nr, n, v, ldv, iwork )
1166*
1167 DO 1117 p = 1, nr
1168 u(p,p) = conjg(u(p,p))
1169 DO 1118 q = p + 1, nr
1170 ctmp = conjg(u(q,p))
1171 u(q,p) = conjg(u(p,q))
1172 u(p,q) = ctmp
1173 1118 CONTINUE
1174 1117 CONTINUE
1175*
1176 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1177 CALL claset('A', m-nr,nr, czero,czero, u(nr+1,1),
1178 $ ldu)
1179 IF ( nr .LT. n1 ) THEN
1180 CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),
1181 $ ldu)
1182 CALL claset( 'A',m-nr,n1-nr,czero,cone,
1183 $ u(nr+1,nr+1), ldu )
1184 END IF
1185 END IF
1186*
1187 ELSE
1188* .. need all N right singular vectors and NR < N
1189* .. copy R**H into [V] and overwrite [V] with the left singular
1190* vectors of R**H
1191* [[The optimal ratio N/NR for using QRF instead of padding
1192* with zeros. Here hard coded to 2; it must be at least
1193* two due to work space constraints.]]
1194* OPTRATIO = ILAENV(6, 'CGESVD', 'S' // 'O', NR,N,0,0)
1195* OPTRATIO = MAX( OPTRATIO, 2 )
1196 optratio = 2
1197 IF ( optratio*nr .GT. n ) THEN
1198 DO 1198 p = 1, nr
1199 DO 1199 q = p, n
1200 v(q,p) = conjg(a(p,q))
1201 1199 CONTINUE
1202 1198 CONTINUE
1203 IF ( nr .GT. 1 )
1204 $ CALL claset('U',nr-1,nr-1, czero,czero, v(1,2),
1205 $ ldv)
1206*
1207 CALL claset('A',n,n-nr,czero,czero,v(1,nr+1),ldv)
1208 CALL cgesvd( 'O', 'A', n, n, v, ldv, s, v, ldv,
1209 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1210*
1211 DO 1113 p = 1, n
1212 v(p,p) = conjg(v(p,p))
1213 DO 1114 q = p + 1, n
1214 ctmp = conjg(v(q,p))
1215 v(q,p) = conjg(v(p,q))
1216 v(p,q) = ctmp
1217 1114 CONTINUE
1218 1113 CONTINUE
1219 CALL clapmt( .false., n, n, v, ldv, iwork )
1220* .. assemble the left singular vector matrix U of dimensions
1221* (M x N1), i.e. (M x N) or (M x M).
1222*
1223 DO 1111 p = 1, n
1224 u(p,p) = conjg(u(p,p))
1225 DO 1112 q = p + 1, n
1226 ctmp = conjg(u(q,p))
1227 u(q,p) = conjg(u(p,q))
1228 u(p,q) = ctmp
1229 1112 CONTINUE
1230 1111 CONTINUE
1231*
1232 IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1233 CALL claset('A',m-n,n,czero,czero,u(n+1,1),ldu)
1234 IF ( n .LT. n1 ) THEN
1235 CALL claset('A',n,n1-n,czero,czero,u(1,n+1),
1236 $ ldu)
1237 CALL claset('A',m-n,n1-n,czero,cone,
1238 $ u(n+1,n+1), ldu )
1239 END IF
1240 END IF
1241 ELSE
1242* .. copy R**H into [U] and overwrite [U] with the right
1243* singular vectors of R
1244 DO 1196 p = 1, nr
1245 DO 1197 q = p, n
1246 u(q,nr+p) = conjg(a(p,q))
1247 1197 CONTINUE
1248 1196 CONTINUE
1249 IF ( nr .GT. 1 )
1250 $ CALL claset('U',nr-1,nr-1,czero,czero,u(1,nr+2),
1251 $ ldu)
1252 CALL cgeqrf( n, nr, u(1,nr+1), ldu, cwork(n+1),
1253 $ cwork(n+nr+1), lcwork-n-nr, ierr )
1254 DO 1143 p = 1, nr
1255 DO 1144 q = 1, n
1256 v(q,p) = conjg(u(p,nr+q))
1257 1144 CONTINUE
1258 1143 CONTINUE
1259 CALL claset('U',nr-1,nr-1,czero,czero,v(1,2),ldv)
1260 CALL cgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1261 $ v,ldv, cwork(n+nr+1),lcwork-n-nr,rwork, info )
1262 CALL claset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1263 CALL claset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1264 CALL claset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),
1265 $ ldv)
1266 CALL cunmqr('R','C', n, n, nr, u(1,nr+1), ldu,
1267 $ cwork(n+1),v,ldv,cwork(n+nr+1),lcwork-n-nr,ierr)
1268 CALL clapmt( .false., n, n, v, ldv, iwork )
1269* .. assemble the left singular vector matrix U of dimensions
1270* (M x NR) or (M x N) or (M x M).
1271 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1272 CALL claset('A',m-nr,nr,czero,czero,u(nr+1,1),
1273 $ ldu)
1274 IF ( nr .LT. n1 ) THEN
1275 CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),
1276 $ ldu)
1277 CALL claset( 'A',m-nr,n1-nr,czero,cone,
1278 $ u(nr+1,nr+1),ldu)
1279 END IF
1280 END IF
1281 END IF
1282 END IF
1283*
1284 ELSE
1285*
1286* .. apply CGESVD to R [[this is the recommended option]]
1287*
1288 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1289* .. copy R into [V] and overwrite V with the right singular vectors
1290 CALL clacpy( 'U', nr, n, a, lda, v, ldv )
1291 IF ( nr .GT. 1 )
1292 $ CALL claset( 'L', nr-1,nr-1, czero,czero, v(2,1),
1293 $ ldv )
1294* .. the right singular vectors of R overwrite [V], the NR left
1295* singular vectors of R stored in [U](1:NR,1:NR)
1296 CALL cgesvd( 'S', 'O', nr, n, v, ldv, s, u, ldu,
1297 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1298 CALL clapmt( .false., nr, n, v, ldv, iwork )
1299* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1300* .. assemble the left singular vector matrix U of dimensions
1301* (M x NR) or (M x N) or (M x M).
1302 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1303 CALL claset('A', m-nr,nr, czero,czero, u(nr+1,1),
1304 $ ldu)
1305 IF ( nr .LT. n1 ) THEN
1306 CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),
1307 $ ldu)
1308 CALL claset( 'A',m-nr,n1-nr,czero,cone,
1309 $ u(nr+1,nr+1), ldu )
1310 END IF
1311 END IF
1312*
1313 ELSE
1314* .. need all N right singular vectors and NR < N
1315* .. the requested number of the left singular vectors
1316* is then N1 (N or M)
1317* [[The optimal ratio N/NR for using LQ instead of padding
1318* with zeros. Here hard coded to 2; it must be at least
1319* two due to work space constraints.]]
1320* OPTRATIO = ILAENV(6, 'CGESVD', 'S' // 'O', NR,N,0,0)
1321* OPTRATIO = MAX( OPTRATIO, 2 )
1322 optratio = 2
1323 IF ( optratio * nr .GT. n ) THEN
1324 CALL clacpy( 'U', nr, n, a, lda, v, ldv )
1325 IF ( nr .GT. 1 )
1326 $ CALL claset('L', nr-1,nr-1, czero,czero, v(2,1),
1327 $ ldv)
1328* .. the right singular vectors of R overwrite [V], the NR left
1329* singular vectors of R stored in [U](1:NR,1:NR)
1330 CALL claset('A', n-nr,n, czero,czero, v(nr+1,1),
1331 $ ldv)
1332 CALL cgesvd( 'S', 'O', n, n, v, ldv, s, u, ldu,
1333 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1334 CALL clapmt( .false., n, n, v, ldv, iwork )
1335* .. now [V] contains the adjoint of the matrix of the right
1336* singular vectors of A. The leading N left singular vectors
1337* are in [U](1:N,1:N)
1338* .. assemble the left singular vector matrix U of dimensions
1339* (M x N1), i.e. (M x N) or (M x M).
1340 IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1341 CALL claset('A',m-n,n,czero,czero,u(n+1,1),ldu)
1342 IF ( n .LT. n1 ) THEN
1343 CALL claset('A',n,n1-n,czero,czero,u(1,n+1),
1344 $ ldu)
1345 CALL claset( 'A',m-n,n1-n,czero,cone,
1346 $ u(n+1,n+1), ldu )
1347 END IF
1348 END IF
1349 ELSE
1350 CALL clacpy( 'U', nr, n, a, lda, u(nr+1,1), ldu )
1351 IF ( nr .GT. 1 )
1352 $ CALL claset('L',nr-1,nr-1,czero,czero,u(nr+2,1),
1353 $ ldu)
1354 CALL cgelqf( nr, n, u(nr+1,1), ldu, cwork(n+1),
1355 $ cwork(n+nr+1), lcwork-n-nr, ierr )
1356 CALL clacpy('L',nr,nr,u(nr+1,1),ldu,v,ldv)
1357 IF ( nr .GT. 1 )
1358 $ CALL claset('U',nr-1,nr-1,czero,czero,v(1,2),ldv)
1359 CALL cgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1360 $ v, ldv, cwork(n+nr+1), lcwork-n-nr, rwork, info )
1361 CALL claset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1362 CALL claset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1363 CALL claset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),
1364 $ ldv)
1365 CALL cunmlq('R','N',n,n,nr,u(nr+1,1),ldu,
1366 $ cwork(n+1),
1367 $ v, ldv, cwork(n+nr+1),lcwork-n-nr,ierr)
1368 CALL clapmt( .false., n, n, v, ldv, iwork )
1369* .. assemble the left singular vector matrix U of dimensions
1370* (M x NR) or (M x N) or (M x M).
1371 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1372 CALL claset('A',m-nr,nr,czero,czero,u(nr+1,1),
1373 $ ldu)
1374 IF ( nr .LT. n1 ) THEN
1375 CALL claset('A',nr,n1-nr,czero,czero,u(1,nr+1),
1376 $ ldu)
1377 CALL claset( 'A',m-nr,n1-nr,czero,cone,
1378 $ u(nr+1,nr+1), ldu )
1379 END IF
1380 END IF
1381 END IF
1382 END IF
1383* .. end of the "R**H or R" branch
1384 END IF
1385*
1386* The Q matrix from the first QRF is built into the left singular
1387* vectors matrix U.
1388*
1389 IF ( .NOT. wntuf )
1390 $ CALL cunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1391 $ ldu, cwork(n+1), lcwork-n, ierr )
1392 IF ( rowprm .AND. .NOT.wntuf )
1393 $ CALL claswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1394*
1395* ... end of the "full SVD" branch
1396 END IF
1397*
1398* Check whether some singular values are returned as zeros, e.g.
1399* due to underflow, and update the numerical rank.
1400 p = nr
1401 DO 4001 q = p, 1, -1
1402 IF ( s(q) .GT. zero ) GO TO 4002
1403 nr = nr - 1
1404 4001 CONTINUE
1405 4002 CONTINUE
1406*
1407* .. if numerical rank deficiency is detected, the truncated
1408* singular values are set to zero.
1409 IF ( nr .LT. n ) CALL slaset( 'G', n-nr,1, zero,zero, s(nr+1),
1410 $ n )
1411* .. undo scaling; this may cause overflow in the largest singular
1412* values.
1413 IF ( ascaled )
1414 $ CALL slascl( 'G',0,0, one,sqrt(real(m)), nr,1, s, n, ierr )
1415 IF ( conda ) rwork(1) = sconda
1416 rwork(2) = real( p - nr )
1417* .. p-NR is the number of singular values that are computed as
1418* exact zeros in CGESVD() applied to the (possibly truncated)
1419* full row rank triangular (trapezoidal) factor of A.
1420 numrank = nr
1421*
1422 RETURN
1423*
1424* End of CGESVDQ
1425*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgelqf
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
Definition cgelqf.f:142
cgeqp3
subroutine cgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
CGEQP3
Definition cgeqp3.f:157
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:144
cgesvd
subroutine cgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition cgesvd.f:213
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:113
clapmt
subroutine clapmt(forwrd, m, n, x, ldx, k)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition clapmt.f:102
clascl
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:142
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:104
claswp
subroutine claswp(n, a, lda, k1, k2, ipiv, incx)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition claswp.f:113
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
scnrm2
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90
cpocon
subroutine cpocon(uplo, n, a, lda, anorm, rcond, work, rwork, info)
CPOCON
Definition cpocon.f:119
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
cunmlq
subroutine cunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMLQ
Definition cunmlq.f:166
cunmqr
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:166
Here is the call graph for this function:
Here is the caller graph for this function: