1*> \brief \b CTGSNA 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CTGSNA + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsna.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsna.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsna.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CTGSNA( 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* REAL DIF( * ), S( * ) 33* COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 34* $ VR( LDVR, * ), WORK( * ) 35* .. 36* 37* 38*> \par Purpose: 39* ============= 40*> 41*> \verbatim 42*> 43*> CTGSNA 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 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 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 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 CTGEVC. 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 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 CTGEVC. 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 REAL 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 REAL 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 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*> \ingroup complexOTHERcomputational 217* 218*> \par Further Details: 219* ===================== 220*> 221*> \verbatim 222*> 223*> The reciprocal of the condition number of the i-th generalized 224*> eigenvalue w = (a, b) is defined as 225*> 226*> S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v)) 227*> 228*> where u and v are the right and left eigenvectors of (A, B) 229*> corresponding to w; |z| denotes the absolute value of the complex 230*> number, and norm(u) denotes the 2-norm of the vector u. The pair 231*> (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the 232*> matrix pair (A, B). If both a and b equal zero, then (A,B) is 233*> singular and S(I) = -1 is returned. 234*> 235*> An approximate error bound on the chordal distance between the i-th 236*> computed generalized eigenvalue w and the corresponding exact 237*> eigenvalue lambda is 238*> 239*> chord(w, lambda) <= EPS * norm(A, B) / S(I), 240*> 241*> where EPS is the machine precision. 242*> 243*> The reciprocal of the condition number of the right eigenvector u 244*> and left eigenvector v corresponding to the generalized eigenvalue w 245*> is defined as follows. Suppose 246*> 247*> (A, B) = ( a * ) ( b * ) 1 248*> ( 0 A22 ),( 0 B22 ) n-1 249*> 1 n-1 1 n-1 250*> 251*> Then the reciprocal condition number DIF(I) is 252*> 253*> Difl[(a, b), (A22, B22)] = sigma-min( Zl ) 254*> 255*> where sigma-min(Zl) denotes the smallest singular value of 256*> 257*> Zl = [ kron(a, In-1) -kron(1, A22) ] 258*> [ kron(b, In-1) -kron(1, B22) ]. 259*> 260*> Here In-1 is the identity matrix of size n-1 and X**H is the conjugate 261*> transpose of X. kron(X, Y) is the Kronecker product between the 262*> matrices X and Y. 263*> 264*> We approximate the smallest singular value of Zl with an upper 265*> bound. This is done by CLATDF. 266*> 267*> An approximate error bound for a computed eigenvector VL(i) or 268*> VR(i) is given by 269*> 270*> EPS * norm(A, B) / DIF(i). 271*> 272*> See ref. [2-3] for more details and further references. 273*> \endverbatim 274* 275*> \par Contributors: 276* ================== 277*> 278*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 279*> Umea University, S-901 87 Umea, Sweden. 280* 281*> \par References: 282* ================ 283*> 284*> \verbatim 285*> 286*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 287*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 288*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 289*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 290*> 291*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 292*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 293*> Estimation: Theory, Algorithms and Software, Report 294*> UMINF - 94.04, Department of Computing Science, Umea University, 295*> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. 296*> To appear in Numerical Algorithms, 1996. 297*> 298*> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 299*> for Solving the Generalized Sylvester Equation and Estimating the 300*> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 301*> Department of Computing Science, Umea University, S-901 87 Umea, 302*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 303*> Note 75. 304*> To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996. 305*> \endverbatim 306*> 307* ===================================================================== 308 SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 309 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 310 $ IWORK, INFO ) 311* 312* -- LAPACK computational routine -- 313* -- LAPACK is a software package provided by Univ. of Tennessee, -- 314* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 315* 316* .. Scalar Arguments .. 317 CHARACTER HOWMNY, JOB 318 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 319* .. 320* .. Array Arguments .. 321 LOGICAL SELECT( * ) 322 INTEGER IWORK( * ) 323 REAL DIF( * ), S( * ) 324 COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 325 $ VR( LDVR, * ), WORK( * ) 326* .. 327* 328* ===================================================================== 329* 330* .. Parameters .. 331 REAL ZERO, ONE 332 INTEGER IDIFJB 333 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, IDIFJB = 3 ) 334* .. 335* .. Local Scalars .. 336 LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS 337 INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2 338 REAL BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM 339 COMPLEX YHAX, YHBX 340* .. 341* .. Local Arrays .. 342 COMPLEX DUMMY( 1 ), DUMMY1( 1 ) 343* .. 344* .. External Functions .. 345 LOGICAL LSAME 346 REAL SCNRM2, SLAMCH, SLAPY2 347 COMPLEX CDOTC 348 EXTERNAL LSAME, SCNRM2, SLAMCH, SLAPY2, CDOTC 349* .. 350* .. External Subroutines .. 351 EXTERNAL CGEMV, CLACPY, CTGEXC, CTGSYL, SLABAD, XERBLA 352* .. 353* .. Intrinsic Functions .. 354 INTRINSIC ABS, CMPLX, MAX 355* .. 356* .. Executable Statements .. 357* 358* Decode and test the input parameters 359* 360 WANTBH = LSAME( JOB, 'B' ) 361 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 362 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH 363* 364 SOMCON = LSAME( HOWMNY, 'S' ) 365* 366 INFO = 0 367 LQUERY = ( LWORK.EQ.-1 ) 368* 369 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN 370 INFO = -1 371 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 372 INFO = -2 373 ELSE IF( N.LT.0 ) THEN 374 INFO = -4 375 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 376 INFO = -6 377 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 378 INFO = -8 379 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN 380 INFO = -10 381 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN 382 INFO = -12 383 ELSE 384* 385* Set M to the number of eigenpairs for which condition numbers 386* are required, and test MM. 387* 388 IF( SOMCON ) THEN 389 M = 0 390 DO 10 K = 1, N 391 IF( SELECT( K ) ) 392 $ M = M + 1 393 10 CONTINUE 394 ELSE 395 M = N 396 END IF 397* 398 IF( N.EQ.0 ) THEN 399 LWMIN = 1 400 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN 401 LWMIN = 2*N*N 402 ELSE 403 LWMIN = N 404 END IF 405 WORK( 1 ) = LWMIN 406* 407 IF( MM.LT.M ) THEN 408 INFO = -15 409 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 410 INFO = -18 411 END IF 412 END IF 413* 414 IF( INFO.NE.0 ) THEN 415 CALL XERBLA( 'CTGSNA', -INFO ) 416 RETURN 417 ELSE IF( LQUERY ) THEN 418 RETURN 419 END IF 420* 421* Quick return if possible 422* 423 IF( N.EQ.0 ) 424 $ RETURN 425* 426* Get machine constants 427* 428 EPS = SLAMCH( 'P' ) 429 SMLNUM = SLAMCH( 'S' ) / EPS 430 BIGNUM = ONE / SMLNUM 431 CALL SLABAD( SMLNUM, BIGNUM ) 432 KS = 0 433 DO 20 K = 1, N 434* 435* Determine whether condition numbers are required for the k-th 436* eigenpair. 437* 438 IF( SOMCON ) THEN 439 IF( .NOT.SELECT( K ) ) 440 $ GO TO 20 441 END IF 442* 443 KS = KS + 1 444* 445 IF( WANTS ) THEN 446* 447* Compute the reciprocal condition number of the k-th 448* eigenvalue. 449* 450 RNRM = SCNRM2( N, VR( 1, KS ), 1 ) 451 LNRM = SCNRM2( N, VL( 1, KS ), 1 ) 452 CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), A, LDA, 453 $ VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 ) 454 YHAX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 ) 455 CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), B, LDB, 456 $ VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 ) 457 YHBX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 ) 458 COND = SLAPY2( ABS( YHAX ), ABS( YHBX ) ) 459 IF( COND.EQ.ZERO ) THEN 460 S( KS ) = -ONE 461 ELSE 462 S( KS ) = COND / ( RNRM*LNRM ) 463 END IF 464 END IF 465* 466 IF( WANTDF ) THEN 467 IF( N.EQ.1 ) THEN 468 DIF( KS ) = SLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) ) 469 ELSE 470* 471* Estimate the reciprocal condition number of the k-th 472* eigenvectors. 473* 474* Copy the matrix (A, B) to the array WORK and move the 475* (k,k)th pair to the (1,1) position. 476* 477 CALL CLACPY( 'Full', N, N, A, LDA, WORK, N ) 478 CALL CLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 479 IFST = K 480 ILST = 1 481* 482 CALL CTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), 483 $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR ) 484* 485 IF( IERR.GT.0 ) THEN 486* 487* Ill-conditioned problem - swap rejected. 488* 489 DIF( KS ) = ZERO 490 ELSE 491* 492* Reordering successful, solve generalized Sylvester 493* equation for R and L, 494* A22 * R - L * A11 = A12 495* B22 * R - L * B11 = B12, 496* and compute estimate of Difl[(A11,B11), (A22, B22)]. 497* 498 N1 = 1 499 N2 = N - N1 500 I = N*N + 1 501 CALL CTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ), 502 $ N, WORK, N, WORK( N1+1 ), N, 503 $ WORK( N*N1+N1+I ), N, WORK( I ), N, 504 $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY, 505 $ 1, IWORK, IERR ) 506 END IF 507 END IF 508 END IF 509* 510 20 CONTINUE 511 WORK( 1 ) = LWMIN 512 RETURN 513* 514* End of CTGSNA 515* 516 END 517