 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

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]].
  Please report all bugs and send interesting examples and/or comments to
drmac@math.hr. Thank you.
References
   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

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  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: