1*> \brief <b> CGESVJ </b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CGESVJ + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvj.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvj.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvj.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, 22* LDV, CWORK, LWORK, RWORK, LRWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, LDA, LDV, LWORK, LRWORK, M, MV, N 26* CHARACTER*1 JOBA, JOBU, JOBV 27* .. 28* .. Array Arguments .. 29* COMPLEX A( LDA, * ), V( LDV, * ), CWORK( LWORK ) 30* REAL RWORK( LRWORK ), SVA( N ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> CGESVJ computes the singular value decomposition (SVD) of a complex 40*> M-by-N matrix A, where M >= N. The SVD of A is written as 41*> [++] [xx] [x0] [xx] 42*> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx] 43*> [++] [xx] 44*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal 45*> matrix, and V is an N-by-N unitary matrix. The diagonal elements 46*> of SIGMA are the singular values of A. The columns of U and V are the 47*> left and the right singular vectors of A, respectively. 48*> \endverbatim 49* 50* Arguments: 51* ========== 52* 53*> \param[in] JOBA 54*> \verbatim 55*> JOBA is CHARACTER*1 56*> Specifies the structure of A. 57*> = 'L': The input matrix A is lower triangular; 58*> = 'U': The input matrix A is upper triangular; 59*> = 'G': The input matrix A is general M-by-N matrix, M >= N. 60*> \endverbatim 61*> 62*> \param[in] JOBU 63*> \verbatim 64*> JOBU is CHARACTER*1 65*> Specifies whether to compute the left singular vectors 66*> (columns of U): 67*> = 'U' or 'F': The left singular vectors corresponding to the nonzero 68*> singular values are computed and returned in the leading 69*> columns of A. See more details in the description of A. 70*> The default numerical orthogonality threshold is set to 71*> approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E'). 72*> = 'C': Analogous to JOBU='U', except that user can control the 73*> level of numerical orthogonality of the computed left 74*> singular vectors. TOL can be set to TOL = CTOL*EPS, where 75*> CTOL is given on input in the array WORK. 76*> No CTOL smaller than ONE is allowed. CTOL greater 77*> than 1 / EPS is meaningless. The option 'C' 78*> can be used if M*EPS is satisfactory orthogonality 79*> of the computed left singular vectors, so CTOL=M could 80*> save few sweeps of Jacobi rotations. 81*> See the descriptions of A and WORK(1). 82*> = 'N': The matrix U is not computed. However, see the 83*> description of A. 84*> \endverbatim 85*> 86*> \param[in] JOBV 87*> \verbatim 88*> JOBV is CHARACTER*1 89*> Specifies whether to compute the right singular vectors, that 90*> is, the matrix V: 91*> = 'V' or 'J': the matrix V is computed and returned in the array V 92*> = 'A': the Jacobi rotations are applied to the MV-by-N 93*> array V. In other words, the right singular vector 94*> matrix V is not computed explicitly; instead it is 95*> applied to an MV-by-N matrix initially stored in the 96*> first MV rows of V. 97*> = 'N': the matrix V is not computed and the array V is not 98*> referenced 99*> \endverbatim 100*> 101*> \param[in] M 102*> \verbatim 103*> M is INTEGER 104*> The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0. 105*> \endverbatim 106*> 107*> \param[in] N 108*> \verbatim 109*> N is INTEGER 110*> The number of columns of the input matrix A. 111*> M >= N >= 0. 112*> \endverbatim 113*> 114*> \param[in,out] A 115*> \verbatim 116*> A is COMPLEX array, dimension (LDA,N) 117*> On entry, the M-by-N matrix A. 118*> On exit, 119*> If JOBU = 'U' .OR. JOBU = 'C': 120*> If INFO = 0 : 121*> RANKA orthonormal columns of U are returned in the 122*> leading RANKA columns of the array A. Here RANKA <= N 123*> is the number of computed singular values of A that are 124*> above the underflow threshold SLAMCH('S'). The singular 125*> vectors corresponding to underflowed or zero singular 126*> values are not computed. The value of RANKA is returned 127*> in the array RWORK as RANKA=NINT(RWORK(2)). Also see the 128*> descriptions of SVA and RWORK. The computed columns of U 129*> are mutually numerically orthogonal up to approximately 130*> TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'), 131*> see the description of JOBU. 132*> If INFO > 0, 133*> the procedure CGESVJ did not converge in the given number 134*> of iterations (sweeps). In that case, the computed 135*> columns of U may not be orthogonal up to TOL. The output 136*> U (stored in A), SIGMA (given by the computed singular 137*> values in SVA(1:N)) and V is still a decomposition of the 138*> input matrix A in the sense that the residual 139*> || A - SCALE * U * SIGMA * V^* ||_2 / ||A||_2 is small. 140*> If JOBU = 'N': 141*> If INFO = 0 : 142*> Note that the left singular vectors are 'for free' in the 143*> one-sided Jacobi SVD algorithm. However, if only the 144*> singular values are needed, the level of numerical 145*> orthogonality of U is not an issue and iterations are 146*> stopped when the columns of the iterated matrix are 147*> numerically orthogonal up to approximately M*EPS. Thus, 148*> on exit, A contains the columns of U scaled with the 149*> corresponding singular values. 150*> If INFO > 0 : 151*> the procedure CGESVJ did not converge in the given number 152*> of iterations (sweeps). 153*> \endverbatim 154*> 155*> \param[in] LDA 156*> \verbatim 157*> LDA is INTEGER 158*> The leading dimension of the array A. LDA >= max(1,M). 159*> \endverbatim 160*> 161*> \param[out] SVA 162*> \verbatim 163*> SVA is REAL array, dimension (N) 164*> On exit, 165*> If INFO = 0 : 166*> depending on the value SCALE = RWORK(1), we have: 167*> If SCALE = ONE: 168*> SVA(1:N) contains the computed singular values of A. 169*> During the computation SVA contains the Euclidean column 170*> norms of the iterated matrices in the array A. 171*> If SCALE .NE. ONE: 172*> The singular values of A are SCALE*SVA(1:N), and this 173*> factored representation is due to the fact that some of the 174*> singular values of A might underflow or overflow. 175*> 176*> If INFO > 0 : 177*> the procedure CGESVJ did not converge in the given number of 178*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. 179*> \endverbatim 180*> 181*> \param[in] MV 182*> \verbatim 183*> MV is INTEGER 184*> If JOBV = 'A', then the product of Jacobi rotations in CGESVJ 185*> is applied to the first MV rows of V. See the description of JOBV. 186*> \endverbatim 187*> 188*> \param[in,out] V 189*> \verbatim 190*> V is COMPLEX array, dimension (LDV,N) 191*> If JOBV = 'V', then V contains on exit the N-by-N matrix of 192*> the right singular vectors; 193*> If JOBV = 'A', then V contains the product of the computed right 194*> singular vector matrix and the initial matrix in 195*> the array V. 196*> If JOBV = 'N', then V is not referenced. 197*> \endverbatim 198*> 199*> \param[in] LDV 200*> \verbatim 201*> LDV is INTEGER 202*> The leading dimension of the array V, LDV >= 1. 203*> If JOBV = 'V', then LDV >= max(1,N). 204*> If JOBV = 'A', then LDV >= max(1,MV) . 205*> \endverbatim 206*> 207*> \param[in,out] CWORK 208*> \verbatim 209*> CWORK is COMPLEX array, dimension (max(1,LWORK)) 210*> Used as workspace. 211*> If on entry LWORK = -1, then a workspace query is assumed and 212*> no computation is done; CWORK(1) is set to the minial (and optimal) 213*> length of CWORK. 214*> \endverbatim 215*> 216*> \param[in] LWORK 217*> \verbatim 218*> LWORK is INTEGER. 219*> Length of CWORK, LWORK >= M+N. 220*> \endverbatim 221*> 222*> \param[in,out] RWORK 223*> \verbatim 224*> RWORK is REAL array, dimension (max(6,LRWORK)) 225*> On entry, 226*> If JOBU = 'C' : 227*> RWORK(1) = CTOL, where CTOL defines the threshold for convergence. 228*> The process stops if all columns of A are mutually 229*> orthogonal up to CTOL*EPS, EPS=SLAMCH('E'). 230*> It is required that CTOL >= ONE, i.e. it is not 231*> allowed to force the routine to obtain orthogonality 232*> below EPSILON. 233*> On exit, 234*> RWORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) 235*> are the computed singular values of A. 236*> (See description of SVA().) 237*> RWORK(2) = NINT(RWORK(2)) is the number of the computed nonzero 238*> singular values. 239*> RWORK(3) = NINT(RWORK(3)) is the number of the computed singular 240*> values that are larger than the underflow threshold. 241*> RWORK(4) = NINT(RWORK(4)) is the number of sweeps of Jacobi 242*> rotations needed for numerical convergence. 243*> RWORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. 244*> This is useful information in cases when CGESVJ did 245*> not converge, as it can be used to estimate whether 246*> the output is still useful and for post festum analysis. 247*> RWORK(6) = the largest absolute value over all sines of the 248*> Jacobi rotation angles in the last sweep. It can be 249*> useful for a post festum analysis. 250*> If on entry LRWORK = -1, then a workspace query is assumed and 251*> no computation is done; RWORK(1) is set to the minial (and optimal) 252*> length of RWORK. 253*> \endverbatim 254*> 255*> \param[in] LRWORK 256*> \verbatim 257*> LRWORK is INTEGER 258*> Length of RWORK, LRWORK >= MAX(6,N). 259*> \endverbatim 260*> 261*> \param[out] INFO 262*> \verbatim 263*> INFO is INTEGER 264*> = 0: successful exit. 265*> < 0: if INFO = -i, then the i-th argument had an illegal value 266*> > 0: CGESVJ did not converge in the maximal allowed number 267*> (NSWEEP=30) of sweeps. The output may still be useful. 268*> See the description of RWORK. 269*> \endverbatim 270*> 271* Authors: 272* ======== 273* 274*> \author Univ. of Tennessee 275*> \author Univ. of California Berkeley 276*> \author Univ. of Colorado Denver 277*> \author NAG Ltd. 278* 279*> \ingroup complexGEcomputational 280* 281*> \par Further Details: 282* ===================== 283*> 284*> \verbatim 285*> 286*> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane 287*> rotations. In the case of underflow of the tangent of the Jacobi angle, a 288*> modified Jacobi transformation of Drmac [3] is used. Pivot strategy uses 289*> column interchanges of de Rijk [1]. The relative accuracy of the computed 290*> singular values and the accuracy of the computed singular vectors (in 291*> angle metric) is as guaranteed by the theory of Demmel and Veselic [2]. 292*> The condition number that determines the accuracy in the full rank case 293*> is essentially min_{D=diag} kappa(A*D), where kappa(.) is the 294*> spectral condition number. The best performance of this Jacobi SVD 295*> procedure is achieved if used in an accelerated version of Drmac and 296*> Veselic [4,5], and it is the kernel routine in the SIGMA library [6]. 297*> Some tuning parameters (marked with [TP]) are available for the 298*> implementer. 299*> The computational range for the nonzero singular values is the machine 300*> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even 301*> denormalized singular values can be computed with the corresponding 302*> gradual loss of accurate digits. 303*> \endverbatim 304* 305*> \par Contributor: 306* ================== 307*> 308*> \verbatim 309*> 310*> ============ 311*> 312*> Zlatko Drmac (Zagreb, Croatia) 313*> 314*> \endverbatim 315* 316*> \par References: 317* ================ 318*> 319*> \verbatim 320*> 321*> [1] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the 322*> singular value decomposition on a vector computer. 323*> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. 324*> [2] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. 325*> [3] Z. Drmac: Implementation of Jacobi rotations for accurate singular 326*> value computation in floating point arithmetic. 327*> SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. 328*> [4] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. 329*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. 330*> LAPACK Working note 169. 331*> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. 332*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. 333*> LAPACK Working note 170. 334*> [6] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, 335*> QSVD, (H,K)-SVD computations. 336*> Department of Mathematics, University of Zagreb, 2008, 2015. 337*> \endverbatim 338* 339*> \par Bugs, examples and comments: 340* ================================= 341*> 342*> \verbatim 343*> =========================== 344*> Please report all bugs and send interesting test examples and comments to 345*> drmac@math.hr. Thank you. 346*> \endverbatim 347*> 348* ===================================================================== 349 SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, 350 $ LDV, CWORK, LWORK, RWORK, LRWORK, INFO ) 351* 352* -- LAPACK computational routine -- 353* -- LAPACK is a software package provided by Univ. of Tennessee, -- 354* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 355* 356 IMPLICIT NONE 357* .. Scalar Arguments .. 358 INTEGER INFO, LDA, LDV, LWORK, LRWORK, M, MV, N 359 CHARACTER*1 JOBA, JOBU, JOBV 360* .. 361* .. Array Arguments .. 362 COMPLEX A( LDA, * ), V( LDV, * ), CWORK( LWORK ) 363 REAL RWORK( LRWORK ), SVA( N ) 364* .. 365* 366* ===================================================================== 367* 368* .. Local Parameters .. 369 REAL ZERO, HALF, ONE 370 PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0) 371 COMPLEX CZERO, CONE 372 PARAMETER ( CZERO = (0.0E0, 0.0E0), CONE = (1.0E0, 0.0E0) ) 373 INTEGER NSWEEP 374 PARAMETER ( NSWEEP = 30 ) 375* .. 376* .. Local Scalars .. 377 COMPLEX AAPQ, OMPQ 378 REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG, 379 $ BIGTHETA, CS, CTOL, EPSLN, MXAAPQ, 380 $ MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL, 381 $ SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL 382 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, 383 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34, 384 $ N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND 385 LOGICAL APPLV, GOSCALE, LOWER, LQUERY, LSVEC, NOSCALE, ROTOK, 386 $ RSVEC, UCTOL, UPPER 387* .. 388* .. 389* .. Intrinsic Functions .. 390 INTRINSIC ABS, MAX, MIN, CONJG, REAL, SIGN, SQRT 391* .. 392* .. External Functions .. 393* .. 394* from BLAS 395 REAL SCNRM2 396 COMPLEX CDOTC 397 EXTERNAL CDOTC, SCNRM2 398 INTEGER ISAMAX 399 EXTERNAL ISAMAX 400* from LAPACK 401 REAL SLAMCH 402 EXTERNAL SLAMCH 403 LOGICAL LSAME 404 EXTERNAL LSAME 405* .. 406* .. External Subroutines .. 407* .. 408* from BLAS 409 EXTERNAL CCOPY, CROT, CSSCAL, CSWAP, CAXPY 410* from LAPACK 411 EXTERNAL CLASCL, CLASET, CLASSQ, SLASCL, XERBLA 412 EXTERNAL CGSVJ0, CGSVJ1 413* .. 414* .. Executable Statements .. 415* 416* Test the input arguments 417* 418 LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' ) 419 UCTOL = LSAME( JOBU, 'C' ) 420 RSVEC = LSAME( JOBV, 'V' ) .OR. LSAME( JOBV, 'J' ) 421 APPLV = LSAME( JOBV, 'A' ) 422 UPPER = LSAME( JOBA, 'U' ) 423 LOWER = LSAME( JOBA, 'L' ) 424* 425 LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 ) 426 IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN 427 INFO = -1 428 ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN 429 INFO = -2 430 ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN 431 INFO = -3 432 ELSE IF( M.LT.0 ) THEN 433 INFO = -4 434 ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN 435 INFO = -5 436 ELSE IF( LDA.LT.M ) THEN 437 INFO = -7 438 ELSE IF( MV.LT.0 ) THEN 439 INFO = -9 440 ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR. 441 $ ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN 442 INFO = -11 443 ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN 444 INFO = -12 445 ELSE IF( LWORK.LT.( M+N ) .AND. ( .NOT.LQUERY ) ) THEN 446 INFO = -13 447 ELSE IF( LRWORK.LT.MAX( N, 6 ) .AND. ( .NOT.LQUERY ) ) THEN 448 INFO = -15 449 ELSE 450 INFO = 0 451 END IF 452* 453* #:( 454 IF( INFO.NE.0 ) THEN 455 CALL XERBLA( 'CGESVJ', -INFO ) 456 RETURN 457 ELSE IF ( LQUERY ) THEN 458 CWORK(1) = M + N 459 RWORK(1) = MAX( N, 6 ) 460 RETURN 461 END IF 462* 463* #:) Quick return for void matrix 464* 465 IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN 466* 467* Set numerical parameters 468* The stopping criterion for Jacobi rotations is 469* 470* max_{i<>j}|A(:,i)^* * A(:,j)| / (||A(:,i)||*||A(:,j)||) < CTOL*EPS 471* 472* where EPS is the round-off and CTOL is defined as follows: 473* 474 IF( UCTOL ) THEN 475* ... user controlled 476 CTOL = RWORK( 1 ) 477 ELSE 478* ... default 479 IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN 480 CTOL = SQRT( REAL( M ) ) 481 ELSE 482 CTOL = REAL( M ) 483 END IF 484 END IF 485* ... and the machine dependent parameters are 486*[!] (Make sure that SLAMCH() works properly on the target machine.) 487* 488 EPSLN = SLAMCH( 'Epsilon' ) 489 ROOTEPS = SQRT( EPSLN ) 490 SFMIN = SLAMCH( 'SafeMinimum' ) 491 ROOTSFMIN = SQRT( SFMIN ) 492 SMALL = SFMIN / EPSLN 493* BIG = SLAMCH( 'Overflow' ) 494 BIG = ONE / SFMIN 495 ROOTBIG = ONE / ROOTSFMIN 496* LARGE = BIG / SQRT( REAL( M*N ) ) 497 BIGTHETA = ONE / ROOTEPS 498* 499 TOL = CTOL*EPSLN 500 ROOTTOL = SQRT( TOL ) 501* 502 IF( REAL( M )*EPSLN.GE.ONE ) THEN 503 INFO = -4 504 CALL XERBLA( 'CGESVJ', -INFO ) 505 RETURN 506 END IF 507* 508* Initialize the right singular vector matrix. 509* 510 IF( RSVEC ) THEN 511 MVL = N 512 CALL CLASET( 'A', MVL, N, CZERO, CONE, V, LDV ) 513 ELSE IF( APPLV ) THEN 514 MVL = MV 515 END IF 516 RSVEC = RSVEC .OR. APPLV 517* 518* Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N ) 519*(!) If necessary, scale A to protect the largest singular value 520* from overflow. It is possible that saving the largest singular 521* value destroys the information about the small ones. 522* This initial scaling is almost minimal in the sense that the 523* goal is to make sure that no column norm overflows, and that 524* SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries 525* in A are detected, the procedure returns with INFO=-6. 526* 527 SKL = ONE / SQRT( REAL( M )*REAL( N ) ) 528 NOSCALE = .TRUE. 529 GOSCALE = .TRUE. 530* 531 IF( LOWER ) THEN 532* the input matrix is M-by-N lower triangular (trapezoidal) 533 DO 1874 p = 1, N 534 AAPP = ZERO 535 AAQQ = ONE 536 CALL CLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ ) 537 IF( AAPP.GT.BIG ) THEN 538 INFO = -6 539 CALL XERBLA( 'CGESVJ', -INFO ) 540 RETURN 541 END IF 542 AAQQ = SQRT( AAQQ ) 543 IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN 544 SVA( p ) = AAPP*AAQQ 545 ELSE 546 NOSCALE = .FALSE. 547 SVA( p ) = AAPP*( AAQQ*SKL ) 548 IF( GOSCALE ) THEN 549 GOSCALE = .FALSE. 550 DO 1873 q = 1, p - 1 551 SVA( q ) = SVA( q )*SKL 552 1873 CONTINUE 553 END IF 554 END IF 555 1874 CONTINUE 556 ELSE IF( UPPER ) THEN 557* the input matrix is M-by-N upper triangular (trapezoidal) 558 DO 2874 p = 1, N 559 AAPP = ZERO 560 AAQQ = ONE 561 CALL CLASSQ( p, A( 1, p ), 1, AAPP, AAQQ ) 562 IF( AAPP.GT.BIG ) THEN 563 INFO = -6 564 CALL XERBLA( 'CGESVJ', -INFO ) 565 RETURN 566 END IF 567 AAQQ = SQRT( AAQQ ) 568 IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN 569 SVA( p ) = AAPP*AAQQ 570 ELSE 571 NOSCALE = .FALSE. 572 SVA( p ) = AAPP*( AAQQ*SKL ) 573 IF( GOSCALE ) THEN 574 GOSCALE = .FALSE. 575 DO 2873 q = 1, p - 1 576 SVA( q ) = SVA( q )*SKL 577 2873 CONTINUE 578 END IF 579 END IF 580 2874 CONTINUE 581 ELSE 582* the input matrix is M-by-N general dense 583 DO 3874 p = 1, N 584 AAPP = ZERO 585 AAQQ = ONE 586 CALL CLASSQ( M, A( 1, p ), 1, AAPP, AAQQ ) 587 IF( AAPP.GT.BIG ) THEN 588 INFO = -6 589 CALL XERBLA( 'CGESVJ', -INFO ) 590 RETURN 591 END IF 592 AAQQ = SQRT( AAQQ ) 593 IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN 594 SVA( p ) = AAPP*AAQQ 595 ELSE 596 NOSCALE = .FALSE. 597 SVA( p ) = AAPP*( AAQQ*SKL ) 598 IF( GOSCALE ) THEN 599 GOSCALE = .FALSE. 600 DO 3873 q = 1, p - 1 601 SVA( q ) = SVA( q )*SKL 602 3873 CONTINUE 603 END IF 604 END IF 605 3874 CONTINUE 606 END IF 607* 608 IF( NOSCALE )SKL = ONE 609* 610* Move the smaller part of the spectrum from the underflow threshold 611*(!) Start by determining the position of the nonzero entries of the 612* array SVA() relative to ( SFMIN, BIG ). 613* 614 AAPP = ZERO 615 AAQQ = BIG 616 DO 4781 p = 1, N 617 IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) ) 618 AAPP = MAX( AAPP, SVA( p ) ) 619 4781 CONTINUE 620* 621* #:) Quick return for zero matrix 622* 623 IF( AAPP.EQ.ZERO ) THEN 624 IF( LSVEC )CALL CLASET( 'G', M, N, CZERO, CONE, A, LDA ) 625 RWORK( 1 ) = ONE 626 RWORK( 2 ) = ZERO 627 RWORK( 3 ) = ZERO 628 RWORK( 4 ) = ZERO 629 RWORK( 5 ) = ZERO 630 RWORK( 6 ) = ZERO 631 RETURN 632 END IF 633* 634* #:) Quick return for one-column matrix 635* 636 IF( N.EQ.1 ) THEN 637 IF( LSVEC )CALL CLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1, 638 $ A( 1, 1 ), LDA, IERR ) 639 RWORK( 1 ) = ONE / SKL 640 IF( SVA( 1 ).GE.SFMIN ) THEN 641 RWORK( 2 ) = ONE 642 ELSE 643 RWORK( 2 ) = ZERO 644 END IF 645 RWORK( 3 ) = ZERO 646 RWORK( 4 ) = ZERO 647 RWORK( 5 ) = ZERO 648 RWORK( 6 ) = ZERO 649 RETURN 650 END IF 651* 652* Protect small singular values from underflow, and try to 653* avoid underflows/overflows in computing Jacobi rotations. 654* 655 SN = SQRT( SFMIN / EPSLN ) 656 TEMP1 = SQRT( BIG / REAL( N ) ) 657 IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR. 658 $ ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN 659 TEMP1 = MIN( BIG, TEMP1 / AAPP ) 660* AAQQ = AAQQ*TEMP1 661* AAPP = AAPP*TEMP1 662 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN 663 TEMP1 = MIN( SN / AAQQ, BIG / ( AAPP*SQRT( REAL( N ) ) ) ) 664* AAQQ = AAQQ*TEMP1 665* AAPP = AAPP*TEMP1 666 ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN 667 TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP ) 668* AAQQ = AAQQ*TEMP1 669* AAPP = AAPP*TEMP1 670 ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN 671 TEMP1 = MIN( SN / AAQQ, BIG / ( SQRT( REAL( N ) )*AAPP ) ) 672* AAQQ = AAQQ*TEMP1 673* AAPP = AAPP*TEMP1 674 ELSE 675 TEMP1 = ONE 676 END IF 677* 678* Scale, if necessary 679* 680 IF( TEMP1.NE.ONE ) THEN 681 CALL SLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR ) 682 END IF 683 SKL = TEMP1*SKL 684 IF( SKL.NE.ONE ) THEN 685 CALL CLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR ) 686 SKL = ONE / SKL 687 END IF 688* 689* Row-cyclic Jacobi SVD algorithm with column pivoting 690* 691 EMPTSW = ( N*( N-1 ) ) / 2 692 NOTROT = 0 693 694 DO 1868 q = 1, N 695 CWORK( q ) = CONE 696 1868 CONTINUE 697* 698* 699* 700 SWBAND = 3 701*[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective 702* if CGESVJ is used as a computational routine in the preconditioned 703* Jacobi SVD algorithm CGEJSV. For sweeps i=1:SWBAND the procedure 704* works on pivots inside a band-like region around the diagonal. 705* The boundaries are determined dynamically, based on the number of 706* pivots above a threshold. 707* 708 KBL = MIN( 8, N ) 709*[TP] KBL is a tuning parameter that defines the tile size in the 710* tiling of the p-q loops of pivot pairs. In general, an optimal 711* value of KBL depends on the matrix dimensions and on the 712* parameters of the computer's memory. 713* 714 NBL = N / KBL 715 IF( ( NBL*KBL ).NE.N )NBL = NBL + 1 716* 717 BLSKIP = KBL**2 718*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. 719* 720 ROWSKIP = MIN( 5, KBL ) 721*[TP] ROWSKIP is a tuning parameter. 722* 723 LKAHEAD = 1 724*[TP] LKAHEAD is a tuning parameter. 725* 726* Quasi block transformations, using the lower (upper) triangular 727* structure of the input matrix. The quasi-block-cycling usually 728* invokes cubic convergence. Big part of this cycle is done inside 729* canonical subspaces of dimensions less than M. 730* 731 IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN 732*[TP] The number of partition levels and the actual partition are 733* tuning parameters. 734 N4 = N / 4 735 N2 = N / 2 736 N34 = 3*N4 737 IF( APPLV ) THEN 738 q = 0 739 ELSE 740 q = 1 741 END IF 742* 743 IF( LOWER ) THEN 744* 745* This works very well on lower triangular matrices, in particular 746* in the framework of the preconditioned Jacobi SVD (xGEJSV). 747* The idea is simple: 748* [+ 0 0 0] Note that Jacobi transformations of [0 0] 749* [+ + 0 0] [0 0] 750* [+ + x 0] actually work on [x 0] [x 0] 751* [+ + x x] [x x]. [x x] 752* 753 CALL CGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA, 754 $ CWORK( N34+1 ), SVA( N34+1 ), MVL, 755 $ V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL, 756 $ 2, CWORK( N+1 ), LWORK-N, IERR ) 757 758 CALL CGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA, 759 $ CWORK( N2+1 ), SVA( N2+1 ), MVL, 760 $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2, 761 $ CWORK( N+1 ), LWORK-N, IERR ) 762 763 CALL CGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA, 764 $ CWORK( N2+1 ), SVA( N2+1 ), MVL, 765 $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, 766 $ CWORK( N+1 ), LWORK-N, IERR ) 767* 768 CALL CGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA, 769 $ CWORK( N4+1 ), SVA( N4+1 ), MVL, 770 $ V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1, 771 $ CWORK( N+1 ), LWORK-N, IERR ) 772* 773 CALL CGSVJ0( JOBV, M, N4, A, LDA, CWORK, SVA, MVL, V, LDV, 774 $ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N, 775 $ IERR ) 776* 777 CALL CGSVJ1( JOBV, M, N2, N4, A, LDA, CWORK, SVA, MVL, V, 778 $ LDV, EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), 779 $ LWORK-N, IERR ) 780* 781* 782 ELSE IF( UPPER ) THEN 783* 784* 785 CALL CGSVJ0( JOBV, N4, N4, A, LDA, CWORK, SVA, MVL, V, LDV, 786 $ EPSLN, SFMIN, TOL, 2, CWORK( N+1 ), LWORK-N, 787 $ IERR ) 788* 789 CALL CGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, CWORK( N4+1 ), 790 $ SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV, 791 $ EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), LWORK-N, 792 $ IERR ) 793* 794 CALL CGSVJ1( JOBV, N2, N2, N4, A, LDA, CWORK, SVA, MVL, V, 795 $ LDV, EPSLN, SFMIN, TOL, 1, CWORK( N+1 ), 796 $ LWORK-N, IERR ) 797* 798 CALL CGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA, 799 $ CWORK( N2+1 ), SVA( N2+1 ), MVL, 800 $ V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1, 801 $ CWORK( N+1 ), LWORK-N, IERR ) 802 803 END IF 804* 805 END IF 806* 807* .. Row-cyclic pivot strategy with de Rijk's pivoting .. 808* 809 DO 1993 i = 1, NSWEEP 810* 811* .. go go go ... 812* 813 MXAAPQ = ZERO 814 MXSINJ = ZERO 815 ISWROT = 0 816* 817 NOTROT = 0 818 PSKIPPED = 0 819* 820* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs 821* 1 <= p < q <= N. This is the first step toward a blocked implementation 822* of the rotations. New implementation, based on block transformations, 823* is under development. 824* 825 DO 2000 ibr = 1, NBL 826* 827 igl = ( ibr-1 )*KBL + 1 828* 829 DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr ) 830* 831 igl = igl + ir1*KBL 832* 833 DO 2001 p = igl, MIN( igl+KBL-1, N-1 ) 834* 835* .. de Rijk's pivoting 836* 837 q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1 838 IF( p.NE.q ) THEN 839 CALL CSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 840 IF( RSVEC )CALL CSWAP( MVL, V( 1, p ), 1, 841 $ V( 1, q ), 1 ) 842 TEMP1 = SVA( p ) 843 SVA( p ) = SVA( q ) 844 SVA( q ) = TEMP1 845 AAPQ = CWORK(p) 846 CWORK(p) = CWORK(q) 847 CWORK(q) = AAPQ 848 END IF 849* 850 IF( ir1.EQ.0 ) THEN 851* 852* Column norms are periodically updated by explicit 853* norm computation. 854*[!] Caveat: 855* Unfortunately, some BLAS implementations compute SCNRM2(M,A(1,p),1) 856* as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to 857* overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to 858* underflow for ||A(:,p)||_2 < SQRT(underflow_threshold). 859* Hence, SCNRM2 cannot be trusted, not even in the case when 860* the true norm is far from the under(over)flow boundaries. 861* If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF 862* below should be replaced with "AAPP = SCNRM2( M, A(1,p), 1 )". 863* 864 IF( ( SVA( p ).LT.ROOTBIG ) .AND. 865 $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN 866 SVA( p ) = SCNRM2( M, A( 1, p ), 1 ) 867 ELSE 868 TEMP1 = ZERO 869 AAPP = ONE 870 CALL CLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) 871 SVA( p ) = TEMP1*SQRT( AAPP ) 872 END IF 873 AAPP = SVA( p ) 874 ELSE 875 AAPP = SVA( p ) 876 END IF 877* 878 IF( AAPP.GT.ZERO ) THEN 879* 880 PSKIPPED = 0 881* 882 DO 2002 q = p + 1, MIN( igl+KBL-1, N ) 883* 884 AAQQ = SVA( q ) 885* 886 IF( AAQQ.GT.ZERO ) THEN 887* 888 AAPP0 = AAPP 889 IF( AAQQ.GE.ONE ) THEN 890 ROTOK = ( SMALL*AAPP ).LE.AAQQ 891 IF( AAPP.LT.( BIG / AAQQ ) ) THEN 892 AAPQ = ( CDOTC( M, A( 1, p ), 1, 893 $ A( 1, q ), 1 ) / AAQQ ) / AAPP 894 ELSE 895 CALL CCOPY( M, A( 1, p ), 1, 896 $ CWORK(N+1), 1 ) 897 CALL CLASCL( 'G', 0, 0, AAPP, ONE, 898 $ M, 1, CWORK(N+1), LDA, IERR ) 899 AAPQ = CDOTC( M, CWORK(N+1), 1, 900 $ A( 1, q ), 1 ) / AAQQ 901 END IF 902 ELSE 903 ROTOK = AAPP.LE.( AAQQ / SMALL ) 904 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 905 AAPQ = ( CDOTC( M, A( 1, p ), 1, 906 $ A( 1, q ), 1 ) / AAPP ) / AAQQ 907 ELSE 908 CALL CCOPY( M, A( 1, q ), 1, 909 $ CWORK(N+1), 1 ) 910 CALL CLASCL( 'G', 0, 0, AAQQ, 911 $ ONE, M, 1, 912 $ CWORK(N+1), LDA, IERR ) 913 AAPQ = CDOTC( M, A(1, p ), 1, 914 $ CWORK(N+1), 1 ) / AAPP 915 END IF 916 END IF 917* 918* AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q) 919 AAPQ1 = -ABS(AAPQ) 920 MXAAPQ = MAX( MXAAPQ, -AAPQ1 ) 921* 922* TO rotate or NOT to rotate, THAT is the question ... 923* 924 IF( ABS( AAPQ1 ).GT.TOL ) THEN 925 OMPQ = AAPQ / ABS(AAPQ) 926* 927* .. rotate 928*[RTD] ROTATED = ROTATED + ONE 929* 930 IF( ir1.EQ.0 ) THEN 931 NOTROT = 0 932 PSKIPPED = 0 933 ISWROT = ISWROT + 1 934 END IF 935* 936 IF( ROTOK ) THEN 937* 938 AQOAP = AAQQ / AAPP 939 APOAQ = AAPP / AAQQ 940 THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1 941* 942 IF( ABS( THETA ).GT.BIGTHETA ) THEN 943* 944 T = HALF / THETA 945 CS = ONE 946 947 CALL CROT( M, A(1,p), 1, A(1,q), 1, 948 $ CS, CONJG(OMPQ)*T ) 949 IF ( RSVEC ) THEN 950 CALL CROT( MVL, V(1,p), 1, 951 $ V(1,q), 1, CS, CONJG(OMPQ)*T ) 952 END IF 953 954 SVA( q ) = AAQQ*SQRT( MAX( ZERO, 955 $ ONE+T*APOAQ*AAPQ1 ) ) 956 AAPP = AAPP*SQRT( MAX( ZERO, 957 $ ONE-T*AQOAP*AAPQ1 ) ) 958 MXSINJ = MAX( MXSINJ, ABS( T ) ) 959* 960 ELSE 961* 962* .. choose correct signum for THETA and rotate 963* 964 THSIGN = -SIGN( ONE, AAPQ1 ) 965 T = ONE / ( THETA+THSIGN* 966 $ SQRT( ONE+THETA*THETA ) ) 967 CS = SQRT( ONE / ( ONE+T*T ) ) 968 SN = T*CS 969* 970 MXSINJ = MAX( MXSINJ, ABS( SN ) ) 971 SVA( q ) = AAQQ*SQRT( MAX( ZERO, 972 $ ONE+T*APOAQ*AAPQ1 ) ) 973 AAPP = AAPP*SQRT( MAX( ZERO, 974 $ ONE-T*AQOAP*AAPQ1 ) ) 975* 976 CALL CROT( M, A(1,p), 1, A(1,q), 1, 977 $ CS, CONJG(OMPQ)*SN ) 978 IF ( RSVEC ) THEN 979 CALL CROT( MVL, V(1,p), 1, 980 $ V(1,q), 1, CS, CONJG(OMPQ)*SN ) 981 END IF 982 END IF 983 CWORK(p) = -CWORK(q) * OMPQ 984* 985 ELSE 986* .. have to use modified Gram-Schmidt like transformation 987 CALL CCOPY( M, A( 1, p ), 1, 988 $ CWORK(N+1), 1 ) 989 CALL CLASCL( 'G', 0, 0, AAPP, ONE, M, 990 $ 1, CWORK(N+1), LDA, 991 $ IERR ) 992 CALL CLASCL( 'G', 0, 0, AAQQ, ONE, M, 993 $ 1, A( 1, q ), LDA, IERR ) 994 CALL CAXPY( M, -AAPQ, CWORK(N+1), 1, 995 $ A( 1, q ), 1 ) 996 CALL CLASCL( 'G', 0, 0, ONE, AAQQ, M, 997 $ 1, A( 1, q ), LDA, IERR ) 998 SVA( q ) = AAQQ*SQRT( MAX( ZERO, 999 $ ONE-AAPQ1*AAPQ1 ) ) 1000 MXSINJ = MAX( MXSINJ, SFMIN ) 1001 END IF 1002* END IF ROTOK THEN ... ELSE 1003* 1004* In the case of cancellation in updating SVA(q), SVA(p) 1005* recompute SVA(q), SVA(p). 1006* 1007 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 1008 $ THEN 1009 IF( ( AAQQ.LT.ROOTBIG ) .AND. 1010 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN 1011 SVA( q ) = SCNRM2( M, A( 1, q ), 1 ) 1012 ELSE 1013 T = ZERO 1014 AAQQ = ONE 1015 CALL CLASSQ( M, A( 1, q ), 1, T, 1016 $ AAQQ ) 1017 SVA( q ) = T*SQRT( AAQQ ) 1018 END IF 1019 END IF 1020 IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN 1021 IF( ( AAPP.LT.ROOTBIG ) .AND. 1022 $ ( AAPP.GT.ROOTSFMIN ) ) THEN 1023 AAPP = SCNRM2( M, A( 1, p ), 1 ) 1024 ELSE 1025 T = ZERO 1026 AAPP = ONE 1027 CALL CLASSQ( M, A( 1, p ), 1, T, 1028 $ AAPP ) 1029 AAPP = T*SQRT( AAPP ) 1030 END IF 1031 SVA( p ) = AAPP 1032 END IF 1033* 1034 ELSE 1035* A(:,p) and A(:,q) already numerically orthogonal 1036 IF( ir1.EQ.0 )NOTROT = NOTROT + 1 1037*[RTD] SKIPPED = SKIPPED + 1 1038 PSKIPPED = PSKIPPED + 1 1039 END IF 1040 ELSE 1041* A(:,q) is zero column 1042 IF( ir1.EQ.0 )NOTROT = NOTROT + 1 1043 PSKIPPED = PSKIPPED + 1 1044 END IF 1045* 1046 IF( ( i.LE.SWBAND ) .AND. 1047 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN 1048 IF( ir1.EQ.0 )AAPP = -AAPP 1049 NOTROT = 0 1050 GO TO 2103 1051 END IF 1052* 1053 2002 CONTINUE 1054* END q-LOOP 1055* 1056 2103 CONTINUE 1057* bailed out of q-loop 1058* 1059 SVA( p ) = AAPP 1060* 1061 ELSE 1062 SVA( p ) = AAPP 1063 IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) 1064 $ NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p 1065 END IF 1066* 1067 2001 CONTINUE 1068* end of the p-loop 1069* end of doing the block ( ibr, ibr ) 1070 1002 CONTINUE 1071* end of ir1-loop 1072* 1073* ... go to the off diagonal blocks 1074* 1075 igl = ( ibr-1 )*KBL + 1 1076* 1077 DO 2010 jbc = ibr + 1, NBL 1078* 1079 jgl = ( jbc-1 )*KBL + 1 1080* 1081* doing the block at ( ibr, jbc ) 1082* 1083 IJBLSK = 0 1084 DO 2100 p = igl, MIN( igl+KBL-1, N ) 1085* 1086 AAPP = SVA( p ) 1087 IF( AAPP.GT.ZERO ) THEN 1088* 1089 PSKIPPED = 0 1090* 1091 DO 2200 q = jgl, MIN( jgl+KBL-1, N ) 1092* 1093 AAQQ = SVA( q ) 1094 IF( AAQQ.GT.ZERO ) THEN 1095 AAPP0 = AAPP 1096* 1097* .. M x 2 Jacobi SVD .. 1098* 1099* Safe Gram matrix computation 1100* 1101 IF( AAQQ.GE.ONE ) THEN 1102 IF( AAPP.GE.AAQQ ) THEN 1103 ROTOK = ( SMALL*AAPP ).LE.AAQQ 1104 ELSE 1105 ROTOK = ( SMALL*AAQQ ).LE.AAPP 1106 END IF 1107 IF( AAPP.LT.( BIG / AAQQ ) ) THEN 1108 AAPQ = ( CDOTC( M, A( 1, p ), 1, 1109 $ A( 1, q ), 1 ) / AAQQ ) / AAPP 1110 ELSE 1111 CALL CCOPY( M, A( 1, p ), 1, 1112 $ CWORK(N+1), 1 ) 1113 CALL CLASCL( 'G', 0, 0, AAPP, 1114 $ ONE, M, 1, 1115 $ CWORK(N+1), LDA, IERR ) 1116 AAPQ = CDOTC( M, CWORK(N+1), 1, 1117 $ A( 1, q ), 1 ) / AAQQ 1118 END IF 1119 ELSE 1120 IF( AAPP.GE.AAQQ ) THEN 1121 ROTOK = AAPP.LE.( AAQQ / SMALL ) 1122 ELSE 1123 ROTOK = AAQQ.LE.( AAPP / SMALL ) 1124 END IF 1125 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 1126 AAPQ = ( CDOTC( M, A( 1, p ), 1, 1127 $ A( 1, q ), 1 ) / MAX(AAQQ,AAPP) ) 1128 $ / MIN(AAQQ,AAPP) 1129 ELSE 1130 CALL CCOPY( M, A( 1, q ), 1, 1131 $ CWORK(N+1), 1 ) 1132 CALL CLASCL( 'G', 0, 0, AAQQ, 1133 $ ONE, M, 1, 1134 $ CWORK(N+1), LDA, IERR ) 1135 AAPQ = CDOTC( M, A( 1, p ), 1, 1136 $ CWORK(N+1), 1 ) / AAPP 1137 END IF 1138 END IF 1139* 1140* AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q) 1141 AAPQ1 = -ABS(AAPQ) 1142 MXAAPQ = MAX( MXAAPQ, -AAPQ1 ) 1143* 1144* TO rotate or NOT to rotate, THAT is the question ... 1145* 1146 IF( ABS( AAPQ1 ).GT.TOL ) THEN 1147 OMPQ = AAPQ / ABS(AAPQ) 1148 NOTROT = 0 1149*[RTD] ROTATED = ROTATED + 1 1150 PSKIPPED = 0 1151 ISWROT = ISWROT + 1 1152* 1153 IF( ROTOK ) THEN 1154* 1155 AQOAP = AAQQ / AAPP 1156 APOAQ = AAPP / AAQQ 1157 THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1 1158 IF( AAQQ.GT.AAPP0 )THETA = -THETA 1159* 1160 IF( ABS( THETA ).GT.BIGTHETA ) THEN 1161 T = HALF / THETA 1162 CS = ONE 1163 CALL CROT( M, A(1,p), 1, A(1,q), 1, 1164 $ CS, CONJG(OMPQ)*T ) 1165 IF( RSVEC ) THEN 1166 CALL CROT( MVL, V(1,p), 1, 1167 $ V(1,q), 1, CS, CONJG(OMPQ)*T ) 1168 END IF 1169 SVA( q ) = AAQQ*SQRT( MAX( ZERO, 1170 $ ONE+T*APOAQ*AAPQ1 ) ) 1171 AAPP = AAPP*SQRT( MAX( ZERO, 1172 $ ONE-T*AQOAP*AAPQ1 ) ) 1173 MXSINJ = MAX( MXSINJ, ABS( T ) ) 1174 ELSE 1175* 1176* .. choose correct signum for THETA and rotate 1177* 1178 THSIGN = -SIGN( ONE, AAPQ1 ) 1179 IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN 1180 T = ONE / ( THETA+THSIGN* 1181 $ SQRT( ONE+THETA*THETA ) ) 1182 CS = SQRT( ONE / ( ONE+T*T ) ) 1183 SN = T*CS 1184 MXSINJ = MAX( MXSINJ, ABS( SN ) ) 1185 SVA( q ) = AAQQ*SQRT( MAX( ZERO, 1186 $ ONE+T*APOAQ*AAPQ1 ) ) 1187 AAPP = AAPP*SQRT( MAX( ZERO, 1188 $ ONE-T*AQOAP*AAPQ1 ) ) 1189* 1190 CALL CROT( M, A(1,p), 1, A(1,q), 1, 1191 $ CS, CONJG(OMPQ)*SN ) 1192 IF( RSVEC ) THEN 1193 CALL CROT( MVL, V(1,p), 1, 1194 $ V(1,q), 1, CS, CONJG(OMPQ)*SN ) 1195 END IF 1196 END IF 1197 CWORK(p) = -CWORK(q) * OMPQ 1198* 1199 ELSE 1200* .. have to use modified Gram-Schmidt like transformation 1201 IF( AAPP.GT.AAQQ ) THEN 1202 CALL CCOPY( M, A( 1, p ), 1, 1203 $ CWORK(N+1), 1 ) 1204 CALL CLASCL( 'G', 0, 0, AAPP, ONE, 1205 $ M, 1, CWORK(N+1),LDA, 1206 $ IERR ) 1207 CALL CLASCL( 'G', 0, 0, AAQQ, ONE, 1208 $ M, 1, A( 1, q ), LDA, 1209 $ IERR ) 1210 CALL CAXPY( M, -AAPQ, CWORK(N+1), 1211 $ 1, A( 1, q ), 1 ) 1212 CALL CLASCL( 'G', 0, 0, ONE, AAQQ, 1213 $ M, 1, A( 1, q ), LDA, 1214 $ IERR ) 1215 SVA( q ) = AAQQ*SQRT( MAX( ZERO, 1216 $ ONE-AAPQ1*AAPQ1 ) ) 1217 MXSINJ = MAX( MXSINJ, SFMIN ) 1218 ELSE 1219 CALL CCOPY( M, A( 1, q ), 1, 1220 $ CWORK(N+1), 1 ) 1221 CALL CLASCL( 'G', 0, 0, AAQQ, ONE, 1222 $ M, 1, CWORK(N+1),LDA, 1223 $ IERR ) 1224 CALL CLASCL( 'G', 0, 0, AAPP, ONE, 1225 $ M, 1, A( 1, p ), LDA, 1226 $ IERR ) 1227 CALL CAXPY( M, -CONJG(AAPQ), 1228 $ CWORK(N+1), 1, A( 1, p ), 1 ) 1229 CALL CLASCL( 'G', 0, 0, ONE, AAPP, 1230 $ M, 1, A( 1, p ), LDA, 1231 $ IERR ) 1232 SVA( p ) = AAPP*SQRT( MAX( ZERO, 1233 $ ONE-AAPQ1*AAPQ1 ) ) 1234 MXSINJ = MAX( MXSINJ, SFMIN ) 1235 END IF 1236 END IF 1237* END IF ROTOK THEN ... ELSE 1238* 1239* In the case of cancellation in updating SVA(q), SVA(p) 1240* .. recompute SVA(q), SVA(p) 1241 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 1242 $ THEN 1243 IF( ( AAQQ.LT.ROOTBIG ) .AND. 1244 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN 1245 SVA( q ) = SCNRM2( M, A( 1, q ), 1) 1246 ELSE 1247 T = ZERO 1248 AAQQ = ONE 1249 CALL CLASSQ( M, A( 1, q ), 1, T, 1250 $ AAQQ ) 1251 SVA( q ) = T*SQRT( AAQQ ) 1252 END IF 1253 END IF 1254 IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN 1255 IF( ( AAPP.LT.ROOTBIG ) .AND. 1256 $ ( AAPP.GT.ROOTSFMIN ) ) THEN 1257 AAPP = SCNRM2( M, A( 1, p ), 1 ) 1258 ELSE 1259 T = ZERO 1260 AAPP = ONE 1261 CALL CLASSQ( M, A( 1, p ), 1, T, 1262 $ AAPP ) 1263 AAPP = T*SQRT( AAPP ) 1264 END IF 1265 SVA( p ) = AAPP 1266 END IF 1267* end of OK rotation 1268 ELSE 1269 NOTROT = NOTROT + 1 1270*[RTD] SKIPPED = SKIPPED + 1 1271 PSKIPPED = PSKIPPED + 1 1272 IJBLSK = IJBLSK + 1 1273 END IF 1274 ELSE 1275 NOTROT = NOTROT + 1 1276 PSKIPPED = PSKIPPED + 1 1277 IJBLSK = IJBLSK + 1 1278 END IF 1279* 1280 IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) 1281 $ THEN 1282 SVA( p ) = AAPP 1283 NOTROT = 0 1284 GO TO 2011 1285 END IF 1286 IF( ( i.LE.SWBAND ) .AND. 1287 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN 1288 AAPP = -AAPP 1289 NOTROT = 0 1290 GO TO 2203 1291 END IF 1292* 1293 2200 CONTINUE 1294* end of the q-loop 1295 2203 CONTINUE 1296* 1297 SVA( p ) = AAPP 1298* 1299 ELSE 1300* 1301 IF( AAPP.EQ.ZERO )NOTROT = NOTROT + 1302 $ MIN( jgl+KBL-1, N ) - jgl + 1 1303 IF( AAPP.LT.ZERO )NOTROT = 0 1304* 1305 END IF 1306* 1307 2100 CONTINUE 1308* end of the p-loop 1309 2010 CONTINUE 1310* end of the jbc-loop 1311 2011 CONTINUE 1312*2011 bailed out of the jbc-loop 1313 DO 2012 p = igl, MIN( igl+KBL-1, N ) 1314 SVA( p ) = ABS( SVA( p ) ) 1315 2012 CONTINUE 1316*** 1317 2000 CONTINUE 1318*2000 :: end of the ibr-loop 1319* 1320* .. update SVA(N) 1321 IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) 1322 $ THEN 1323 SVA( N ) = SCNRM2( M, A( 1, N ), 1 ) 1324 ELSE 1325 T = ZERO 1326 AAPP = ONE 1327 CALL CLASSQ( M, A( 1, N ), 1, T, AAPP ) 1328 SVA( N ) = T*SQRT( AAPP ) 1329 END IF 1330* 1331* Additional steering devices 1332* 1333 IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. 1334 $ ( ISWROT.LE.N ) ) )SWBAND = i 1335* 1336 IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( REAL( N ) )* 1337 $ TOL ) .AND. ( REAL( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN 1338 GO TO 1994 1339 END IF 1340* 1341 IF( NOTROT.GE.EMPTSW )GO TO 1994 1342* 1343 1993 CONTINUE 1344* end i=1:NSWEEP loop 1345* 1346* #:( Reaching this point means that the procedure has not converged. 1347 INFO = NSWEEP - 1 1348 GO TO 1995 1349* 1350 1994 CONTINUE 1351* #:) Reaching this point means numerical convergence after the i-th 1352* sweep. 1353* 1354 INFO = 0 1355* #:) INFO = 0 confirms successful iterations. 1356 1995 CONTINUE 1357* 1358* Sort the singular values and find how many are above 1359* the underflow threshold. 1360* 1361 N2 = 0 1362 N4 = 0 1363 DO 5991 p = 1, N - 1 1364 q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1 1365 IF( p.NE.q ) THEN 1366 TEMP1 = SVA( p ) 1367 SVA( p ) = SVA( q ) 1368 SVA( q ) = TEMP1 1369 CALL CSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 1370 IF( RSVEC )CALL CSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 ) 1371 END IF 1372 IF( SVA( p ).NE.ZERO ) THEN 1373 N4 = N4 + 1 1374 IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1 1375 END IF 1376 5991 CONTINUE 1377 IF( SVA( N ).NE.ZERO ) THEN 1378 N4 = N4 + 1 1379 IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1 1380 END IF 1381* 1382* Normalize the left singular vectors. 1383* 1384 IF( LSVEC .OR. UCTOL ) THEN 1385 DO 1998 p = 1, N4 1386* CALL CSSCAL( M, ONE / SVA( p ), A( 1, p ), 1 ) 1387 CALL CLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR ) 1388 1998 CONTINUE 1389 END IF 1390* 1391* Scale the product of Jacobi rotations. 1392* 1393 IF( RSVEC ) THEN 1394 DO 2399 p = 1, N 1395 TEMP1 = ONE / SCNRM2( MVL, V( 1, p ), 1 ) 1396 CALL CSSCAL( MVL, TEMP1, V( 1, p ), 1 ) 1397 2399 CONTINUE 1398 END IF 1399* 1400* Undo scaling, if necessary (and possible). 1401 IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL ) ) ) 1402 $ .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT. 1403 $ ( SFMIN / SKL ) ) ) ) THEN 1404 DO 2400 p = 1, N 1405 SVA( P ) = SKL*SVA( P ) 1406 2400 CONTINUE 1407 SKL = ONE 1408 END IF 1409* 1410 RWORK( 1 ) = SKL 1411* The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE 1412* then some of the singular values may overflow or underflow and 1413* the spectrum is given in this factored representation. 1414* 1415 RWORK( 2 ) = REAL( N4 ) 1416* N4 is the number of computed nonzero singular values of A. 1417* 1418 RWORK( 3 ) = REAL( N2 ) 1419* N2 is the number of singular values of A greater than SFMIN. 1420* If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers 1421* that may carry some information. 1422* 1423 RWORK( 4 ) = REAL( i ) 1424* i is the index of the last sweep before declaring convergence. 1425* 1426 RWORK( 5 ) = MXAAPQ 1427* MXAAPQ is the largest absolute value of scaled pivots in the 1428* last sweep 1429* 1430 RWORK( 6 ) = MXSINJ 1431* MXSINJ is the largest absolute value of the sines of Jacobi angles 1432* in the last sweep 1433* 1434 RETURN 1435* .. 1436* .. END OF CGESVJ 1437* .. 1438 END 1439* 1440