LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ cgesvdq()

 subroutine cgesvdq ( character JOBA, character JOBP, character JOBR, character JOBU, character JOBV, integer M, integer N, complex, dimension( lda, * ) A, integer LDA, real, dimension( * ) S, complex, dimension( ldu, * ) U, integer LDU, complex, dimension( ldv, * ) V, integer LDV, integer NUMRANK, integer, dimension( * ) IWORK, integer LIWORK, complex, dimension( * ) CWORK, integer LCWORK, real, dimension( * ) RWORK, integer LRWORK, integer INFO )

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

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