LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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

Definition at line 410 of file zgesvdq.f.

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