LAPACK 3.12.0
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 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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
Definition cgelqf.f:143
subroutine cgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
CGEQP3
Definition cgeqp3.f:159
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
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
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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 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 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 slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:143
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
subroutine 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 claswp(n, a, lda, k1, k2, ipiv, incx)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition claswp.f:115
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90
subroutine cpocon(uplo, n, a, lda, anorm, rcond, work, rwork, info)
CPOCON
Definition cpocon.f:121
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
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
Here is the call graph for this function:
Here is the caller graph for this function: