Annotation of rpl/lapack/lapack/ztgsna.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZTGSNA
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZTGSNA + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsna.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsna.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsna.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
! 22: * LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
! 23: * IWORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER HOWMNY, JOB
! 27: * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * LOGICAL SELECT( * )
! 31: * INTEGER IWORK( * )
! 32: * DOUBLE PRECISION DIF( * ), S( * )
! 33: * COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
! 34: * $ VR( LDVR, * ), WORK( * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: *> ZTGSNA estimates reciprocal condition numbers for specified
! 44: *> eigenvalues and/or eigenvectors of a matrix pair (A, B).
! 45: *>
! 46: *> (A, B) must be in generalized Schur canonical form, that is, A and
! 47: *> B are both upper triangular.
! 48: *> \endverbatim
! 49: *
! 50: * Arguments:
! 51: * ==========
! 52: *
! 53: *> \param[in] JOB
! 54: *> \verbatim
! 55: *> JOB is CHARACTER*1
! 56: *> Specifies whether condition numbers are required for
! 57: *> eigenvalues (S) or eigenvectors (DIF):
! 58: *> = 'E': for eigenvalues only (S);
! 59: *> = 'V': for eigenvectors only (DIF);
! 60: *> = 'B': for both eigenvalues and eigenvectors (S and DIF).
! 61: *> \endverbatim
! 62: *>
! 63: *> \param[in] HOWMNY
! 64: *> \verbatim
! 65: *> HOWMNY is CHARACTER*1
! 66: *> = 'A': compute condition numbers for all eigenpairs;
! 67: *> = 'S': compute condition numbers for selected eigenpairs
! 68: *> specified by the array SELECT.
! 69: *> \endverbatim
! 70: *>
! 71: *> \param[in] SELECT
! 72: *> \verbatim
! 73: *> SELECT is LOGICAL array, dimension (N)
! 74: *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
! 75: *> condition numbers are required. To select condition numbers
! 76: *> for the corresponding j-th eigenvalue and/or eigenvector,
! 77: *> SELECT(j) must be set to .TRUE..
! 78: *> If HOWMNY = 'A', SELECT is not referenced.
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] N
! 82: *> \verbatim
! 83: *> N is INTEGER
! 84: *> The order of the square matrix pair (A, B). N >= 0.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in] A
! 88: *> \verbatim
! 89: *> A is COMPLEX*16 array, dimension (LDA,N)
! 90: *> The upper triangular matrix A in the pair (A,B).
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[in] LDA
! 94: *> \verbatim
! 95: *> LDA is INTEGER
! 96: *> The leading dimension of the array A. LDA >= max(1,N).
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] B
! 100: *> \verbatim
! 101: *> B is COMPLEX*16 array, dimension (LDB,N)
! 102: *> The upper triangular matrix B in the pair (A, B).
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] LDB
! 106: *> \verbatim
! 107: *> LDB is INTEGER
! 108: *> The leading dimension of the array B. LDB >= max(1,N).
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[in] VL
! 112: *> \verbatim
! 113: *> VL is COMPLEX*16 array, dimension (LDVL,M)
! 114: *> IF JOB = 'E' or 'B', VL must contain left eigenvectors of
! 115: *> (A, B), corresponding to the eigenpairs specified by HOWMNY
! 116: *> and SELECT. The eigenvectors must be stored in consecutive
! 117: *> columns of VL, as returned by ZTGEVC.
! 118: *> If JOB = 'V', VL is not referenced.
! 119: *> \endverbatim
! 120: *>
! 121: *> \param[in] LDVL
! 122: *> \verbatim
! 123: *> LDVL is INTEGER
! 124: *> The leading dimension of the array VL. LDVL >= 1; and
! 125: *> If JOB = 'E' or 'B', LDVL >= N.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in] VR
! 129: *> \verbatim
! 130: *> VR is COMPLEX*16 array, dimension (LDVR,M)
! 131: *> IF JOB = 'E' or 'B', VR must contain right eigenvectors of
! 132: *> (A, B), corresponding to the eigenpairs specified by HOWMNY
! 133: *> and SELECT. The eigenvectors must be stored in consecutive
! 134: *> columns of VR, as returned by ZTGEVC.
! 135: *> If JOB = 'V', VR is not referenced.
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[in] LDVR
! 139: *> \verbatim
! 140: *> LDVR is INTEGER
! 141: *> The leading dimension of the array VR. LDVR >= 1;
! 142: *> If JOB = 'E' or 'B', LDVR >= N.
! 143: *> \endverbatim
! 144: *>
! 145: *> \param[out] S
! 146: *> \verbatim
! 147: *> S is DOUBLE PRECISION array, dimension (MM)
! 148: *> If JOB = 'E' or 'B', the reciprocal condition numbers of the
! 149: *> selected eigenvalues, stored in consecutive elements of the
! 150: *> array.
! 151: *> If JOB = 'V', S is not referenced.
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[out] DIF
! 155: *> \verbatim
! 156: *> DIF is DOUBLE PRECISION array, dimension (MM)
! 157: *> If JOB = 'V' or 'B', the estimated reciprocal condition
! 158: *> numbers of the selected eigenvectors, stored in consecutive
! 159: *> elements of the array.
! 160: *> If the eigenvalues cannot be reordered to compute DIF(j),
! 161: *> DIF(j) is set to 0; this can only occur when the true value
! 162: *> would be very small anyway.
! 163: *> For each eigenvalue/vector specified by SELECT, DIF stores
! 164: *> a Frobenius norm-based estimate of Difl.
! 165: *> If JOB = 'E', DIF is not referenced.
! 166: *> \endverbatim
! 167: *>
! 168: *> \param[in] MM
! 169: *> \verbatim
! 170: *> MM is INTEGER
! 171: *> The number of elements in the arrays S and DIF. MM >= M.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[out] M
! 175: *> \verbatim
! 176: *> M is INTEGER
! 177: *> The number of elements of the arrays S and DIF used to store
! 178: *> the specified condition numbers; for each selected eigenvalue
! 179: *> one element is used. If HOWMNY = 'A', M is set to N.
! 180: *> \endverbatim
! 181: *>
! 182: *> \param[out] WORK
! 183: *> \verbatim
! 184: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 185: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 186: *> \endverbatim
! 187: *>
! 188: *> \param[in] LWORK
! 189: *> \verbatim
! 190: *> LWORK is INTEGER
! 191: *> The dimension of the array WORK. LWORK >= max(1,N).
! 192: *> If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
! 193: *> \endverbatim
! 194: *>
! 195: *> \param[out] IWORK
! 196: *> \verbatim
! 197: *> IWORK is INTEGER array, dimension (N+2)
! 198: *> If JOB = 'E', IWORK is not referenced.
! 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: *> \endverbatim
! 207: *
! 208: * Authors:
! 209: * ========
! 210: *
! 211: *> \author Univ. of Tennessee
! 212: *> \author Univ. of California Berkeley
! 213: *> \author Univ. of Colorado Denver
! 214: *> \author NAG Ltd.
! 215: *
! 216: *> \date November 2011
! 217: *
! 218: *> \ingroup complex16OTHERcomputational
! 219: *
! 220: *> \par Further Details:
! 221: * =====================
! 222: *>
! 223: *> \verbatim
! 224: *>
! 225: *> The reciprocal of the condition number of the i-th generalized
! 226: *> eigenvalue w = (a, b) is defined as
! 227: *>
! 228: *> S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
! 229: *>
! 230: *> where u and v are the right and left eigenvectors of (A, B)
! 231: *> corresponding to w; |z| denotes the absolute value of the complex
! 232: *> number, and norm(u) denotes the 2-norm of the vector u. The pair
! 233: *> (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
! 234: *> matrix pair (A, B). If both a and b equal zero, then (A,B) is
! 235: *> singular and S(I) = -1 is returned.
! 236: *>
! 237: *> An approximate error bound on the chordal distance between the i-th
! 238: *> computed generalized eigenvalue w and the corresponding exact
! 239: *> eigenvalue lambda is
! 240: *>
! 241: *> chord(w, lambda) <= EPS * norm(A, B) / S(I),
! 242: *>
! 243: *> where EPS is the machine precision.
! 244: *>
! 245: *> The reciprocal of the condition number of the right eigenvector u
! 246: *> and left eigenvector v corresponding to the generalized eigenvalue w
! 247: *> is defined as follows. Suppose
! 248: *>
! 249: *> (A, B) = ( a * ) ( b * ) 1
! 250: *> ( 0 A22 ),( 0 B22 ) n-1
! 251: *> 1 n-1 1 n-1
! 252: *>
! 253: *> Then the reciprocal condition number DIF(I) is
! 254: *>
! 255: *> Difl[(a, b), (A22, B22)] = sigma-min( Zl )
! 256: *>
! 257: *> where sigma-min(Zl) denotes the smallest singular value of
! 258: *>
! 259: *> Zl = [ kron(a, In-1) -kron(1, A22) ]
! 260: *> [ kron(b, In-1) -kron(1, B22) ].
! 261: *>
! 262: *> Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
! 263: *> transpose of X. kron(X, Y) is the Kronecker product between the
! 264: *> matrices X and Y.
! 265: *>
! 266: *> We approximate the smallest singular value of Zl with an upper
! 267: *> bound. This is done by ZLATDF.
! 268: *>
! 269: *> An approximate error bound for a computed eigenvector VL(i) or
! 270: *> VR(i) is given by
! 271: *>
! 272: *> EPS * norm(A, B) / DIF(i).
! 273: *>
! 274: *> See ref. [2-3] for more details and further references.
! 275: *> \endverbatim
! 276: *
! 277: *> \par Contributors:
! 278: * ==================
! 279: *>
! 280: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
! 281: *> Umea University, S-901 87 Umea, Sweden.
! 282: *
! 283: *> \par References:
! 284: * ================
! 285: *>
! 286: *> \verbatim
! 287: *>
! 288: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
! 289: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
! 290: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
! 291: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
! 292: *>
! 293: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
! 294: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
! 295: *> Estimation: Theory, Algorithms and Software, Report
! 296: *> UMINF - 94.04, Department of Computing Science, Umea University,
! 297: *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
! 298: *> To appear in Numerical Algorithms, 1996.
! 299: *>
! 300: *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
! 301: *> for Solving the Generalized Sylvester Equation and Estimating the
! 302: *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
! 303: *> Department of Computing Science, Umea University, S-901 87 Umea,
! 304: *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
! 305: *> Note 75.
! 306: *> To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
! 307: *> \endverbatim
! 308: *>
! 309: * =====================================================================
1.1 bertrand 310: SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
311: $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
312: $ IWORK, INFO )
313: *
1.9 ! bertrand 314: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 315: * -- LAPACK is a software package provided by Univ. of Tennessee, --
316: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 317: * November 2011
1.1 bertrand 318: *
319: * .. Scalar Arguments ..
320: CHARACTER HOWMNY, JOB
321: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
322: * ..
323: * .. Array Arguments ..
324: LOGICAL SELECT( * )
325: INTEGER IWORK( * )
326: DOUBLE PRECISION DIF( * ), S( * )
327: COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
328: $ VR( LDVR, * ), WORK( * )
329: * ..
330: *
331: * =====================================================================
332: *
333: * .. Parameters ..
334: DOUBLE PRECISION ZERO, ONE
335: INTEGER IDIFJB
336: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
337: * ..
338: * .. Local Scalars ..
339: LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS
340: INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
341: DOUBLE PRECISION BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
342: COMPLEX*16 YHAX, YHBX
343: * ..
344: * .. Local Arrays ..
345: COMPLEX*16 DUMMY( 1 ), DUMMY1( 1 )
346: * ..
347: * .. External Functions ..
348: LOGICAL LSAME
349: DOUBLE PRECISION DLAMCH, DLAPY2, DZNRM2
350: COMPLEX*16 ZDOTC
351: EXTERNAL LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
352: * ..
353: * .. External Subroutines ..
354: EXTERNAL DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
355: * ..
356: * .. Intrinsic Functions ..
357: INTRINSIC ABS, DCMPLX, MAX
358: * ..
359: * .. Executable Statements ..
360: *
361: * Decode and test the input parameters
362: *
363: WANTBH = LSAME( JOB, 'B' )
364: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
365: WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
366: *
367: SOMCON = LSAME( HOWMNY, 'S' )
368: *
369: INFO = 0
370: LQUERY = ( LWORK.EQ.-1 )
371: *
372: IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
373: INFO = -1
374: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
375: INFO = -2
376: ELSE IF( N.LT.0 ) THEN
377: INFO = -4
378: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
379: INFO = -6
380: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
381: INFO = -8
382: ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
383: INFO = -10
384: ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
385: INFO = -12
386: ELSE
387: *
388: * Set M to the number of eigenpairs for which condition numbers
389: * are required, and test MM.
390: *
391: IF( SOMCON ) THEN
392: M = 0
393: DO 10 K = 1, N
394: IF( SELECT( K ) )
395: $ M = M + 1
396: 10 CONTINUE
397: ELSE
398: M = N
399: END IF
400: *
401: IF( N.EQ.0 ) THEN
402: LWMIN = 1
403: ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
404: LWMIN = 2*N*N
405: ELSE
406: LWMIN = N
407: END IF
408: WORK( 1 ) = LWMIN
409: *
410: IF( MM.LT.M ) THEN
411: INFO = -15
412: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
413: INFO = -18
414: END IF
415: END IF
416: *
417: IF( INFO.NE.0 ) THEN
418: CALL XERBLA( 'ZTGSNA', -INFO )
419: RETURN
420: ELSE IF( LQUERY ) THEN
421: RETURN
422: END IF
423: *
424: * Quick return if possible
425: *
426: IF( N.EQ.0 )
427: $ RETURN
428: *
429: * Get machine constants
430: *
431: EPS = DLAMCH( 'P' )
432: SMLNUM = DLAMCH( 'S' ) / EPS
433: BIGNUM = ONE / SMLNUM
434: CALL DLABAD( SMLNUM, BIGNUM )
435: KS = 0
436: DO 20 K = 1, N
437: *
438: * Determine whether condition numbers are required for the k-th
439: * eigenpair.
440: *
441: IF( SOMCON ) THEN
442: IF( .NOT.SELECT( K ) )
443: $ GO TO 20
444: END IF
445: *
446: KS = KS + 1
447: *
448: IF( WANTS ) THEN
449: *
450: * Compute the reciprocal condition number of the k-th
451: * eigenvalue.
452: *
453: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
454: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
455: CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
456: $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
457: YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
458: CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
459: $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
460: YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
461: COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
462: IF( COND.EQ.ZERO ) THEN
463: S( KS ) = -ONE
464: ELSE
465: S( KS ) = COND / ( RNRM*LNRM )
466: END IF
467: END IF
468: *
469: IF( WANTDF ) THEN
470: IF( N.EQ.1 ) THEN
471: DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
472: ELSE
473: *
474: * Estimate the reciprocal condition number of the k-th
475: * eigenvectors.
476: *
477: * Copy the matrix (A, B) to the array WORK and move the
478: * (k,k)th pair to the (1,1) position.
479: *
480: CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
481: CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
482: IFST = K
483: ILST = 1
484: *
485: CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
486: $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
487: *
488: IF( IERR.GT.0 ) THEN
489: *
490: * Ill-conditioned problem - swap rejected.
491: *
492: DIF( KS ) = ZERO
493: ELSE
494: *
495: * Reordering successful, solve generalized Sylvester
496: * equation for R and L,
497: * A22 * R - L * A11 = A12
498: * B22 * R - L * B11 = B12,
499: * and compute estimate of Difl[(A11,B11), (A22, B22)].
500: *
501: N1 = 1
502: N2 = N - N1
503: I = N*N + 1
504: CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
505: $ N, WORK, N, WORK( N1+1 ), N,
506: $ WORK( N*N1+N1+I ), N, WORK( I ), N,
507: $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
508: $ 1, IWORK, IERR )
509: END IF
510: END IF
511: END IF
512: *
513: 20 CONTINUE
514: WORK( 1 ) = LWMIN
515: RETURN
516: *
517: * End of ZTGSNA
518: *
519: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack