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