 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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
   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  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)
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: