SRC/ztgsna.f

   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 *  =====================================================================
 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 *
 314 *  -- LAPACK computational routine (version 3.4.0) --
 315 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 316 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 317 *     November 2011
 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