LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

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

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

Definition at line 410 of file 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
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 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 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
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 zlaswp(N, A, LDA, K1, K2, IPIV, INCX)
ZLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: zlaswp.f:115
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
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 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 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
subroutine zpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
ZPOCON
Definition: zpocon.f:121
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition: dznrm2.f90:90
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
Here is the call graph for this function:
Here is the caller graph for this function: