LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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, enlosed in [[double brackets]].
Bugs, examples and comments
  Please report all bugs and send interesting examples and/or comments to
  drmac@math.hr. Thank you.
References
  [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
      Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
      44(1): 11:1-11:30 (2017)

  SIGMA library, xGESVDQ section updated February 2016.
  Developed and coded by Zlatko Drmac, Department of Mathematics
  University of Zagreb, Croatia, drmac@math.hr
Contributors:
 Developed and coded by Zlatko Drmac, Department of Mathematics
  University of Zagreb, Croatia, drmac@math.hr
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 410 of file cgesvdq.f.

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