LAPACK 3.12.1
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 408 of file zgesvdq.f.

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