1*> \brief \b DBDSVDX 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DBDSVDX + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, 22* $ NS, S, Z, LDZ, WORK, IWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBZ, RANGE, UPLO 26* INTEGER IL, INFO, IU, LDZ, N, NS 27* DOUBLE PRECISION VL, VU 28* .. 29* .. Array Arguments .. 30* INTEGER IWORK( * ) 31* DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), 32* Z( LDZ, * ) 33* .. 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> DBDSVDX computes the singular value decomposition (SVD) of a real 41*> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT, 42*> where S is a diagonal matrix with non-negative diagonal elements 43*> (the singular values of B), and U and VT are orthogonal matrices 44*> of left and right singular vectors, respectively. 45*> 46*> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ] 47*> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the 48*> singular value decompositon of B through the eigenvalues and 49*> eigenvectors of the N*2-by-N*2 tridiagonal matrix 50*> 51*> | 0 d_1 | 52*> | d_1 0 e_1 | 53*> TGK = | e_1 0 d_2 | 54*> | d_2 . . | 55*> | . . . | 56*> 57*> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then 58*> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) / 59*> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and 60*> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ]. 61*> 62*> Given a TGK matrix, one can either a) compute -s,-v and change signs 63*> so that the singular values (and corresponding vectors) are already in 64*> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder 65*> the values (and corresponding vectors). DBDSVDX implements a) by 66*> calling DSTEVX (bisection plus inverse iteration, to be replaced 67*> with a version of the Multiple Relative Robust Representation 68*> algorithm. (See P. Willems and B. Lang, A framework for the MR^3 69*> algorithm: theory and implementation, SIAM J. Sci. Comput., 70*> 35:740-766, 2013.) 71*> \endverbatim 72* 73* Arguments: 74* ========== 75* 76*> \param[in] UPLO 77*> \verbatim 78*> UPLO is CHARACTER*1 79*> = 'U': B is upper bidiagonal; 80*> = 'L': B is lower bidiagonal. 81*> \endverbatim 82*> 83*> \param[in] JOBZ 84*> \verbatim 85*> JOBZ is CHARACTER*1 86*> = 'N': Compute singular values only; 87*> = 'V': Compute singular values and singular vectors. 88*> \endverbatim 89*> 90*> \param[in] RANGE 91*> \verbatim 92*> RANGE is CHARACTER*1 93*> = 'A': all singular values will be found. 94*> = 'V': all singular values in the half-open interval [VL,VU) 95*> will be found. 96*> = 'I': the IL-th through IU-th singular values will be found. 97*> \endverbatim 98*> 99*> \param[in] N 100*> \verbatim 101*> N is INTEGER 102*> The order of the bidiagonal matrix. N >= 0. 103*> \endverbatim 104*> 105*> \param[in] D 106*> \verbatim 107*> D is DOUBLE PRECISION array, dimension (N) 108*> The n diagonal elements of the bidiagonal matrix B. 109*> \endverbatim 110*> 111*> \param[in] E 112*> \verbatim 113*> E is DOUBLE PRECISION array, dimension (max(1,N-1)) 114*> The (n-1) superdiagonal elements of the bidiagonal matrix 115*> B in elements 1 to N-1. 116*> \endverbatim 117*> 118*> \param[in] VL 119*> \verbatim 120*> VL is DOUBLE PRECISION 121*> If RANGE='V', the lower bound of the interval to 122*> be searched for singular values. VU > VL. 123*> Not referenced if RANGE = 'A' or 'I'. 124*> \endverbatim 125*> 126*> \param[in] VU 127*> \verbatim 128*> VU is DOUBLE PRECISION 129*> If RANGE='V', the upper bound of the interval to 130*> be searched for singular values. VU > VL. 131*> Not referenced if RANGE = 'A' or 'I'. 132*> \endverbatim 133*> 134*> \param[in] IL 135*> \verbatim 136*> IL is INTEGER 137*> If RANGE='I', the index of the 138*> smallest singular value to be returned. 139*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. 140*> Not referenced if RANGE = 'A' or 'V'. 141*> \endverbatim 142*> 143*> \param[in] IU 144*> \verbatim 145*> IU is INTEGER 146*> If RANGE='I', the index of the 147*> largest singular value to be returned. 148*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0. 149*> Not referenced if RANGE = 'A' or 'V'. 150*> \endverbatim 151*> 152*> \param[out] NS 153*> \verbatim 154*> NS is INTEGER 155*> The total number of singular values found. 0 <= NS <= N. 156*> If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1. 157*> \endverbatim 158*> 159*> \param[out] S 160*> \verbatim 161*> S is DOUBLE PRECISION array, dimension (N) 162*> The first NS elements contain the selected singular values in 163*> ascending order. 164*> \endverbatim 165*> 166*> \param[out] Z 167*> \verbatim 168*> Z is DOUBLE PRECISION array, dimension (2*N,K) 169*> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z 170*> contain the singular vectors of the matrix B corresponding to 171*> the selected singular values, with U in rows 1 to N and V 172*> in rows N+1 to N*2, i.e. 173*> Z = [ U ] 174*> [ V ] 175*> If JOBZ = 'N', then Z is not referenced. 176*> Note: The user must ensure that at least K = NS+1 columns are 177*> supplied in the array Z; if RANGE = 'V', the exact value of 178*> NS is not known in advance and an upper bound must be used. 179*> \endverbatim 180*> 181*> \param[in] LDZ 182*> \verbatim 183*> LDZ is INTEGER 184*> The leading dimension of the array Z. LDZ >= 1, and if 185*> JOBZ = 'V', LDZ >= max(2,N*2). 186*> \endverbatim 187*> 188*> \param[out] WORK 189*> \verbatim 190*> WORK is DOUBLE PRECISION array, dimension (14*N) 191*> \endverbatim 192*> 193*> \param[out] IWORK 194*> \verbatim 195*> IWORK is INTEGER array, dimension (12*N) 196*> If JOBZ = 'V', then if INFO = 0, the first NS elements of 197*> IWORK are zero. If INFO > 0, then IWORK contains the indices 198*> of the eigenvectors that failed to converge in DSTEVX. 199*> \endverbatim 200*> 201*> \param[out] INFO 202*> \verbatim 203*> INFO is INTEGER 204*> = 0: successful exit 205*> < 0: if INFO = -i, the i-th argument had an illegal value 206*> > 0: if INFO = i, then i eigenvectors failed to converge 207*> in DSTEVX. The indices of the eigenvectors 208*> (as returned by DSTEVX) are stored in the 209*> array IWORK. 210*> if INFO = N*2 + 1, an internal error occurred. 211*> \endverbatim 212* 213* Authors: 214* ======== 215* 216*> \author Univ. of Tennessee 217*> \author Univ. of California Berkeley 218*> \author Univ. of Colorado Denver 219*> \author NAG Ltd. 220* 221*> \ingroup doubleOTHEReigen 222* 223* ===================================================================== 224 SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, 225 $ NS, S, Z, LDZ, WORK, IWORK, INFO) 226* 227* -- LAPACK driver routine -- 228* -- LAPACK is a software package provided by Univ. of Tennessee, -- 229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 230* 231* .. Scalar Arguments .. 232 CHARACTER JOBZ, RANGE, UPLO 233 INTEGER IL, INFO, IU, LDZ, N, NS 234 DOUBLE PRECISION VL, VU 235* .. 236* .. Array Arguments .. 237 INTEGER IWORK( * ) 238 DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ), 239 $ Z( LDZ, * ) 240* .. 241* 242* ===================================================================== 243* 244* .. Parameters .. 245 DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH 246 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0, 247 $ HNDRD = 100.0D0, MEIGTH = -0.1250D0 ) 248 DOUBLE PRECISION FUDGE 249 PARAMETER ( FUDGE = 2.0D0 ) 250* .. 251* .. Local Scalars .. 252 CHARACTER RNGVX 253 LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ 254 INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR, 255 $ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV, 256 $ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K, 257 $ NTGK, NRU, NRV, NSL 258 DOUBLE PRECISION ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX, 259 $ SMIN, SQRT2, THRESH, TOL, ULP, 260 $ VLTGK, VUTGK, ZJTJI 261* .. 262* .. External Functions .. 263 LOGICAL LSAME 264 INTEGER IDAMAX 265 DOUBLE PRECISION DDOT, DLAMCH, DNRM2 266 EXTERNAL IDAMAX, LSAME, DAXPY, DDOT, DLAMCH, DNRM2 267* .. 268* .. External Subroutines .. 269 EXTERNAL DSTEVX, DCOPY, DLASET, DSCAL, DSWAP, XERBLA 270* .. 271* .. Intrinsic Functions .. 272 INTRINSIC ABS, DBLE, SIGN, SQRT 273* .. 274* .. Executable Statements .. 275* 276* Test the input parameters. 277* 278 ALLSV = LSAME( RANGE, 'A' ) 279 VALSV = LSAME( RANGE, 'V' ) 280 INDSV = LSAME( RANGE, 'I' ) 281 WANTZ = LSAME( JOBZ, 'V' ) 282 LOWER = LSAME( UPLO, 'L' ) 283* 284 INFO = 0 285 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LOWER ) THEN 286 INFO = -1 287 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 288 INFO = -2 289 ELSE IF( .NOT.( ALLSV .OR. VALSV .OR. INDSV ) ) THEN 290 INFO = -3 291 ELSE IF( N.LT.0 ) THEN 292 INFO = -4 293 ELSE IF( N.GT.0 ) THEN 294 IF( VALSV ) THEN 295 IF( VL.LT.ZERO ) THEN 296 INFO = -7 297 ELSE IF( VU.LE.VL ) THEN 298 INFO = -8 299 END IF 300 ELSE IF( INDSV ) THEN 301 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 302 INFO = -9 303 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 304 INFO = -10 305 END IF 306 END IF 307 END IF 308 IF( INFO.EQ.0 ) THEN 309 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N*2 ) ) INFO = -14 310 END IF 311* 312 IF( INFO.NE.0 ) THEN 313 CALL XERBLA( 'DBDSVDX', -INFO ) 314 RETURN 315 END IF 316* 317* Quick return if possible (N.LE.1) 318* 319 NS = 0 320 IF( N.EQ.0 ) RETURN 321* 322 IF( N.EQ.1 ) THEN 323 IF( ALLSV .OR. INDSV ) THEN 324 NS = 1 325 S( 1 ) = ABS( D( 1 ) ) 326 ELSE 327 IF( VL.LT.ABS( D( 1 ) ) .AND. VU.GE.ABS( D( 1 ) ) ) THEN 328 NS = 1 329 S( 1 ) = ABS( D( 1 ) ) 330 END IF 331 END IF 332 IF( WANTZ ) THEN 333 Z( 1, 1 ) = SIGN( ONE, D( 1 ) ) 334 Z( 2, 1 ) = ONE 335 END IF 336 RETURN 337 END IF 338* 339 ABSTOL = 2*DLAMCH( 'Safe Minimum' ) 340 ULP = DLAMCH( 'Precision' ) 341 EPS = DLAMCH( 'Epsilon' ) 342 SQRT2 = SQRT( 2.0D0 ) 343 ORTOL = SQRT( ULP ) 344* 345* Criterion for splitting is taken from DBDSQR when singular 346* values are computed to relative accuracy TOL. (See J. Demmel and 347* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM 348* J. Sci. and Stat. Comput., 11:873–912, 1990.) 349* 350 TOL = MAX( TEN, MIN( HNDRD, EPS**MEIGTH ) )*EPS 351* 352* Compute approximate maximum, minimum singular values. 353* 354 I = IDAMAX( N, D, 1 ) 355 SMAX = ABS( D( I ) ) 356 I = IDAMAX( N-1, E, 1 ) 357 SMAX = MAX( SMAX, ABS( E( I ) ) ) 358* 359* Compute threshold for neglecting D's and E's. 360* 361 SMIN = ABS( D( 1 ) ) 362 IF( SMIN.NE.ZERO ) THEN 363 MU = SMIN 364 DO I = 2, N 365 MU = ABS( D( I ) )*( MU / ( MU+ABS( E( I-1 ) ) ) ) 366 SMIN = MIN( SMIN, MU ) 367 IF( SMIN.EQ.ZERO ) EXIT 368 END DO 369 END IF 370 SMIN = SMIN / SQRT( DBLE( N ) ) 371 THRESH = TOL*SMIN 372* 373* Check for zeros in D and E (splits), i.e. submatrices. 374* 375 DO I = 1, N-1 376 IF( ABS( D( I ) ).LE.THRESH ) D( I ) = ZERO 377 IF( ABS( E( I ) ).LE.THRESH ) E( I ) = ZERO 378 END DO 379 IF( ABS( D( N ) ).LE.THRESH ) D( N ) = ZERO 380* 381* Pointers for arrays used by DSTEVX. 382* 383 IDTGK = 1 384 IETGK = IDTGK + N*2 385 ITEMP = IETGK + N*2 386 IIFAIL = 1 387 IIWORK = IIFAIL + N*2 388* 389* Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode. 390* VL,VU or IL,IU are redefined to conform to implementation a) 391* described in the leading comments. 392* 393 ILTGK = 0 394 IUTGK = 0 395 VLTGK = ZERO 396 VUTGK = ZERO 397* 398 IF( ALLSV ) THEN 399* 400* All singular values will be found. We aim at -s (see 401* leading comments) with RNGVX = 'I'. IL and IU are set 402* later (as ILTGK and IUTGK) according to the dimension 403* of the active submatrix. 404* 405 RNGVX = 'I' 406 IF( WANTZ ) CALL DLASET( 'F', N*2, N+1, ZERO, ZERO, Z, LDZ ) 407 ELSE IF( VALSV ) THEN 408* 409* Find singular values in a half-open interval. We aim 410* at -s (see leading comments) and we swap VL and VU 411* (as VUTGK and VLTGK), changing their signs. 412* 413 RNGVX = 'V' 414 VLTGK = -VU 415 VUTGK = -VL 416 WORK( IDTGK:IDTGK+2*N-1 ) = ZERO 417 CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) 418 CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) 419 CALL DSTEVX( 'N', 'V', N*2, WORK( IDTGK ), WORK( IETGK ), 420 $ VLTGK, VUTGK, ILTGK, ILTGK, ABSTOL, NS, S, 421 $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), 422 $ IWORK( IIFAIL ), INFO ) 423 IF( NS.EQ.0 ) THEN 424 RETURN 425 ELSE 426 IF( WANTZ ) CALL DLASET( 'F', N*2, NS, ZERO, ZERO, Z, LDZ ) 427 END IF 428 ELSE IF( INDSV ) THEN 429* 430* Find the IL-th through the IU-th singular values. We aim 431* at -s (see leading comments) and indices are mapped into 432* values, therefore mimicking DSTEBZ, where 433* 434* GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN 435* GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN 436* 437 ILTGK = IL 438 IUTGK = IU 439 RNGVX = 'V' 440 WORK( IDTGK:IDTGK+2*N-1 ) = ZERO 441 CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) 442 CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) 443 CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), 444 $ VLTGK, VLTGK, ILTGK, ILTGK, ABSTOL, NS, S, 445 $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), 446 $ IWORK( IIFAIL ), INFO ) 447 VLTGK = S( 1 ) - FUDGE*SMAX*ULP*N 448 WORK( IDTGK:IDTGK+2*N-1 ) = ZERO 449 CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) 450 CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) 451 CALL DSTEVX( 'N', 'I', N*2, WORK( IDTGK ), WORK( IETGK ), 452 $ VUTGK, VUTGK, IUTGK, IUTGK, ABSTOL, NS, S, 453 $ Z, LDZ, WORK( ITEMP ), IWORK( IIWORK ), 454 $ IWORK( IIFAIL ), INFO ) 455 VUTGK = S( 1 ) + FUDGE*SMAX*ULP*N 456 VUTGK = MIN( VUTGK, ZERO ) 457* 458* If VLTGK=VUTGK, DSTEVX returns an error message, 459* so if needed we change VUTGK slightly. 460* 461 IF( VLTGK.EQ.VUTGK ) VLTGK = VLTGK - TOL 462* 463 IF( WANTZ ) CALL DLASET( 'F', N*2, IU-IL+1, ZERO, ZERO, Z, LDZ) 464 END IF 465* 466* Initialize variables and pointers for S, Z, and WORK. 467* 468* NRU, NRV: number of rows in U and V for the active submatrix 469* IDBEG, ISBEG: offsets for the entries of D and S 470* IROWZ, ICOLZ: offsets for the rows and columns of Z 471* IROWU, IROWV: offsets for the rows of U and V 472* 473 NS = 0 474 NRU = 0 475 NRV = 0 476 IDBEG = 1 477 ISBEG = 1 478 IROWZ = 1 479 ICOLZ = 1 480 IROWU = 2 481 IROWV = 1 482 SPLIT = .FALSE. 483 SVEQ0 = .FALSE. 484* 485* Form the tridiagonal TGK matrix. 486* 487 S( 1:N ) = ZERO 488 WORK( IETGK+2*N-1 ) = ZERO 489 WORK( IDTGK:IDTGK+2*N-1 ) = ZERO 490 CALL DCOPY( N, D, 1, WORK( IETGK ), 2 ) 491 CALL DCOPY( N-1, E, 1, WORK( IETGK+1 ), 2 ) 492* 493* 494* Check for splits in two levels, outer level 495* in E and inner level in D. 496* 497 DO IEPTR = 2, N*2, 2 498 IF( WORK( IETGK+IEPTR-1 ).EQ.ZERO ) THEN 499* 500* Split in E (this piece of B is square) or bottom 501* of the (input bidiagonal) matrix. 502* 503 ISPLT = IDBEG 504 IDEND = IEPTR - 1 505 DO IDPTR = IDBEG, IDEND, 2 506 IF( WORK( IETGK+IDPTR-1 ).EQ.ZERO ) THEN 507* 508* Split in D (rectangular submatrix). Set the number 509* of rows in U and V (NRU and NRV) accordingly. 510* 511 IF( IDPTR.EQ.IDBEG ) THEN 512* 513* D=0 at the top. 514* 515 SVEQ0 = .TRUE. 516 IF( IDBEG.EQ.IDEND) THEN 517 NRU = 1 518 NRV = 1 519 END IF 520 ELSE IF( IDPTR.EQ.IDEND ) THEN 521* 522* D=0 at the bottom. 523* 524 SVEQ0 = .TRUE. 525 NRU = (IDEND-ISPLT)/2 + 1 526 NRV = NRU 527 IF( ISPLT.NE.IDBEG ) THEN 528 NRU = NRU + 1 529 END IF 530 ELSE 531 IF( ISPLT.EQ.IDBEG ) THEN 532* 533* Split: top rectangular submatrix. 534* 535 NRU = (IDPTR-IDBEG)/2 536 NRV = NRU + 1 537 ELSE 538* 539* Split: middle square submatrix. 540* 541 NRU = (IDPTR-ISPLT)/2 + 1 542 NRV = NRU 543 END IF 544 END IF 545 ELSE IF( IDPTR.EQ.IDEND ) THEN 546* 547* Last entry of D in the active submatrix. 548* 549 IF( ISPLT.EQ.IDBEG ) THEN 550* 551* No split (trivial case). 552* 553 NRU = (IDEND-IDBEG)/2 + 1 554 NRV = NRU 555 ELSE 556* 557* Split: bottom rectangular submatrix. 558* 559 NRV = (IDEND-ISPLT)/2 + 1 560 NRU = NRV + 1 561 END IF 562 END IF 563* 564 NTGK = NRU + NRV 565* 566 IF( NTGK.GT.0 ) THEN 567* 568* Compute eigenvalues/vectors of the active 569* submatrix according to RANGE: 570* if RANGE='A' (ALLSV) then RNGVX = 'I' 571* if RANGE='V' (VALSV) then RNGVX = 'V' 572* if RANGE='I' (INDSV) then RNGVX = 'V' 573* 574 ILTGK = 1 575 IUTGK = NTGK / 2 576 IF( ALLSV .OR. VUTGK.EQ.ZERO ) THEN 577 IF( SVEQ0 .OR. 578 $ SMIN.LT.EPS .OR. 579 $ MOD(NTGK,2).GT.0 ) THEN 580* Special case: eigenvalue equal to zero or very 581* small, additional eigenvector is needed. 582 IUTGK = IUTGK + 1 583 END IF 584 END IF 585* 586* Workspace needed by DSTEVX: 587* WORK( ITEMP: ): 2*5*NTGK 588* IWORK( 1: ): 2*6*NTGK 589* 590 CALL DSTEVX( JOBZ, RNGVX, NTGK, WORK( IDTGK+ISPLT-1 ), 591 $ WORK( IETGK+ISPLT-1 ), VLTGK, VUTGK, 592 $ ILTGK, IUTGK, ABSTOL, NSL, S( ISBEG ), 593 $ Z( IROWZ,ICOLZ ), LDZ, WORK( ITEMP ), 594 $ IWORK( IIWORK ), IWORK( IIFAIL ), 595 $ INFO ) 596 IF( INFO.NE.0 ) THEN 597* Exit with the error code from DSTEVX. 598 RETURN 599 END IF 600 EMIN = ABS( MAXVAL( S( ISBEG:ISBEG+NSL-1 ) ) ) 601* 602 IF( NSL.GT.0 .AND. WANTZ ) THEN 603* 604* Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:), 605* changing the sign of v as discussed in the leading 606* comments. The norms of u and v may be (slightly) 607* different from 1/sqrt(2) if the corresponding 608* eigenvalues are very small or too close. We check 609* those norms and, if needed, reorthogonalize the 610* vectors. 611* 612 IF( NSL.GT.1 .AND. 613 $ VUTGK.EQ.ZERO .AND. 614 $ MOD(NTGK,2).EQ.0 .AND. 615 $ EMIN.EQ.0 .AND. .NOT.SPLIT ) THEN 616* 617* D=0 at the top or bottom of the active submatrix: 618* one eigenvalue is equal to zero; concatenate the 619* eigenvectors corresponding to the two smallest 620* eigenvalues. 621* 622 Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) = 623 $ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-2 ) + 624 $ Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) 625 Z( IROWZ:IROWZ+NTGK-1,ICOLZ+NSL-1 ) = 626 $ ZERO 627* IF( IUTGK*2.GT.NTGK ) THEN 628* Eigenvalue equal to zero or very small. 629* NSL = NSL - 1 630* END IF 631 END IF 632* 633 DO I = 0, MIN( NSL-1, NRU-1 ) 634 NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 ) 635 IF( NRMU.EQ.ZERO ) THEN 636 INFO = N*2 + 1 637 RETURN 638 END IF 639 CALL DSCAL( NRU, ONE/NRMU, 640 $ Z( IROWU,ICOLZ+I ), 2 ) 641 IF( NRMU.NE.ONE .AND. 642 $ ABS( NRMU-ORTOL )*SQRT2.GT.ONE ) 643 $ THEN 644 DO J = 0, I-1 645 ZJTJI = -DDOT( NRU, Z( IROWU, ICOLZ+J ), 646 $ 2, Z( IROWU, ICOLZ+I ), 2 ) 647 CALL DAXPY( NRU, ZJTJI, 648 $ Z( IROWU, ICOLZ+J ), 2, 649 $ Z( IROWU, ICOLZ+I ), 2 ) 650 END DO 651 NRMU = DNRM2( NRU, Z( IROWU, ICOLZ+I ), 2 ) 652 CALL DSCAL( NRU, ONE/NRMU, 653 $ Z( IROWU,ICOLZ+I ), 2 ) 654 END IF 655 END DO 656 DO I = 0, MIN( NSL-1, NRV-1 ) 657 NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 ) 658 IF( NRMV.EQ.ZERO ) THEN 659 INFO = N*2 + 1 660 RETURN 661 END IF 662 CALL DSCAL( NRV, -ONE/NRMV, 663 $ Z( IROWV,ICOLZ+I ), 2 ) 664 IF( NRMV.NE.ONE .AND. 665 $ ABS( NRMV-ORTOL )*SQRT2.GT.ONE ) 666 $ THEN 667 DO J = 0, I-1 668 ZJTJI = -DDOT( NRV, Z( IROWV, ICOLZ+J ), 669 $ 2, Z( IROWV, ICOLZ+I ), 2 ) 670 CALL DAXPY( NRU, ZJTJI, 671 $ Z( IROWV, ICOLZ+J ), 2, 672 $ Z( IROWV, ICOLZ+I ), 2 ) 673 END DO 674 NRMV = DNRM2( NRV, Z( IROWV, ICOLZ+I ), 2 ) 675 CALL DSCAL( NRV, ONE/NRMV, 676 $ Z( IROWV,ICOLZ+I ), 2 ) 677 END IF 678 END DO 679 IF( VUTGK.EQ.ZERO .AND. 680 $ IDPTR.LT.IDEND .AND. 681 $ MOD(NTGK,2).GT.0 ) THEN 682* 683* D=0 in the middle of the active submatrix (one 684* eigenvalue is equal to zero): save the corresponding 685* eigenvector for later use (when bottom of the 686* active submatrix is reached). 687* 688 SPLIT = .TRUE. 689 Z( IROWZ:IROWZ+NTGK-1,N+1 ) = 690 $ Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) 691 Z( IROWZ:IROWZ+NTGK-1,NS+NSL ) = 692 $ ZERO 693 END IF 694 END IF !** WANTZ **! 695* 696 NSL = MIN( NSL, NRU ) 697 SVEQ0 = .FALSE. 698* 699* Absolute values of the eigenvalues of TGK. 700* 701 DO I = 0, NSL-1 702 S( ISBEG+I ) = ABS( S( ISBEG+I ) ) 703 END DO 704* 705* Update pointers for TGK, S and Z. 706* 707 ISBEG = ISBEG + NSL 708 IROWZ = IROWZ + NTGK 709 ICOLZ = ICOLZ + NSL 710 IROWU = IROWZ 711 IROWV = IROWZ + 1 712 ISPLT = IDPTR + 1 713 NS = NS + NSL 714 NRU = 0 715 NRV = 0 716 END IF !** NTGK.GT.0 **! 717 IF( IROWZ.LT.N*2 .AND. WANTZ ) THEN 718 Z( 1:IROWZ-1, ICOLZ ) = ZERO 719 END IF 720 END DO !** IDPTR loop **! 721 IF( SPLIT .AND. WANTZ ) THEN 722* 723* Bring back eigenvector corresponding 724* to eigenvalue equal to zero. 725* 726 Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) = 727 $ Z( IDBEG:IDEND-NTGK+1,ISBEG-1 ) + 728 $ Z( IDBEG:IDEND-NTGK+1,N+1 ) 729 Z( IDBEG:IDEND-NTGK+1,N+1 ) = 0 730 END IF 731 IROWV = IROWV - 1 732 IROWU = IROWU + 1 733 IDBEG = IEPTR + 1 734 SVEQ0 = .FALSE. 735 SPLIT = .FALSE. 736 END IF !** Check for split in E **! 737 END DO !** IEPTR loop **! 738* 739* Sort the singular values into decreasing order (insertion sort on 740* singular values, but only one transposition per singular vector) 741* 742 DO I = 1, NS-1 743 K = 1 744 SMIN = S( 1 ) 745 DO J = 2, NS + 1 - I 746 IF( S( J ).LE.SMIN ) THEN 747 K = J 748 SMIN = S( J ) 749 END IF 750 END DO 751 IF( K.NE.NS+1-I ) THEN 752 S( K ) = S( NS+1-I ) 753 S( NS+1-I ) = SMIN 754 IF( WANTZ ) CALL DSWAP( N*2, Z( 1,K ), 1, Z( 1,NS+1-I ), 1 ) 755 END IF 756 END DO 757* 758* If RANGE=I, check for singular values/vectors to be discarded. 759* 760 IF( INDSV ) THEN 761 K = IU - IL + 1 762 IF( K.LT.NS ) THEN 763 S( K+1:NS ) = ZERO 764 IF( WANTZ ) Z( 1:N*2,K+1:NS ) = ZERO 765 NS = K 766 END IF 767 END IF 768* 769* Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ). 770* If B is a lower diagonal, swap U and V. 771* 772 IF( WANTZ ) THEN 773 DO I = 1, NS 774 CALL DCOPY( N*2, Z( 1,I ), 1, WORK, 1 ) 775 IF( LOWER ) THEN 776 CALL DCOPY( N, WORK( 2 ), 2, Z( N+1,I ), 1 ) 777 CALL DCOPY( N, WORK( 1 ), 2, Z( 1 ,I ), 1 ) 778 ELSE 779 CALL DCOPY( N, WORK( 2 ), 2, Z( 1 ,I ), 1 ) 780 CALL DCOPY( N, WORK( 1 ), 2, Z( N+1,I ), 1 ) 781 END IF 782 END DO 783 END IF 784* 785 RETURN 786* 787* End of DBDSVDX 788* 789 END 790