Annotation of rpl/lapack/lapack/dtrsna.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
! 2: $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
! 3: $ INFO )
! 4: *
! 5: * -- LAPACK routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * November 2006
! 9: *
! 10: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
! 11: *
! 12: * .. Scalar Arguments ..
! 13: CHARACTER HOWMNY, JOB
! 14: INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
! 15: * ..
! 16: * .. Array Arguments ..
! 17: LOGICAL SELECT( * )
! 18: INTEGER IWORK( * )
! 19: DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
! 20: $ VR( LDVR, * ), WORK( LDWORK, * )
! 21: * ..
! 22: *
! 23: * Purpose
! 24: * =======
! 25: *
! 26: * DTRSNA estimates reciprocal condition numbers for specified
! 27: * eigenvalues and/or right eigenvectors of a real upper
! 28: * quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
! 29: * orthogonal).
! 30: *
! 31: * T must be in Schur canonical form (as returned by DHSEQR), that is,
! 32: * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
! 33: * 2-by-2 diagonal block has its diagonal elements equal and its
! 34: * off-diagonal elements of opposite sign.
! 35: *
! 36: * Arguments
! 37: * =========
! 38: *
! 39: * JOB (input) CHARACTER*1
! 40: * Specifies whether condition numbers are required for
! 41: * eigenvalues (S) or eigenvectors (SEP):
! 42: * = 'E': for eigenvalues only (S);
! 43: * = 'V': for eigenvectors only (SEP);
! 44: * = 'B': for both eigenvalues and eigenvectors (S and SEP).
! 45: *
! 46: * HOWMNY (input) CHARACTER*1
! 47: * = 'A': compute condition numbers for all eigenpairs;
! 48: * = 'S': compute condition numbers for selected eigenpairs
! 49: * specified by the array SELECT.
! 50: *
! 51: * SELECT (input) LOGICAL array, dimension (N)
! 52: * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
! 53: * condition numbers are required. To select condition numbers
! 54: * for the eigenpair corresponding to a real eigenvalue w(j),
! 55: * SELECT(j) must be set to .TRUE.. To select condition numbers
! 56: * corresponding to a complex conjugate pair of eigenvalues w(j)
! 57: * and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
! 58: * set to .TRUE..
! 59: * If HOWMNY = 'A', SELECT is not referenced.
! 60: *
! 61: * N (input) INTEGER
! 62: * The order of the matrix T. N >= 0.
! 63: *
! 64: * T (input) DOUBLE PRECISION array, dimension (LDT,N)
! 65: * The upper quasi-triangular matrix T, in Schur canonical form.
! 66: *
! 67: * LDT (input) INTEGER
! 68: * The leading dimension of the array T. LDT >= max(1,N).
! 69: *
! 70: * VL (input) DOUBLE PRECISION array, dimension (LDVL,M)
! 71: * If JOB = 'E' or 'B', VL must contain left eigenvectors of T
! 72: * (or of any Q*T*Q**T with Q orthogonal), corresponding to the
! 73: * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
! 74: * must be stored in consecutive columns of VL, as returned by
! 75: * DHSEIN or DTREVC.
! 76: * If JOB = 'V', VL is not referenced.
! 77: *
! 78: * LDVL (input) INTEGER
! 79: * The leading dimension of the array VL.
! 80: * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
! 81: *
! 82: * VR (input) DOUBLE PRECISION array, dimension (LDVR,M)
! 83: * If JOB = 'E' or 'B', VR must contain right eigenvectors of T
! 84: * (or of any Q*T*Q**T with Q orthogonal), corresponding to the
! 85: * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
! 86: * must be stored in consecutive columns of VR, as returned by
! 87: * DHSEIN or DTREVC.
! 88: * If JOB = 'V', VR is not referenced.
! 89: *
! 90: * LDVR (input) INTEGER
! 91: * The leading dimension of the array VR.
! 92: * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
! 93: *
! 94: * S (output) DOUBLE PRECISION array, dimension (MM)
! 95: * If JOB = 'E' or 'B', the reciprocal condition numbers of the
! 96: * selected eigenvalues, stored in consecutive elements of the
! 97: * array. For a complex conjugate pair of eigenvalues two
! 98: * consecutive elements of S are set to the same value. Thus
! 99: * S(j), SEP(j), and the j-th columns of VL and VR all
! 100: * correspond to the same eigenpair (but not in general the
! 101: * j-th eigenpair, unless all eigenpairs are selected).
! 102: * If JOB = 'V', S is not referenced.
! 103: *
! 104: * SEP (output) DOUBLE PRECISION array, dimension (MM)
! 105: * If JOB = 'V' or 'B', the estimated reciprocal condition
! 106: * numbers of the selected eigenvectors, stored in consecutive
! 107: * elements of the array. For a complex eigenvector two
! 108: * consecutive elements of SEP are set to the same value. If
! 109: * the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
! 110: * is set to 0; this can only occur when the true value would be
! 111: * very small anyway.
! 112: * If JOB = 'E', SEP is not referenced.
! 113: *
! 114: * MM (input) INTEGER
! 115: * The number of elements in the arrays S (if JOB = 'E' or 'B')
! 116: * and/or SEP (if JOB = 'V' or 'B'). MM >= M.
! 117: *
! 118: * M (output) INTEGER
! 119: * The number of elements of the arrays S and/or SEP actually
! 120: * used to store the estimated condition numbers.
! 121: * If HOWMNY = 'A', M is set to N.
! 122: *
! 123: * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,N+6)
! 124: * If JOB = 'E', WORK is not referenced.
! 125: *
! 126: * LDWORK (input) INTEGER
! 127: * The leading dimension of the array WORK.
! 128: * LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
! 129: *
! 130: * IWORK (workspace) INTEGER array, dimension (2*(N-1))
! 131: * If JOB = 'E', IWORK is not referenced.
! 132: *
! 133: * INFO (output) INTEGER
! 134: * = 0: successful exit
! 135: * < 0: if INFO = -i, the i-th argument had an illegal value
! 136: *
! 137: * Further Details
! 138: * ===============
! 139: *
! 140: * The reciprocal of the condition number of an eigenvalue lambda is
! 141: * defined as
! 142: *
! 143: * S(lambda) = |v'*u| / (norm(u)*norm(v))
! 144: *
! 145: * where u and v are the right and left eigenvectors of T corresponding
! 146: * to lambda; v' denotes the conjugate-transpose of v, and norm(u)
! 147: * denotes the Euclidean norm. These reciprocal condition numbers always
! 148: * lie between zero (very badly conditioned) and one (very well
! 149: * conditioned). If n = 1, S(lambda) is defined to be 1.
! 150: *
! 151: * An approximate error bound for a computed eigenvalue W(i) is given by
! 152: *
! 153: * EPS * norm(T) / S(i)
! 154: *
! 155: * where EPS is the machine precision.
! 156: *
! 157: * The reciprocal of the condition number of the right eigenvector u
! 158: * corresponding to lambda is defined as follows. Suppose
! 159: *
! 160: * T = ( lambda c )
! 161: * ( 0 T22 )
! 162: *
! 163: * Then the reciprocal condition number is
! 164: *
! 165: * SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
! 166: *
! 167: * where sigma-min denotes the smallest singular value. We approximate
! 168: * the smallest singular value by the reciprocal of an estimate of the
! 169: * one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
! 170: * defined to be abs(T(1,1)).
! 171: *
! 172: * An approximate error bound for a computed right eigenvector VR(i)
! 173: * is given by
! 174: *
! 175: * EPS * norm(T) / SEP(i)
! 176: *
! 177: * =====================================================================
! 178: *
! 179: * .. Parameters ..
! 180: DOUBLE PRECISION ZERO, ONE, TWO
! 181: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
! 182: * ..
! 183: * .. Local Scalars ..
! 184: LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
! 185: INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
! 186: DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
! 187: $ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
! 188: * ..
! 189: * .. Local Arrays ..
! 190: INTEGER ISAVE( 3 )
! 191: DOUBLE PRECISION DUMMY( 1 )
! 192: * ..
! 193: * .. External Functions ..
! 194: LOGICAL LSAME
! 195: DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
! 196: EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
! 197: * ..
! 198: * .. External Subroutines ..
! 199: EXTERNAL DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
! 200: * ..
! 201: * .. Intrinsic Functions ..
! 202: INTRINSIC ABS, MAX, SQRT
! 203: * ..
! 204: * .. Executable Statements ..
! 205: *
! 206: * Decode and test the input parameters
! 207: *
! 208: WANTBH = LSAME( JOB, 'B' )
! 209: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
! 210: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
! 211: *
! 212: SOMCON = LSAME( HOWMNY, 'S' )
! 213: *
! 214: INFO = 0
! 215: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
! 216: INFO = -1
! 217: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
! 218: INFO = -2
! 219: ELSE IF( N.LT.0 ) THEN
! 220: INFO = -4
! 221: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
! 222: INFO = -6
! 223: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
! 224: INFO = -8
! 225: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
! 226: INFO = -10
! 227: ELSE
! 228: *
! 229: * Set M to the number of eigenpairs for which condition numbers
! 230: * are required, and test MM.
! 231: *
! 232: IF( SOMCON ) THEN
! 233: M = 0
! 234: PAIR = .FALSE.
! 235: DO 10 K = 1, N
! 236: IF( PAIR ) THEN
! 237: PAIR = .FALSE.
! 238: ELSE
! 239: IF( K.LT.N ) THEN
! 240: IF( T( K+1, K ).EQ.ZERO ) THEN
! 241: IF( SELECT( K ) )
! 242: $ M = M + 1
! 243: ELSE
! 244: PAIR = .TRUE.
! 245: IF( SELECT( K ) .OR. SELECT( K+1 ) )
! 246: $ M = M + 2
! 247: END IF
! 248: ELSE
! 249: IF( SELECT( N ) )
! 250: $ M = M + 1
! 251: END IF
! 252: END IF
! 253: 10 CONTINUE
! 254: ELSE
! 255: M = N
! 256: END IF
! 257: *
! 258: IF( MM.LT.M ) THEN
! 259: INFO = -13
! 260: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
! 261: INFO = -16
! 262: END IF
! 263: END IF
! 264: IF( INFO.NE.0 ) THEN
! 265: CALL XERBLA( 'DTRSNA', -INFO )
! 266: RETURN
! 267: END IF
! 268: *
! 269: * Quick return if possible
! 270: *
! 271: IF( N.EQ.0 )
! 272: $ RETURN
! 273: *
! 274: IF( N.EQ.1 ) THEN
! 275: IF( SOMCON ) THEN
! 276: IF( .NOT.SELECT( 1 ) )
! 277: $ RETURN
! 278: END IF
! 279: IF( WANTS )
! 280: $ S( 1 ) = ONE
! 281: IF( WANTSP )
! 282: $ SEP( 1 ) = ABS( T( 1, 1 ) )
! 283: RETURN
! 284: END IF
! 285: *
! 286: * Get machine constants
! 287: *
! 288: EPS = DLAMCH( 'P' )
! 289: SMLNUM = DLAMCH( 'S' ) / EPS
! 290: BIGNUM = ONE / SMLNUM
! 291: CALL DLABAD( SMLNUM, BIGNUM )
! 292: *
! 293: KS = 0
! 294: PAIR = .FALSE.
! 295: DO 60 K = 1, N
! 296: *
! 297: * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
! 298: *
! 299: IF( PAIR ) THEN
! 300: PAIR = .FALSE.
! 301: GO TO 60
! 302: ELSE
! 303: IF( K.LT.N )
! 304: $ PAIR = T( K+1, K ).NE.ZERO
! 305: END IF
! 306: *
! 307: * Determine whether condition numbers are required for the k-th
! 308: * eigenpair.
! 309: *
! 310: IF( SOMCON ) THEN
! 311: IF( PAIR ) THEN
! 312: IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
! 313: $ GO TO 60
! 314: ELSE
! 315: IF( .NOT.SELECT( K ) )
! 316: $ GO TO 60
! 317: END IF
! 318: END IF
! 319: *
! 320: KS = KS + 1
! 321: *
! 322: IF( WANTS ) THEN
! 323: *
! 324: * Compute the reciprocal condition number of the k-th
! 325: * eigenvalue.
! 326: *
! 327: IF( .NOT.PAIR ) THEN
! 328: *
! 329: * Real eigenvalue.
! 330: *
! 331: PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
! 332: RNRM = DNRM2( N, VR( 1, KS ), 1 )
! 333: LNRM = DNRM2( N, VL( 1, KS ), 1 )
! 334: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
! 335: ELSE
! 336: *
! 337: * Complex eigenvalue.
! 338: *
! 339: PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
! 340: PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
! 341: $ 1 )
! 342: PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
! 343: PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
! 344: $ 1 )
! 345: RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
! 346: $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
! 347: LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
! 348: $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
! 349: COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
! 350: S( KS ) = COND
! 351: S( KS+1 ) = COND
! 352: END IF
! 353: END IF
! 354: *
! 355: IF( WANTSP ) THEN
! 356: *
! 357: * Estimate the reciprocal condition number of the k-th
! 358: * eigenvector.
! 359: *
! 360: * Copy the matrix T to the array WORK and swap the diagonal
! 361: * block beginning at T(k,k) to the (1,1) position.
! 362: *
! 363: CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
! 364: IFST = K
! 365: ILST = 1
! 366: CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
! 367: $ WORK( 1, N+1 ), IERR )
! 368: *
! 369: IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
! 370: *
! 371: * Could not swap because blocks not well separated
! 372: *
! 373: SCALE = ONE
! 374: EST = BIGNUM
! 375: ELSE
! 376: *
! 377: * Reordering successful
! 378: *
! 379: IF( WORK( 2, 1 ).EQ.ZERO ) THEN
! 380: *
! 381: * Form C = T22 - lambda*I in WORK(2:N,2:N).
! 382: *
! 383: DO 20 I = 2, N
! 384: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
! 385: 20 CONTINUE
! 386: N2 = 1
! 387: NN = N - 1
! 388: ELSE
! 389: *
! 390: * Triangularize the 2 by 2 block by unitary
! 391: * transformation U = [ cs i*ss ]
! 392: * [ i*ss cs ].
! 393: * such that the (1,1) position of WORK is complex
! 394: * eigenvalue lambda with positive imaginary part. (2,2)
! 395: * position of WORK is the complex eigenvalue lambda
! 396: * with negative imaginary part.
! 397: *
! 398: MU = SQRT( ABS( WORK( 1, 2 ) ) )*
! 399: $ SQRT( ABS( WORK( 2, 1 ) ) )
! 400: DELTA = DLAPY2( MU, WORK( 2, 1 ) )
! 401: CS = MU / DELTA
! 402: SN = -WORK( 2, 1 ) / DELTA
! 403: *
! 404: * Form
! 405: *
! 406: * C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
! 407: * [ mu ]
! 408: * [ .. ]
! 409: * [ .. ]
! 410: * [ mu ]
! 411: * where C' is conjugate transpose of complex matrix C,
! 412: * and RWORK is stored starting in the N+1-st column of
! 413: * WORK.
! 414: *
! 415: DO 30 J = 3, N
! 416: WORK( 2, J ) = CS*WORK( 2, J )
! 417: WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
! 418: 30 CONTINUE
! 419: WORK( 2, 2 ) = ZERO
! 420: *
! 421: WORK( 1, N+1 ) = TWO*MU
! 422: DO 40 I = 2, N - 1
! 423: WORK( I, N+1 ) = SN*WORK( 1, I+1 )
! 424: 40 CONTINUE
! 425: N2 = 2
! 426: NN = 2*( N-1 )
! 427: END IF
! 428: *
! 429: * Estimate norm(inv(C'))
! 430: *
! 431: EST = ZERO
! 432: KASE = 0
! 433: 50 CONTINUE
! 434: CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
! 435: $ EST, KASE, ISAVE )
! 436: IF( KASE.NE.0 ) THEN
! 437: IF( KASE.EQ.1 ) THEN
! 438: IF( N2.EQ.1 ) THEN
! 439: *
! 440: * Real eigenvalue: solve C'*x = scale*c.
! 441: *
! 442: CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
! 443: $ LDWORK, DUMMY, DUMM, SCALE,
! 444: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
! 445: $ IERR )
! 446: ELSE
! 447: *
! 448: * Complex eigenvalue: solve
! 449: * C'*(p+iq) = scale*(c+id) in real arithmetic.
! 450: *
! 451: CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
! 452: $ LDWORK, WORK( 1, N+1 ), MU, SCALE,
! 453: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
! 454: $ IERR )
! 455: END IF
! 456: ELSE
! 457: IF( N2.EQ.1 ) THEN
! 458: *
! 459: * Real eigenvalue: solve C*x = scale*c.
! 460: *
! 461: CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
! 462: $ LDWORK, DUMMY, DUMM, SCALE,
! 463: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
! 464: $ IERR )
! 465: ELSE
! 466: *
! 467: * Complex eigenvalue: solve
! 468: * C*(p+iq) = scale*(c+id) in real arithmetic.
! 469: *
! 470: CALL DLAQTR( .FALSE., .FALSE., N-1,
! 471: $ WORK( 2, 2 ), LDWORK,
! 472: $ WORK( 1, N+1 ), MU, SCALE,
! 473: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
! 474: $ IERR )
! 475: *
! 476: END IF
! 477: END IF
! 478: *
! 479: GO TO 50
! 480: END IF
! 481: END IF
! 482: *
! 483: SEP( KS ) = SCALE / MAX( EST, SMLNUM )
! 484: IF( PAIR )
! 485: $ SEP( KS+1 ) = SEP( KS )
! 486: END IF
! 487: *
! 488: IF( PAIR )
! 489: $ KS = KS + 1
! 490: *
! 491: 60 CONTINUE
! 492: RETURN
! 493: *
! 494: * End of DTRSNA
! 495: *
! 496: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack