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

◆ zgesvdq()

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

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

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

Purpose:
 ZCGESVDQ 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 = DLAMCH('Epsilon'). This authorises ZGESVDQ 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, ZGESVD 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*16 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 ZGEQP3. 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 DOUBLE PRECISION array of dimension N.
          The singular values of A, ordered so that S(i) >= S(i+1).
[out]U
          U is COMPLEX*16 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*16 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*12 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
          ZGEQP3.

          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 DOUBLE PRECISION array, dimension (max(1, LRWORK)).
          On exit,
          1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
          number of column scaled A. If A = C * D where D is diagonal and C
          has unit columns in the Euclidean norm, then, assuming full column rank,
          N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
          Otherwise, RWORK(1) = -1.
          2. RWORK(2) contains the number of singular values computed as
          exact zeros in ZGESVD applied to the upper triangular or trapezoidal
          R (from the initial QR factorization). In case of early exit (no call to
          ZGESVD, 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 ZBDSQR did not converge, INFO specifies how many superdiagonals
          of an intermediate bidiagonal form B (computed in ZGESVD) 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 zgesvdq.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*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
421 DOUBLE PRECISION S( * ), RWORK( * )
422 INTEGER IWORK( * )
423*
424* =====================================================================
425*
426* .. Parameters ..
427 DOUBLE PRECISION ZERO, ONE
428 parameter( zero = 0.0d0, one = 1.0d0 )
429 COMPLEX*16 CZERO, CONE
430 parameter( czero = (0.0d0,0.0d0), cone = (1.0d0,0.0d0) )
431* ..
432* .. Local Scalars ..
433 INTEGER IERR, NR, N1, OPTRATIO, p, q
434 INTEGER LWCON, LWQP3, LWRK_ZGELQF, LWRK_ZGESVD, LWRK_ZGESVD2,
435 $ LWRK_ZGEQP3, LWRK_ZGEQRF, LWRK_ZUNMLQ, LWRK_ZUNMQR,
436 $ LWRK_ZUNMQR2, 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 DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN
443 COMPLEX*16 CTMP
444* ..
445* .. Local Arrays
446 COMPLEX*16 CDUMMY(1)
447 DOUBLE PRECISION RDUMMY(1)
448* ..
449* .. External Subroutines (BLAS, LAPACK)
450 EXTERNAL zgelqf, zgeqp3, zgeqrf, zgesvd, zlacpy, zlapmt,
453* ..
454* .. External Functions (BLAS, LAPACK)
455 LOGICAL LSAME
456 INTEGER IDAMAX
457 DOUBLE PRECISION ZLANGE, DZNRM2, DLAMCH
458 EXTERNAL lsame, zlange, idamax, dznrm2, dlamch
459* ..
460* .. Intrinsic Functions ..
461 INTRINSIC abs, conjg, max, min, dble, 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* .. compute the minimal and the optimal workspace lengths
528* [[The expressions for computing the minimal and the optimal
529* values of LCWORK are written with a lot of redundancy and
530* can be simplified. However, this detailed form is easier for
531* maintenance and modifications of the code.]]
532*
533* .. minimal workspace length for ZGEQP3 of an M x N matrix
534 lwqp3 = n+1
535* .. minimal workspace length for ZUNMQR to build left singular vectors
536 IF ( wntus .OR. wntur ) THEN
537 lwunq = max( n , 1 )
538 ELSE IF ( wntua ) THEN
539 lwunq = max( m , 1 )
540 END IF
541* .. minimal workspace length for ZPOCON of an N x N matrix
542 lwcon = 2 * n
543* .. ZGESVD of an N x N matrix
544 lwsvd = max( 3 * n, 1 )
545 IF ( lquery ) THEN
546 CALL zgeqp3( m, n, a, lda, iwork, cdummy, cdummy, -1,
547 $ rdummy, ierr )
548 lwrk_zgeqp3 = int( cdummy(1) )
549 IF ( wntus .OR. wntur ) THEN
550 CALL zunmqr( 'L', 'N', m, n, n, a, lda, cdummy, u,
551 $ ldu, cdummy, -1, ierr )
552 lwrk_zunmqr = int( cdummy(1) )
553 ELSE IF ( wntua ) THEN
554 CALL zunmqr( 'L', 'N', m, m, n, a, lda, cdummy, u,
555 $ ldu, cdummy, -1, ierr )
556 lwrk_zunmqr = int( cdummy(1) )
557 ELSE
558 lwrk_zunmqr = 0
559 END IF
560 END IF
561 minwrk = 2
562 optwrk = 2
563 IF ( .NOT. (lsvec .OR. rsvec ) ) THEN
564* .. minimal and optimal sizes of the complex workspace if
565* only the singular values are requested
566 IF ( conda ) THEN
567 minwrk = max( n+lwqp3, lwcon, lwsvd )
568 ELSE
569 minwrk = max( n+lwqp3, lwsvd )
570 END IF
571 IF ( lquery ) THEN
572 CALL zgesvd( 'N', 'N', n, n, a, lda, s, u, ldu,
573 $ v, ldv, cdummy, -1, rdummy, ierr )
574 lwrk_zgesvd = int( cdummy(1) )
575 IF ( conda ) THEN
576 optwrk = max( n+lwrk_zgeqp3, n+lwcon, lwrk_zgesvd )
577 ELSE
578 optwrk = max( n+lwrk_zgeqp3, lwrk_zgesvd )
579 END IF
580 END IF
581 ELSE IF ( lsvec .AND. (.NOT.rsvec) ) THEN
582* .. minimal and optimal sizes of the complex workspace if the
583* singular values and the left singular vectors are requested
584 IF ( conda ) THEN
585 minwrk = n + max( lwqp3, lwcon, lwsvd, lwunq )
586 ELSE
587 minwrk = n + max( lwqp3, lwsvd, lwunq )
588 END IF
589 IF ( lquery ) THEN
590 IF ( rtrans ) THEN
591 CALL zgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
592 $ v, ldv, cdummy, -1, rdummy, ierr )
593 ELSE
594 CALL zgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
595 $ v, ldv, cdummy, -1, rdummy, ierr )
596 END IF
597 lwrk_zgesvd = int( cdummy(1) )
598 IF ( conda ) THEN
599 optwrk = n + max( lwrk_zgeqp3, lwcon, lwrk_zgesvd,
600 $ lwrk_zunmqr )
601 ELSE
602 optwrk = n + max( lwrk_zgeqp3, lwrk_zgesvd,
603 $ lwrk_zunmqr )
604 END IF
605 END IF
606 ELSE IF ( rsvec .AND. (.NOT.lsvec) ) THEN
607* .. minimal and optimal sizes of the complex workspace if the
608* singular values and the right singular vectors are requested
609 IF ( conda ) THEN
610 minwrk = n + max( lwqp3, lwcon, lwsvd )
611 ELSE
612 minwrk = n + max( lwqp3, lwsvd )
613 END IF
614 IF ( lquery ) THEN
615 IF ( rtrans ) THEN
616 CALL zgesvd( 'O', 'N', n, n, a, lda, s, u, ldu,
617 $ v, ldv, cdummy, -1, rdummy, ierr )
618 ELSE
619 CALL zgesvd( 'N', 'O', n, n, a, lda, s, u, ldu,
620 $ v, ldv, cdummy, -1, rdummy, ierr )
621 END IF
622 lwrk_zgesvd = int( cdummy(1) )
623 IF ( conda ) THEN
624 optwrk = n + max( lwrk_zgeqp3, lwcon, lwrk_zgesvd )
625 ELSE
626 optwrk = n + max( lwrk_zgeqp3, lwrk_zgesvd )
627 END IF
628 END IF
629 ELSE
630* .. minimal and optimal sizes of the complex workspace if the
631* full SVD is requested
632 IF ( rtrans ) THEN
633 minwrk = max( lwqp3, lwsvd, lwunq )
634 IF ( conda ) minwrk = max( minwrk, lwcon )
635 minwrk = minwrk + n
636 IF ( wntva ) THEN
637* .. minimal workspace length for N x N/2 ZGEQRF
638 lwqrf = max( n/2, 1 )
639* .. minimal workspace length for N/2 x N/2 ZGESVD
640 lwsvd2 = max( 3 * (n/2), 1 )
641 lwunq2 = max( n, 1 )
642 minwrk2 = max( lwqp3, n/2+lwqrf, n/2+lwsvd2,
643 $ n/2+lwunq2, lwunq )
644 IF ( conda ) minwrk2 = max( minwrk2, lwcon )
645 minwrk2 = n + minwrk2
646 minwrk = max( minwrk, minwrk2 )
647 END IF
648 ELSE
649 minwrk = max( lwqp3, lwsvd, lwunq )
650 IF ( conda ) minwrk = max( minwrk, lwcon )
651 minwrk = minwrk + n
652 IF ( wntva ) THEN
653* .. minimal workspace length for N/2 x N ZGELQF
654 lwlqf = max( n/2, 1 )
655 lwsvd2 = max( 3 * (n/2), 1 )
656 lwunlq = max( n , 1 )
657 minwrk2 = max( lwqp3, n/2+lwlqf, n/2+lwsvd2,
658 $ n/2+lwunlq, lwunq )
659 IF ( conda ) minwrk2 = max( minwrk2, lwcon )
660 minwrk2 = n + minwrk2
661 minwrk = max( minwrk, minwrk2 )
662 END IF
663 END IF
664 IF ( lquery ) THEN
665 IF ( rtrans ) THEN
666 CALL zgesvd( 'O', 'A', n, n, a, lda, s, u, ldu,
667 $ v, ldv, cdummy, -1, rdummy, ierr )
668 lwrk_zgesvd = int( cdummy(1) )
669 optwrk = max(lwrk_zgeqp3,lwrk_zgesvd,lwrk_zunmqr)
670 IF ( conda ) optwrk = max( optwrk, lwcon )
671 optwrk = n + optwrk
672 IF ( wntva ) THEN
673 CALL zgeqrf(n,n/2,u,ldu,cdummy,cdummy,-1,ierr)
674 lwrk_zgeqrf = int( cdummy(1) )
675 CALL zgesvd( 'S', 'O', n/2,n/2, v,ldv, s, u,ldu,
676 $ v, ldv, cdummy, -1, rdummy, ierr )
677 lwrk_zgesvd2 = int( cdummy(1) )
678 CALL zunmqr( 'R', 'C', n, n, n/2, u, ldu, cdummy,
679 $ v, ldv, cdummy, -1, ierr )
680 lwrk_zunmqr2 = int( cdummy(1) )
681 optwrk2 = max( lwrk_zgeqp3, n/2+lwrk_zgeqrf,
682 $ n/2+lwrk_zgesvd2, n/2+lwrk_zunmqr2 )
683 IF ( conda ) optwrk2 = max( optwrk2, lwcon )
684 optwrk2 = n + optwrk2
685 optwrk = max( optwrk, optwrk2 )
686 END IF
687 ELSE
688 CALL zgesvd( 'S', 'O', n, n, a, lda, s, u, ldu,
689 $ v, ldv, cdummy, -1, rdummy, ierr )
690 lwrk_zgesvd = int( cdummy(1) )
691 optwrk = max(lwrk_zgeqp3,lwrk_zgesvd,lwrk_zunmqr)
692 IF ( conda ) optwrk = max( optwrk, lwcon )
693 optwrk = n + optwrk
694 IF ( wntva ) THEN
695 CALL zgelqf(n/2,n,u,ldu,cdummy,cdummy,-1,ierr)
696 lwrk_zgelqf = int( cdummy(1) )
697 CALL zgesvd( 'S','O', n/2,n/2, v, ldv, s, u, ldu,
698 $ v, ldv, cdummy, -1, rdummy, ierr )
699 lwrk_zgesvd2 = int( cdummy(1) )
700 CALL zunmlq( 'R', 'N', n, n, n/2, u, ldu, cdummy,
701 $ v, ldv, cdummy,-1,ierr )
702 lwrk_zunmlq = int( cdummy(1) )
703 optwrk2 = max( lwrk_zgeqp3, n/2+lwrk_zgelqf,
704 $ n/2+lwrk_zgesvd2, n/2+lwrk_zunmlq )
705 IF ( conda ) optwrk2 = max( optwrk2, lwcon )
706 optwrk2 = n + optwrk2
707 optwrk = max( optwrk, optwrk2 )
708 END IF
709 END IF
710 END IF
711 END IF
712*
713 minwrk = max( 2, minwrk )
714 optwrk = max( 2, optwrk )
715 IF ( lcwork .LT. minwrk .AND. (.NOT.lquery) ) info = -19
716*
717 END IF
718*
719 IF (info .EQ. 0 .AND. lrwork .LT. rminwrk .AND. .NOT. lquery) THEN
720 info = -21
721 END IF
722 IF( info.NE.0 ) THEN
723 CALL xerbla( 'ZGESVDQ', -info )
724 RETURN
725 ELSE IF ( lquery ) THEN
726*
727* Return optimal workspace
728*
729 iwork(1) = iminwrk
730 cwork(1) = optwrk
731 cwork(2) = minwrk
732 rwork(1) = rminwrk
733 RETURN
734 END IF
735*
736* Quick return if the matrix is void.
737*
738 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) ) THEN
739* .. all output is void.
740 RETURN
741 END IF
742*
743 big = dlamch('O')
744 ascaled = .false.
745 IF ( rowprm ) THEN
746* .. reordering the rows in decreasing sequence in the
747* ell-infinity norm - this enhances numerical robustness in
748* the case of differently scaled rows.
749 DO 1904 p = 1, m
750* RWORK(p) = ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
751* [[ZLANGE will return NaN if an entry of the p-th row is Nan]]
752 rwork(p) = zlange( 'M', 1, n, a(p,1), lda, rdummy )
753* .. check for NaN's and Inf's
754 IF ( ( rwork(p) .NE. rwork(p) ) .OR.
755 $ ( (rwork(p)*zero) .NE. zero ) ) THEN
756 info = -8
757 CALL xerbla( 'ZGESVDQ', -info )
758 RETURN
759 END IF
760 1904 CONTINUE
761 DO 1952 p = 1, m - 1
762 q = idamax( m-p+1, rwork(p), 1 ) + p - 1
763 iwork(n+p) = q
764 IF ( p .NE. q ) THEN
765 rtmp = rwork(p)
766 rwork(p) = rwork(q)
767 rwork(q) = rtmp
768 END IF
769 1952 CONTINUE
770*
771 IF ( rwork(1) .EQ. zero ) THEN
772* Quick return: A is the M x N zero matrix.
773 numrank = 0
774 CALL dlaset( 'G', n, 1, zero, zero, s, n )
775 IF ( wntus ) CALL zlaset('G', m, n, czero, cone, u, ldu)
776 IF ( wntua ) CALL zlaset('G', m, m, czero, cone, u, ldu)
777 IF ( wntva ) CALL zlaset('G', n, n, czero, cone, v, ldv)
778 IF ( wntuf ) THEN
779 CALL zlaset( 'G', n, 1, czero, czero, cwork, n )
780 CALL zlaset( 'G', m, n, czero, cone, u, ldu )
781 END IF
782 DO 5001 p = 1, n
783 iwork(p) = p
784 5001 CONTINUE
785 IF ( rowprm ) THEN
786 DO 5002 p = n + 1, n + m - 1
787 iwork(p) = p - n
788 5002 CONTINUE
789 END IF
790 IF ( conda ) rwork(1) = -1
791 rwork(2) = -1
792 RETURN
793 END IF
794*
795 IF ( rwork(1) .GT. big / sqrt(dble(m)) ) THEN
796* .. to prevent overflow in the QR factorization, scale the
797* matrix by 1/sqrt(M) if too large entry detected
798 CALL zlascl('G',0,0,sqrt(dble(m)),one, m,n, a,lda, ierr)
799 ascaled = .true.
800 END IF
801 CALL zlaswp( n, a, lda, 1, m-1, iwork(n+1), 1 )
802 END IF
803*
804* .. At this stage, preemptive scaling is done only to avoid column
805* norms overflows during the QR factorization. The SVD procedure should
806* have its own scaling to save the singular values from overflows and
807* underflows. That depends on the SVD procedure.
808*
809 IF ( .NOT.rowprm ) THEN
810 rtmp = zlange( 'M', m, n, a, lda, rwork )
811 IF ( ( rtmp .NE. rtmp ) .OR.
812 $ ( (rtmp*zero) .NE. zero ) ) THEN
813 info = -8
814 CALL xerbla( 'ZGESVDQ', -info )
815 RETURN
816 END IF
817 IF ( rtmp .GT. big / sqrt(dble(m)) ) THEN
818* .. to prevent overflow in the QR factorization, scale the
819* matrix by 1/sqrt(M) if too large entry detected
820 CALL zlascl('G',0,0, sqrt(dble(m)),one, m,n, a,lda, ierr)
821 ascaled = .true.
822 END IF
823 END IF
824*
825* .. QR factorization with column pivoting
826*
827* A * P = Q * [ R ]
828* [ 0 ]
829*
830 DO 1963 p = 1, n
831* .. all columns are free columns
832 iwork(p) = 0
833 1963 CONTINUE
834 CALL zgeqp3( m, n, a, lda, iwork, cwork, cwork(n+1), lcwork-n,
835 $ rwork, ierr )
836*
837* If the user requested accuracy level allows truncation in the
838* computed upper triangular factor, the matrix R is examined and,
839* if possible, replaced with its leading upper trapezoidal part.
840*
841 epsln = dlamch('E')
842 sfmin = dlamch('S')
843* SMALL = SFMIN / EPSLN
844 nr = n
845*
846 IF ( accla ) THEN
847*
848* Standard absolute error bound suffices. All sigma_i with
849* sigma_i < N*EPS*||A||_F are flushed to zero. This is an
850* aggressive enforcement of lower numerical rank by introducing a
851* backward error of the order of N*EPS*||A||_F.
852 nr = 1
853 rtmp = sqrt(dble(n))*epsln
854 DO 3001 p = 2, n
855 IF ( abs(a(p,p)) .LT. (rtmp*abs(a(1,1))) ) GO TO 3002
856 nr = nr + 1
857 3001 CONTINUE
858 3002 CONTINUE
859*
860 ELSEIF ( acclm ) THEN
861* .. similarly as above, only slightly more gentle (less aggressive).
862* Sudden drop on the diagonal of R is used as the criterion for being
863* close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E').
864* [[This can be made more flexible by replacing this hard-coded value
865* with a user specified threshold.]] Also, the values that underflow
866* will be truncated.
867 nr = 1
868 DO 3401 p = 2, n
869 IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
870 $ ( abs(a(p,p)) .LT. sfmin ) ) GO TO 3402
871 nr = nr + 1
872 3401 CONTINUE
873 3402 CONTINUE
874*
875 ELSE
876* .. RRQR not authorized to determine numerical rank except in the
877* obvious case of zero pivots.
878* .. inspect R for exact zeros on the diagonal;
879* R(i,i)=0 => R(i:N,i:N)=0.
880 nr = 1
881 DO 3501 p = 2, n
882 IF ( abs(a(p,p)) .EQ. zero ) GO TO 3502
883 nr = nr + 1
884 3501 CONTINUE
885 3502 CONTINUE
886*
887 IF ( conda ) THEN
888* Estimate the scaled condition number of A. Use the fact that it is
889* the same as the scaled condition number of R.
890* .. V is used as workspace
891 CALL zlacpy( 'U', n, n, a, lda, v, ldv )
892* Only the leading NR x NR submatrix of the triangular factor
893* is considered. Only if NR=N will this give a reliable error
894* bound. However, even for NR < N, this can be used on an
895* expert level and obtain useful information in the sense of
896* perturbation theory.
897 DO 3053 p = 1, nr
898 rtmp = dznrm2( p, v(1,p), 1 )
899 CALL zdscal( p, one/rtmp, v(1,p), 1 )
900 3053 CONTINUE
901 IF ( .NOT. ( lsvec .OR. rsvec ) ) THEN
902 CALL zpocon( 'U', nr, v, ldv, one, rtmp,
903 $ cwork, rwork, ierr )
904 ELSE
905 CALL zpocon( 'U', nr, v, ldv, one, rtmp,
906 $ cwork(n+1), rwork, ierr )
907 END IF
908 sconda = one / sqrt(rtmp)
909* For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
910* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
911* See the reference [1] for more details.
912 END IF
913*
914 ENDIF
915*
916 IF ( wntur ) THEN
917 n1 = nr
918 ELSE IF ( wntus .OR. wntuf) THEN
919 n1 = n
920 ELSE IF ( wntua ) THEN
921 n1 = m
922 END IF
923*
924 IF ( .NOT. ( rsvec .OR. lsvec ) ) THEN
925*.......................................................................
926* .. only the singular values are requested
927*.......................................................................
928 IF ( rtrans ) THEN
929*
930* .. compute the singular values of R**H = [A](1:NR,1:N)**H
931* .. set the lower triangle of [A] to [A](1:NR,1:N)**H and
932* the upper triangle of [A] to zero.
933 DO 1146 p = 1, min( n, nr )
934 a(p,p) = conjg(a(p,p))
935 DO 1147 q = p + 1, n
936 a(q,p) = conjg(a(p,q))
937 IF ( q .LE. nr ) a(p,q) = czero
938 1147 CONTINUE
939 1146 CONTINUE
940*
941 CALL zgesvd( 'N', 'N', n, nr, a, lda, s, u, ldu,
942 $ v, ldv, cwork, lcwork, rwork, info )
943*
944 ELSE
945*
946* .. compute the singular values of R = [A](1:NR,1:N)
947*
948 IF ( nr .GT. 1 )
949 $ CALL zlaset( 'L', nr-1,nr-1, czero,czero, a(2,1), lda )
950 CALL zgesvd( 'N', 'N', nr, n, a, lda, s, u, ldu,
951 $ v, ldv, cwork, lcwork, rwork, info )
952*
953 END IF
954*
955 ELSE IF ( lsvec .AND. ( .NOT. rsvec) ) THEN
956*.......................................................................
957* .. the singular values and the left singular vectors requested
958*.......................................................................""""""""
959 IF ( rtrans ) THEN
960* .. apply ZGESVD to R**H
961* .. copy R**H into [U] and overwrite [U] with the right singular
962* vectors of R
963 DO 1192 p = 1, nr
964 DO 1193 q = p, n
965 u(q,p) = conjg(a(p,q))
966 1193 CONTINUE
967 1192 CONTINUE
968 IF ( nr .GT. 1 )
969 $ CALL zlaset( 'U', nr-1,nr-1, czero,czero, u(1,2), ldu )
970* .. the left singular vectors not computed, the NR right singular
971* vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These
972* will be pre-multiplied by Q to build the left singular vectors of A.
973 CALL zgesvd( 'N', 'O', n, nr, u, ldu, s, u, ldu,
974 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
975*
976 DO 1119 p = 1, nr
977 u(p,p) = conjg(u(p,p))
978 DO 1120 q = p + 1, nr
979 ctmp = conjg(u(q,p))
980 u(q,p) = conjg(u(p,q))
981 u(p,q) = ctmp
982 1120 CONTINUE
983 1119 CONTINUE
984*
985 ELSE
986* .. apply ZGESVD to R
987* .. copy R into [U] and overwrite [U] with the left singular vectors
988 CALL zlacpy( 'U', nr, n, a, lda, u, ldu )
989 IF ( nr .GT. 1 )
990 $ CALL zlaset( 'L', nr-1, nr-1, czero, czero, u(2,1), ldu )
991* .. the right singular vectors not computed, the NR left singular
992* vectors overwrite [U](1:NR,1:NR)
993 CALL zgesvd( 'O', 'N', nr, n, u, ldu, s, u, ldu,
994 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
995* .. now [U](1:NR,1:NR) contains the NR left singular vectors of
996* R. These will be pre-multiplied by Q to build the left singular
997* vectors of A.
998 END IF
999*
1000* .. assemble the left singular vector matrix U of dimensions
1001* (M x NR) or (M x N) or (M x M).
1002 IF ( ( nr .LT. m ) .AND. ( .NOT.wntuf ) ) THEN
1003 CALL zlaset('A', m-nr, nr, czero, czero, u(nr+1,1), ldu)
1004 IF ( nr .LT. n1 ) THEN
1005 CALL zlaset( 'A',nr,n1-nr,czero,czero,u(1,nr+1), ldu )
1006 CALL zlaset( 'A',m-nr,n1-nr,czero,cone,
1007 $ u(nr+1,nr+1), ldu )
1008 END IF
1009 END IF
1010*
1011* The Q matrix from the first QRF is built into the left singular
1012* vectors matrix U.
1013*
1014 IF ( .NOT.wntuf )
1015 $ CALL zunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1016 $ ldu, cwork(n+1), lcwork-n, ierr )
1017 IF ( rowprm .AND. .NOT.wntuf )
1018 $ CALL zlaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1019*
1020 ELSE IF ( rsvec .AND. ( .NOT. lsvec ) ) THEN
1021*.......................................................................
1022* .. the singular values and the right singular vectors requested
1023*.......................................................................
1024 IF ( rtrans ) THEN
1025* .. apply ZGESVD to R**H
1026* .. copy R**H into V and overwrite V with the left singular vectors
1027 DO 1165 p = 1, nr
1028 DO 1166 q = p, n
1029 v(q,p) = conjg(a(p,q))
1030 1166 CONTINUE
1031 1165 CONTINUE
1032 IF ( nr .GT. 1 )
1033 $ CALL zlaset( 'U', nr-1,nr-1, czero,czero, v(1,2), ldv )
1034* .. the left singular vectors of R**H overwrite V, the right singular
1035* vectors not computed
1036 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1037 CALL zgesvd( 'O', 'N', n, nr, v, ldv, s, u, ldu,
1038 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1039*
1040 DO 1121 p = 1, nr
1041 v(p,p) = conjg(v(p,p))
1042 DO 1122 q = p + 1, nr
1043 ctmp = conjg(v(q,p))
1044 v(q,p) = conjg(v(p,q))
1045 v(p,q) = ctmp
1046 1122 CONTINUE
1047 1121 CONTINUE
1048*
1049 IF ( nr .LT. n ) THEN
1050 DO 1103 p = 1, nr
1051 DO 1104 q = nr + 1, n
1052 v(p,q) = conjg(v(q,p))
1053 1104 CONTINUE
1054 1103 CONTINUE
1055 END IF
1056 CALL zlapmt( .false., nr, n, v, ldv, iwork )
1057 ELSE
1058* .. need all N right singular vectors and NR < N
1059* [!] This is simple implementation that augments [V](1:N,1:NR)
1060* by padding a zero block. In the case NR << N, a more efficient
1061* way is to first use the QR factorization. For more details
1062* how to implement this, see the " FULL SVD " branch.
1063 CALL zlaset('G', n, n-nr, czero, czero, v(1,nr+1), ldv)
1064 CALL zgesvd( 'O', 'N', n, n, v, ldv, s, u, ldu,
1065 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1066*
1067 DO 1123 p = 1, n
1068 v(p,p) = conjg(v(p,p))
1069 DO 1124 q = p + 1, n
1070 ctmp = conjg(v(q,p))
1071 v(q,p) = conjg(v(p,q))
1072 v(p,q) = ctmp
1073 1124 CONTINUE
1074 1123 CONTINUE
1075 CALL zlapmt( .false., n, n, v, ldv, iwork )
1076 END IF
1077*
1078 ELSE
1079* .. aply ZGESVD to R
1080* .. copy R into V and overwrite V with the right singular vectors
1081 CALL zlacpy( 'U', nr, n, a, lda, v, ldv )
1082 IF ( nr .GT. 1 )
1083 $ CALL zlaset( 'L', nr-1, nr-1, czero, czero, v(2,1), ldv )
1084* .. the right singular vectors overwrite V, the NR left singular
1085* vectors stored in U(1:NR,1:NR)
1086 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1087 CALL zgesvd( 'N', 'O', nr, n, v, ldv, s, u, ldu,
1088 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1089 CALL zlapmt( .false., nr, n, v, ldv, iwork )
1090* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1091 ELSE
1092* .. need all N right singular vectors and NR < N
1093* [!] This is simple implementation that augments [V](1:NR,1:N)
1094* by padding a zero block. In the case NR << N, a more efficient
1095* way is to first use the LQ factorization. For more details
1096* how to implement this, see the " FULL SVD " branch.
1097 CALL zlaset('G', n-nr, n, czero,czero, v(nr+1,1), ldv)
1098 CALL zgesvd( 'N', 'O', n, n, v, ldv, s, u, ldu,
1099 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1100 CALL zlapmt( .false., n, n, v, ldv, iwork )
1101 END IF
1102* .. now [V] contains the adjoint of the matrix of the right singular
1103* vectors of A.
1104 END IF
1105*
1106 ELSE
1107*.......................................................................
1108* .. FULL SVD requested
1109*.......................................................................
1110 IF ( rtrans ) THEN
1111*
1112* .. apply ZGESVD to R**H [[this option is left for R&D&T]]
1113*
1114 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1115* .. copy R**H into [V] and overwrite [V] with the left singular
1116* vectors of R**H
1117 DO 1168 p = 1, nr
1118 DO 1169 q = p, n
1119 v(q,p) = conjg(a(p,q))
1120 1169 CONTINUE
1121 1168 CONTINUE
1122 IF ( nr .GT. 1 )
1123 $ CALL zlaset( 'U', nr-1,nr-1, czero,czero, v(1,2), ldv )
1124*
1125* .. the left singular vectors of R**H overwrite [V], the NR right
1126* singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate
1127* transposed
1128 CALL zgesvd( 'O', 'A', n, nr, v, ldv, s, v, ldv,
1129 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1130* .. assemble V
1131 DO 1115 p = 1, nr
1132 v(p,p) = conjg(v(p,p))
1133 DO 1116 q = p + 1, nr
1134 ctmp = conjg(v(q,p))
1135 v(q,p) = conjg(v(p,q))
1136 v(p,q) = ctmp
1137 1116 CONTINUE
1138 1115 CONTINUE
1139 IF ( nr .LT. n ) THEN
1140 DO 1101 p = 1, nr
1141 DO 1102 q = nr+1, n
1142 v(p,q) = conjg(v(q,p))
1143 1102 CONTINUE
1144 1101 CONTINUE
1145 END IF
1146 CALL zlapmt( .false., nr, n, v, ldv, iwork )
1147*
1148 DO 1117 p = 1, nr
1149 u(p,p) = conjg(u(p,p))
1150 DO 1118 q = p + 1, nr
1151 ctmp = conjg(u(q,p))
1152 u(q,p) = conjg(u(p,q))
1153 u(p,q) = ctmp
1154 1118 CONTINUE
1155 1117 CONTINUE
1156*
1157 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1158 CALL zlaset('A', m-nr,nr, czero,czero, u(nr+1,1), ldu)
1159 IF ( nr .LT. n1 ) THEN
1160 CALL zlaset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1161 CALL zlaset( 'A',m-nr,n1-nr,czero,cone,
1162 $ u(nr+1,nr+1), ldu )
1163 END IF
1164 END IF
1165*
1166 ELSE
1167* .. need all N right singular vectors and NR < N
1168* .. copy R**H into [V] and overwrite [V] with the left singular
1169* vectors of R**H
1170* [[The optimal ratio N/NR for using QRF instead of padding
1171* with zeros. Here hard coded to 2; it must be at least
1172* two due to work space constraints.]]
1173* OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0)
1174* OPTRATIO = MAX( OPTRATIO, 2 )
1175 optratio = 2
1176 IF ( optratio*nr .GT. n ) THEN
1177 DO 1198 p = 1, nr
1178 DO 1199 q = p, n
1179 v(q,p) = conjg(a(p,q))
1180 1199 CONTINUE
1181 1198 CONTINUE
1182 IF ( nr .GT. 1 )
1183 $ CALL zlaset('U',nr-1,nr-1, czero,czero, v(1,2),ldv)
1184*
1185 CALL zlaset('A',n,n-nr,czero,czero,v(1,nr+1),ldv)
1186 CALL zgesvd( 'O', 'A', n, n, v, ldv, s, v, ldv,
1187 $ u, ldu, cwork(n+1), lcwork-n, rwork, info )
1188*
1189 DO 1113 p = 1, n
1190 v(p,p) = conjg(v(p,p))
1191 DO 1114 q = p + 1, n
1192 ctmp = conjg(v(q,p))
1193 v(q,p) = conjg(v(p,q))
1194 v(p,q) = ctmp
1195 1114 CONTINUE
1196 1113 CONTINUE
1197 CALL zlapmt( .false., n, n, v, ldv, iwork )
1198* .. assemble the left singular vector matrix U of dimensions
1199* (M x N1), i.e. (M x N) or (M x M).
1200*
1201 DO 1111 p = 1, n
1202 u(p,p) = conjg(u(p,p))
1203 DO 1112 q = p + 1, n
1204 ctmp = conjg(u(q,p))
1205 u(q,p) = conjg(u(p,q))
1206 u(p,q) = ctmp
1207 1112 CONTINUE
1208 1111 CONTINUE
1209*
1210 IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1211 CALL zlaset('A',m-n,n,czero,czero,u(n+1,1),ldu)
1212 IF ( n .LT. n1 ) THEN
1213 CALL zlaset('A',n,n1-n,czero,czero,u(1,n+1),ldu)
1214 CALL zlaset('A',m-n,n1-n,czero,cone,
1215 $ u(n+1,n+1), ldu )
1216 END IF
1217 END IF
1218 ELSE
1219* .. copy R**H into [U] and overwrite [U] with the right
1220* singular vectors of R
1221 DO 1196 p = 1, nr
1222 DO 1197 q = p, n
1223 u(q,nr+p) = conjg(a(p,q))
1224 1197 CONTINUE
1225 1196 CONTINUE
1226 IF ( nr .GT. 1 )
1227 $ CALL zlaset('U',nr-1,nr-1,czero,czero,u(1,nr+2),ldu)
1228 CALL zgeqrf( n, nr, u(1,nr+1), ldu, cwork(n+1),
1229 $ cwork(n+nr+1), lcwork-n-nr, ierr )
1230 DO 1143 p = 1, nr
1231 DO 1144 q = 1, n
1232 v(q,p) = conjg(u(p,nr+q))
1233 1144 CONTINUE
1234 1143 CONTINUE
1235 CALL zlaset('U',nr-1,nr-1,czero,czero,v(1,2),ldv)
1236 CALL zgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1237 $ v,ldv, cwork(n+nr+1),lcwork-n-nr,rwork, info )
1238 CALL zlaset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1239 CALL zlaset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1240 CALL zlaset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1241 CALL zunmqr('R','C', n, n, nr, u(1,nr+1), ldu,
1242 $ cwork(n+1),v,ldv,cwork(n+nr+1),lcwork-n-nr,ierr)
1243 CALL zlapmt( .false., n, n, v, ldv, iwork )
1244* .. assemble the left singular vector matrix U of dimensions
1245* (M x NR) or (M x N) or (M x M).
1246 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1247 CALL zlaset('A',m-nr,nr,czero,czero,u(nr+1,1),ldu)
1248 IF ( nr .LT. n1 ) THEN
1249 CALL zlaset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1250 CALL zlaset( 'A',m-nr,n1-nr,czero,cone,
1251 $ u(nr+1,nr+1),ldu)
1252 END IF
1253 END IF
1254 END IF
1255 END IF
1256*
1257 ELSE
1258*
1259* .. apply ZGESVD to R [[this is the recommended option]]
1260*
1261 IF ( wntvr .OR. ( nr .EQ. n ) ) THEN
1262* .. copy R into [V] and overwrite V with the right singular vectors
1263 CALL zlacpy( 'U', nr, n, a, lda, v, ldv )
1264 IF ( nr .GT. 1 )
1265 $ CALL zlaset( 'L', nr-1,nr-1, czero,czero, v(2,1), ldv )
1266* .. the right singular vectors of R overwrite [V], the NR left
1267* singular vectors of R stored in [U](1:NR,1:NR)
1268 CALL zgesvd( 'S', 'O', nr, n, v, ldv, s, u, ldu,
1269 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1270 CALL zlapmt( .false., nr, n, v, ldv, iwork )
1271* .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
1272* .. assemble the left singular vector matrix U of dimensions
1273* (M x NR) or (M x N) or (M x M).
1274 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1275 CALL zlaset('A', m-nr,nr, czero,czero, u(nr+1,1), ldu)
1276 IF ( nr .LT. n1 ) THEN
1277 CALL zlaset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1278 CALL zlaset( 'A',m-nr,n1-nr,czero,cone,
1279 $ u(nr+1,nr+1), ldu )
1280 END IF
1281 END IF
1282*
1283 ELSE
1284* .. need all N right singular vectors and NR < N
1285* .. the requested number of the left singular vectors
1286* is then N1 (N or M)
1287* [[The optimal ratio N/NR for using LQ instead of padding
1288* with zeros. Here hard coded to 2; it must be at least
1289* two due to work space constraints.]]
1290* OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0)
1291* OPTRATIO = MAX( OPTRATIO, 2 )
1292 optratio = 2
1293 IF ( optratio * nr .GT. n ) THEN
1294 CALL zlacpy( 'U', nr, n, a, lda, v, ldv )
1295 IF ( nr .GT. 1 )
1296 $ CALL zlaset('L', nr-1,nr-1, czero,czero, v(2,1),ldv)
1297* .. the right singular vectors of R overwrite [V], the NR left
1298* singular vectors of R stored in [U](1:NR,1:NR)
1299 CALL zlaset('A', n-nr,n, czero,czero, v(nr+1,1),ldv)
1300 CALL zgesvd( 'S', 'O', n, n, v, ldv, s, u, ldu,
1301 $ v, ldv, cwork(n+1), lcwork-n, rwork, info )
1302 CALL zlapmt( .false., n, n, v, ldv, iwork )
1303* .. now [V] contains the adjoint of the matrix of the right
1304* singular vectors of A. The leading N left singular vectors
1305* are in [U](1:N,1:N)
1306* .. assemble the left singular vector matrix U of dimensions
1307* (M x N1), i.e. (M x N) or (M x M).
1308 IF ( ( n .LT. m ) .AND. .NOT.(wntuf)) THEN
1309 CALL zlaset('A',m-n,n,czero,czero,u(n+1,1),ldu)
1310 IF ( n .LT. n1 ) THEN
1311 CALL zlaset('A',n,n1-n,czero,czero,u(1,n+1),ldu)
1312 CALL zlaset( 'A',m-n,n1-n,czero,cone,
1313 $ u(n+1,n+1), ldu )
1314 END IF
1315 END IF
1316 ELSE
1317 CALL zlacpy( 'U', nr, n, a, lda, u(nr+1,1), ldu )
1318 IF ( nr .GT. 1 )
1319 $ CALL zlaset('L',nr-1,nr-1,czero,czero,u(nr+2,1),ldu)
1320 CALL zgelqf( nr, n, u(nr+1,1), ldu, cwork(n+1),
1321 $ cwork(n+nr+1), lcwork-n-nr, ierr )
1322 CALL zlacpy('L',nr,nr,u(nr+1,1),ldu,v,ldv)
1323 IF ( nr .GT. 1 )
1324 $ CALL zlaset('U',nr-1,nr-1,czero,czero,v(1,2),ldv)
1325 CALL zgesvd( 'S', 'O', nr, nr, v, ldv, s, u, ldu,
1326 $ v, ldv, cwork(n+nr+1), lcwork-n-nr, rwork, info )
1327 CALL zlaset('A',n-nr,nr,czero,czero,v(nr+1,1),ldv)
1328 CALL zlaset('A',nr,n-nr,czero,czero,v(1,nr+1),ldv)
1329 CALL zlaset('A',n-nr,n-nr,czero,cone,v(nr+1,nr+1),ldv)
1330 CALL zunmlq('R','N',n,n,nr,u(nr+1,1),ldu,cwork(n+1),
1331 $ v, ldv, cwork(n+nr+1),lcwork-n-nr,ierr)
1332 CALL zlapmt( .false., n, n, v, ldv, iwork )
1333* .. assemble the left singular vector matrix U of dimensions
1334* (M x NR) or (M x N) or (M x M).
1335 IF ( ( nr .LT. m ) .AND. .NOT.(wntuf)) THEN
1336 CALL zlaset('A',m-nr,nr,czero,czero,u(nr+1,1),ldu)
1337 IF ( nr .LT. n1 ) THEN
1338 CALL zlaset('A',nr,n1-nr,czero,czero,u(1,nr+1),ldu)
1339 CALL zlaset( 'A',m-nr,n1-nr,czero,cone,
1340 $ u(nr+1,nr+1), ldu )
1341 END IF
1342 END IF
1343 END IF
1344 END IF
1345* .. end of the "R**H or R" branch
1346 END IF
1347*
1348* The Q matrix from the first QRF is built into the left singular
1349* vectors matrix U.
1350*
1351 IF ( .NOT. wntuf )
1352 $ CALL zunmqr( 'L', 'N', m, n1, n, a, lda, cwork, u,
1353 $ ldu, cwork(n+1), lcwork-n, ierr )
1354 IF ( rowprm .AND. .NOT.wntuf )
1355 $ CALL zlaswp( n1, u, ldu, 1, m-1, iwork(n+1), -1 )
1356*
1357* ... end of the "full SVD" branch
1358 END IF
1359*
1360* Check whether some singular values are returned as zeros, e.g.
1361* due to underflow, and update the numerical rank.
1362 p = nr
1363 DO 4001 q = p, 1, -1
1364 IF ( s(q) .GT. zero ) GO TO 4002
1365 nr = nr - 1
1366 4001 CONTINUE
1367 4002 CONTINUE
1368*
1369* .. if numerical rank deficiency is detected, the truncated
1370* singular values are set to zero.
1371 IF ( nr .LT. n ) CALL dlaset( 'G', n-nr,1, zero,zero, s(nr+1), n )
1372* .. undo scaling; this may cause overflow in the largest singular
1373* values.
1374 IF ( ascaled )
1375 $ CALL dlascl( 'G',0,0, one,sqrt(dble(m)), nr,1, s, n, ierr )
1376 IF ( conda ) rwork(1) = sconda
1377 rwork(2) = p - nr
1378* .. p-NR is the number of singular values that are computed as
1379* exact zeros in ZGESVD() applied to the (possibly truncated)
1380* full row rank triangular (trapezoidal) factor of A.
1381 numrank = nr
1382*
1383 RETURN
1384*
1385* End of ZGESVDQ
1386*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgelqf(m, n, a, lda, tau, work, lwork, info)
ZGELQF
Definition zgelqf.f:143
subroutine zgeqp3(m, n, a, lda, jpvt, tau, work, lwork, rwork, info)
ZGEQP3
Definition zgeqp3.f:159
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:146
subroutine zgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, rwork, info)
ZGESVD computes the singular value decomposition (SVD) for GE matrices
Definition zgesvd.f:214
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:115
subroutine zlapmt(forwrd, m, n, x, ldx, k)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition zlapmt.f:104
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:143
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zlaswp(n, a, lda, k1, k2, ipiv, incx)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition zlaswp.f:115
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
subroutine zpocon(uplo, n, a, lda, anorm, rcond, work, rwork, info)
ZPOCON
Definition zpocon.f:121
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine zunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMLQ
Definition zunmlq.f:167
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: