LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stgsna.f
Go to the documentation of this file.
00001 *> \brief \b STGSNA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STGSNA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsna.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsna.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsna.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
00022 *                          LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
00023 *                          IWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          HOWMNY, JOB
00027 *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       LOGICAL            SELECT( * )
00031 *       INTEGER            IWORK( * )
00032 *       REAL               A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
00033 *      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> STGSNA estimates reciprocal condition numbers for specified
00043 *> eigenvalues and/or eigenvectors of a matrix pair (A, B) in
00044 *> generalized real Schur canonical form (or of any matrix pair
00045 *> (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
00046 *> Z**T denotes the transpose of Z.
00047 *>
00048 *> (A, B) must be in generalized real Schur form (as returned by SGGES),
00049 *> i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
00050 *> blocks. B is upper triangular.
00051 *>
00052 *> \endverbatim
00053 *
00054 *  Arguments:
00055 *  ==========
00056 *
00057 *> \param[in] JOB
00058 *> \verbatim
00059 *>          JOB is CHARACTER*1
00060 *>          Specifies whether condition numbers are required for
00061 *>          eigenvalues (S) or eigenvectors (DIF):
00062 *>          = 'E': for eigenvalues only (S);
00063 *>          = 'V': for eigenvectors only (DIF);
00064 *>          = 'B': for both eigenvalues and eigenvectors (S and DIF).
00065 *> \endverbatim
00066 *>
00067 *> \param[in] HOWMNY
00068 *> \verbatim
00069 *>          HOWMNY is CHARACTER*1
00070 *>          = 'A': compute condition numbers for all eigenpairs;
00071 *>          = 'S': compute condition numbers for selected eigenpairs
00072 *>                 specified by the array SELECT.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] SELECT
00076 *> \verbatim
00077 *>          SELECT is LOGICAL array, dimension (N)
00078 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00079 *>          condition numbers are required. To select condition numbers
00080 *>          for the eigenpair corresponding to a real eigenvalue w(j),
00081 *>          SELECT(j) must be set to .TRUE.. To select condition numbers
00082 *>          corresponding to a complex conjugate pair of eigenvalues w(j)
00083 *>          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
00084 *>          set to .TRUE..
00085 *>          If HOWMNY = 'A', SELECT is not referenced.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] N
00089 *> \verbatim
00090 *>          N is INTEGER
00091 *>          The order of the square matrix pair (A, B). N >= 0.
00092 *> \endverbatim
00093 *>
00094 *> \param[in] A
00095 *> \verbatim
00096 *>          A is REAL array, dimension (LDA,N)
00097 *>          The upper quasi-triangular matrix A in the pair (A,B).
00098 *> \endverbatim
00099 *>
00100 *> \param[in] LDA
00101 *> \verbatim
00102 *>          LDA is INTEGER
00103 *>          The leading dimension of the array A. LDA >= max(1,N).
00104 *> \endverbatim
00105 *>
00106 *> \param[in] B
00107 *> \verbatim
00108 *>          B is REAL array, dimension (LDB,N)
00109 *>          The upper triangular matrix B in the pair (A,B).
00110 *> \endverbatim
00111 *>
00112 *> \param[in] LDB
00113 *> \verbatim
00114 *>          LDB is INTEGER
00115 *>          The leading dimension of the array B. LDB >= max(1,N).
00116 *> \endverbatim
00117 *>
00118 *> \param[in] VL
00119 *> \verbatim
00120 *>          VL is REAL array, dimension (LDVL,M)
00121 *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of
00122 *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
00123 *>          and SELECT. The eigenvectors must be stored in consecutive
00124 *>          columns of VL, as returned by STGEVC.
00125 *>          If JOB = 'V', VL is not referenced.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] LDVL
00129 *> \verbatim
00130 *>          LDVL is INTEGER
00131 *>          The leading dimension of the array VL. LDVL >= 1.
00132 *>          If JOB = 'E' or 'B', LDVL >= N.
00133 *> \endverbatim
00134 *>
00135 *> \param[in] VR
00136 *> \verbatim
00137 *>          VR is REAL array, dimension (LDVR,M)
00138 *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of
00139 *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
00140 *>          and SELECT. The eigenvectors must be stored in consecutive
00141 *>          columns ov VR, as returned by STGEVC.
00142 *>          If JOB = 'V', VR is not referenced.
00143 *> \endverbatim
00144 *>
00145 *> \param[in] LDVR
00146 *> \verbatim
00147 *>          LDVR is INTEGER
00148 *>          The leading dimension of the array VR. LDVR >= 1.
00149 *>          If JOB = 'E' or 'B', LDVR >= N.
00150 *> \endverbatim
00151 *>
00152 *> \param[out] S
00153 *> \verbatim
00154 *>          S is REAL array, dimension (MM)
00155 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00156 *>          selected eigenvalues, stored in consecutive elements of the
00157 *>          array. For a complex conjugate pair of eigenvalues two
00158 *>          consecutive elements of S are set to the same value. Thus
00159 *>          S(j), DIF(j), and the j-th columns of VL and VR all
00160 *>          correspond to the same eigenpair (but not in general the
00161 *>          j-th eigenpair, unless all eigenpairs are selected).
00162 *>          If JOB = 'V', S is not referenced.
00163 *> \endverbatim
00164 *>
00165 *> \param[out] DIF
00166 *> \verbatim
00167 *>          DIF is REAL array, dimension (MM)
00168 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
00169 *>          numbers of the selected eigenvectors, stored in consecutive
00170 *>          elements of the array. For a complex eigenvector two
00171 *>          consecutive elements of DIF are set to the same value. If
00172 *>          the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
00173 *>          is set to 0; this can only occur when the true value would be
00174 *>          very small anyway.
00175 *>          If JOB = 'E', DIF is not referenced.
00176 *> \endverbatim
00177 *>
00178 *> \param[in] MM
00179 *> \verbatim
00180 *>          MM is INTEGER
00181 *>          The number of elements in the arrays S and DIF. MM >= M.
00182 *> \endverbatim
00183 *>
00184 *> \param[out] M
00185 *> \verbatim
00186 *>          M is INTEGER
00187 *>          The number of elements of the arrays S and DIF used to store
00188 *>          the specified condition numbers; for each selected real
00189 *>          eigenvalue one element is used, and for each selected complex
00190 *>          conjugate pair of eigenvalues, two elements are used.
00191 *>          If HOWMNY = 'A', M is set to N.
00192 *> \endverbatim
00193 *>
00194 *> \param[out] WORK
00195 *> \verbatim
00196 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00197 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00198 *> \endverbatim
00199 *>
00200 *> \param[in] LWORK
00201 *> \verbatim
00202 *>          LWORK is INTEGER
00203 *>          The dimension of the array WORK. LWORK >= max(1,N).
00204 *>          If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
00205 *>
00206 *>          If LWORK = -1, then a workspace query is assumed; the routine
00207 *>          only calculates the optimal size of the WORK array, returns
00208 *>          this value as the first entry of the WORK array, and no error
00209 *>          message related to LWORK is issued by XERBLA.
00210 *> \endverbatim
00211 *>
00212 *> \param[out] IWORK
00213 *> \verbatim
00214 *>          IWORK is INTEGER array, dimension (N + 6)
00215 *>          If JOB = 'E', IWORK is not referenced.
00216 *> \endverbatim
00217 *>
00218 *> \param[out] INFO
00219 *> \verbatim
00220 *>          INFO is INTEGER
00221 *>          =0: Successful exit
00222 *>          <0: If INFO = -i, the i-th argument had an illegal value
00223 *> \endverbatim
00224 *
00225 *  Authors:
00226 *  ========
00227 *
00228 *> \author Univ. of Tennessee 
00229 *> \author Univ. of California Berkeley 
00230 *> \author Univ. of Colorado Denver 
00231 *> \author NAG Ltd. 
00232 *
00233 *> \date November 2011
00234 *
00235 *> \ingroup realOTHERcomputational
00236 *
00237 *> \par Further Details:
00238 *  =====================
00239 *>
00240 *> \verbatim
00241 *>
00242 *>  The reciprocal of the condition number of a generalized eigenvalue
00243 *>  w = (a, b) is defined as
00244 *>
00245 *>       S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
00246 *>
00247 *>  where u and v are the left and right eigenvectors of (A, B)
00248 *>  corresponding to w; |z| denotes the absolute value of the complex
00249 *>  number, and norm(u) denotes the 2-norm of the vector u.
00250 *>  The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
00251 *>  of the matrix pair (A, B). If both a and b equal zero, then (A B) is
00252 *>  singular and S(I) = -1 is returned.
00253 *>
00254 *>  An approximate error bound on the chordal distance between the i-th
00255 *>  computed generalized eigenvalue w and the corresponding exact
00256 *>  eigenvalue lambda is
00257 *>
00258 *>       chord(w, lambda) <= EPS * norm(A, B) / S(I)
00259 *>
00260 *>  where EPS is the machine precision.
00261 *>
00262 *>  The reciprocal of the condition number DIF(i) of right eigenvector u
00263 *>  and left eigenvector v corresponding to the generalized eigenvalue w
00264 *>  is defined as follows:
00265 *>
00266 *>  a) If the i-th eigenvalue w = (a,b) is real
00267 *>
00268 *>     Suppose U and V are orthogonal transformations such that
00269 *>
00270 *>              U**T*(A, B)*V  = (S, T) = ( a   *  ) ( b  *  )  1
00271 *>                                        ( 0  S22 ),( 0 T22 )  n-1
00272 *>                                          1  n-1     1 n-1
00273 *>
00274 *>     Then the reciprocal condition number DIF(i) is
00275 *>
00276 *>                Difl((a, b), (S22, T22)) = sigma-min( Zl ),
00277 *>
00278 *>     where sigma-min(Zl) denotes the smallest singular value of the
00279 *>     2(n-1)-by-2(n-1) matrix
00280 *>
00281 *>         Zl = [ kron(a, In-1)  -kron(1, S22) ]
00282 *>              [ kron(b, In-1)  -kron(1, T22) ] .
00283 *>
00284 *>     Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
00285 *>     Kronecker product between the matrices X and Y.
00286 *>
00287 *>     Note that if the default method for computing DIF(i) is wanted
00288 *>     (see SLATDF), then the parameter DIFDRI (see below) should be
00289 *>     changed from 3 to 4 (routine SLATDF(IJOB = 2 will be used)).
00290 *>     See STGSYL for more details.
00291 *>
00292 *>  b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
00293 *>
00294 *>     Suppose U and V are orthogonal transformations such that
00295 *>
00296 *>              U**T*(A, B)*V = (S, T) = ( S11  *   ) ( T11  *  )  2
00297 *>                                       ( 0    S22 ),( 0    T22) n-2
00298 *>                                         2    n-2     2    n-2
00299 *>
00300 *>     and (S11, T11) corresponds to the complex conjugate eigenvalue
00301 *>     pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
00302 *>     that
00303 *>
00304 *>       U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
00305 *>                      (  0  s22 )                    (  0  t22 )
00306 *>
00307 *>     where the generalized eigenvalues w = s11/t11 and
00308 *>     conjg(w) = s22/t22.
00309 *>
00310 *>     Then the reciprocal condition number DIF(i) is bounded by
00311 *>
00312 *>         min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
00313 *>
00314 *>     where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
00315 *>     Z1 is the complex 2-by-2 matrix
00316 *>
00317 *>              Z1 =  [ s11  -s22 ]
00318 *>                    [ t11  -t22 ],
00319 *>
00320 *>     This is done by computing (using real arithmetic) the
00321 *>     roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
00322 *>     where Z1**T denotes the transpose of Z1 and det(X) denotes
00323 *>     the determinant of X.
00324 *>
00325 *>     and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
00326 *>     upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
00327 *>
00328 *>              Z2 = [ kron(S11**T, In-2)  -kron(I2, S22) ]
00329 *>                   [ kron(T11**T, In-2)  -kron(I2, T22) ]
00330 *>
00331 *>     Note that if the default method for computing DIF is wanted (see
00332 *>     SLATDF), then the parameter DIFDRI (see below) should be changed
00333 *>     from 3 to 4 (routine SLATDF(IJOB = 2 will be used)). See STGSYL
00334 *>     for more details.
00335 *>
00336 *>  For each eigenvalue/vector specified by SELECT, DIF stores a
00337 *>  Frobenius norm-based estimate of Difl.
00338 *>
00339 *>  An approximate error bound for the i-th computed eigenvector VL(i) or
00340 *>  VR(i) is given by
00341 *>
00342 *>             EPS * norm(A, B) / DIF(i).
00343 *>
00344 *>  See ref. [2-3] for more details and further references.
00345 *> \endverbatim
00346 *
00347 *> \par Contributors:
00348 *  ==================
00349 *>
00350 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00351 *>     Umea University, S-901 87 Umea, Sweden.
00352 *
00353 *> \par References:
00354 *  ================
00355 *>
00356 *> \verbatim
00357 *>
00358 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00359 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00360 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00361 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00362 *>
00363 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
00364 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
00365 *>      Estimation: Theory, Algorithms and Software,
00366 *>      Report UMINF - 94.04, Department of Computing Science, Umea
00367 *>      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
00368 *>      Note 87. To appear in Numerical Algorithms, 1996.
00369 *>
00370 *>  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
00371 *>      for Solving the Generalized Sylvester Equation and Estimating the
00372 *>      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
00373 *>      Department of Computing Science, Umea University, S-901 87 Umea,
00374 *>      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
00375 *>      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
00376 *>      No 1, 1996.
00377 *> \endverbatim
00378 *>
00379 *  =====================================================================
00380       SUBROUTINE STGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
00381      $                   LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
00382      $                   IWORK, INFO )
00383 *
00384 *  -- LAPACK computational routine (version 3.4.0) --
00385 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00386 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00387 *     November 2011
00388 *
00389 *     .. Scalar Arguments ..
00390       CHARACTER          HOWMNY, JOB
00391       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
00392 *     ..
00393 *     .. Array Arguments ..
00394       LOGICAL            SELECT( * )
00395       INTEGER            IWORK( * )
00396       REAL               A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
00397      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
00398 *     ..
00399 *
00400 *  =====================================================================
00401 *
00402 *     .. Parameters ..
00403       INTEGER            DIFDRI
00404       PARAMETER          ( DIFDRI = 3 )
00405       REAL               ZERO, ONE, TWO, FOUR
00406       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
00407      $                   FOUR = 4.0E+0 )
00408 *     ..
00409 *     .. Local Scalars ..
00410       LOGICAL            LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
00411       INTEGER            I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
00412       REAL               ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
00413      $                   EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM,
00414      $                   TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV,
00415      $                   UHBVI
00416 *     ..
00417 *     .. Local Arrays ..
00418       REAL               DUMMY( 1 ), DUMMY1( 1 )
00419 *     ..
00420 *     .. External Functions ..
00421       LOGICAL            LSAME
00422       REAL               SDOT, SLAMCH, SLAPY2, SNRM2
00423       EXTERNAL           LSAME, SDOT, SLAMCH, SLAPY2, SNRM2
00424 *     ..
00425 *     .. External Subroutines ..
00426       EXTERNAL           SGEMV, SLACPY, SLAG2, STGEXC, STGSYL, XERBLA
00427 *     ..
00428 *     .. Intrinsic Functions ..
00429       INTRINSIC          MAX, MIN, SQRT
00430 *     ..
00431 *     .. Executable Statements ..
00432 *
00433 *     Decode and test the input parameters
00434 *
00435       WANTBH = LSAME( JOB, 'B' )
00436       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00437       WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
00438 *
00439       SOMCON = LSAME( HOWMNY, 'S' )
00440 *
00441       INFO = 0
00442       LQUERY = ( LWORK.EQ.-1 )
00443 *
00444       IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
00445          INFO = -1
00446       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00447          INFO = -2
00448       ELSE IF( N.LT.0 ) THEN
00449          INFO = -4
00450       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00451          INFO = -6
00452       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00453          INFO = -8
00454       ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
00455          INFO = -10
00456       ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
00457          INFO = -12
00458       ELSE
00459 *
00460 *        Set M to the number of eigenpairs for which condition numbers
00461 *        are required, and test MM.
00462 *
00463          IF( SOMCON ) THEN
00464             M = 0
00465             PAIR = .FALSE.
00466             DO 10 K = 1, N
00467                IF( PAIR ) THEN
00468                   PAIR = .FALSE.
00469                ELSE
00470                   IF( K.LT.N ) THEN
00471                      IF( A( K+1, K ).EQ.ZERO ) THEN
00472                         IF( SELECT( K ) )
00473      $                     M = M + 1
00474                      ELSE
00475                         PAIR = .TRUE.
00476                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
00477      $                     M = M + 2
00478                      END IF
00479                   ELSE
00480                      IF( SELECT( N ) )
00481      $                  M = M + 1
00482                   END IF
00483                END IF
00484    10       CONTINUE
00485          ELSE
00486             M = N
00487          END IF
00488 *
00489          IF( N.EQ.0 ) THEN
00490             LWMIN = 1
00491          ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
00492             LWMIN = 2*N*( N + 2 ) + 16
00493          ELSE
00494             LWMIN = N
00495          END IF
00496          WORK( 1 ) = LWMIN
00497 *
00498          IF( MM.LT.M ) THEN
00499             INFO = -15
00500          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00501             INFO = -18
00502          END IF
00503       END IF
00504 *
00505       IF( INFO.NE.0 ) THEN
00506          CALL XERBLA( 'STGSNA', -INFO )
00507          RETURN
00508       ELSE IF( LQUERY ) THEN
00509          RETURN
00510       END IF
00511 *
00512 *     Quick return if possible
00513 *
00514       IF( N.EQ.0 )
00515      $   RETURN
00516 *
00517 *     Get machine constants
00518 *
00519       EPS = SLAMCH( 'P' )
00520       SMLNUM = SLAMCH( 'S' ) / EPS
00521       KS = 0
00522       PAIR = .FALSE.
00523 *
00524       DO 20 K = 1, N
00525 *
00526 *        Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block.
00527 *
00528          IF( PAIR ) THEN
00529             PAIR = .FALSE.
00530             GO TO 20
00531          ELSE
00532             IF( K.LT.N )
00533      $         PAIR = A( K+1, K ).NE.ZERO
00534          END IF
00535 *
00536 *        Determine whether condition numbers are required for the k-th
00537 *        eigenpair.
00538 *
00539          IF( SOMCON ) THEN
00540             IF( PAIR ) THEN
00541                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
00542      $            GO TO 20
00543             ELSE
00544                IF( .NOT.SELECT( K ) )
00545      $            GO TO 20
00546             END IF
00547          END IF
00548 *
00549          KS = KS + 1
00550 *
00551          IF( WANTS ) THEN
00552 *
00553 *           Compute the reciprocal condition number of the k-th
00554 *           eigenvalue.
00555 *
00556             IF( PAIR ) THEN
00557 *
00558 *              Complex eigenvalue pair.
00559 *
00560                RNRM = SLAPY2( SNRM2( N, VR( 1, KS ), 1 ),
00561      $                SNRM2( N, VR( 1, KS+1 ), 1 ) )
00562                LNRM = SLAPY2( SNRM2( N, VL( 1, KS ), 1 ),
00563      $                SNRM2( N, VL( 1, KS+1 ), 1 ) )
00564                CALL SGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
00565      $                     WORK, 1 )
00566                TMPRR = SDOT( N, WORK, 1, VL( 1, KS ), 1 )
00567                TMPRI = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
00568                CALL SGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1,
00569      $                     ZERO, WORK, 1 )
00570                TMPII = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
00571                TMPIR = SDOT( N, WORK, 1, VL( 1, KS ), 1 )
00572                UHAV = TMPRR + TMPII
00573                UHAVI = TMPIR - TMPRI
00574                CALL SGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
00575      $                     WORK, 1 )
00576                TMPRR = SDOT( N, WORK, 1, VL( 1, KS ), 1 )
00577                TMPRI = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
00578                CALL SGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1,
00579      $                     ZERO, WORK, 1 )
00580                TMPII = SDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
00581                TMPIR = SDOT( N, WORK, 1, VL( 1, KS ), 1 )
00582                UHBV = TMPRR + TMPII
00583                UHBVI = TMPIR - TMPRI
00584                UHAV = SLAPY2( UHAV, UHAVI )
00585                UHBV = SLAPY2( UHBV, UHBVI )
00586                COND = SLAPY2( UHAV, UHBV )
00587                S( KS ) = COND / ( RNRM*LNRM )
00588                S( KS+1 ) = S( KS )
00589 *
00590             ELSE
00591 *
00592 *              Real eigenvalue.
00593 *
00594                RNRM = SNRM2( N, VR( 1, KS ), 1 )
00595                LNRM = SNRM2( N, VL( 1, KS ), 1 )
00596                CALL SGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
00597      $                     WORK, 1 )
00598                UHAV = SDOT( N, WORK, 1, VL( 1, KS ), 1 )
00599                CALL SGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
00600      $                     WORK, 1 )
00601                UHBV = SDOT( N, WORK, 1, VL( 1, KS ), 1 )
00602                COND = SLAPY2( UHAV, UHBV )
00603                IF( COND.EQ.ZERO ) THEN
00604                   S( KS ) = -ONE
00605                ELSE
00606                   S( KS ) = COND / ( RNRM*LNRM )
00607                END IF
00608             END IF
00609          END IF
00610 *
00611          IF( WANTDF ) THEN
00612             IF( N.EQ.1 ) THEN
00613                DIF( KS ) = SLAPY2( A( 1, 1 ), B( 1, 1 ) )
00614                GO TO 20
00615             END IF
00616 *
00617 *           Estimate the reciprocal condition number of the k-th
00618 *           eigenvectors.
00619             IF( PAIR ) THEN
00620 *
00621 *              Copy the  2-by 2 pencil beginning at (A(k,k), B(k, k)).
00622 *              Compute the eigenvalue(s) at position K.
00623 *
00624                WORK( 1 ) = A( K, K )
00625                WORK( 2 ) = A( K+1, K )
00626                WORK( 3 ) = A( K, K+1 )
00627                WORK( 4 ) = A( K+1, K+1 )
00628                WORK( 5 ) = B( K, K )
00629                WORK( 6 ) = B( K+1, K )
00630                WORK( 7 ) = B( K, K+1 )
00631                WORK( 8 ) = B( K+1, K+1 )
00632                CALL SLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA,
00633      $                     DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI )
00634                ALPRQT = ONE
00635                C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA )
00636                C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI
00637                ROOT1 = C1 + SQRT( C1*C1-4.0*C2 )
00638                ROOT2 = C2 / ROOT1
00639                ROOT1 = ROOT1 / TWO
00640                COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) )
00641             END IF
00642 *
00643 *           Copy the matrix (A, B) to the array WORK and swap the
00644 *           diagonal block beginning at A(k,k) to the (1,1) position.
00645 *
00646             CALL SLACPY( 'Full', N, N, A, LDA, WORK, N )
00647             CALL SLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
00648             IFST = K
00649             ILST = 1
00650 *
00651             CALL STGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N,
00652      $                   DUMMY, 1, DUMMY1, 1, IFST, ILST,
00653      $                   WORK( N*N*2+1 ), LWORK-2*N*N, IERR )
00654 *
00655             IF( IERR.GT.0 ) THEN
00656 *
00657 *              Ill-conditioned problem - swap rejected.
00658 *
00659                DIF( KS ) = ZERO
00660             ELSE
00661 *
00662 *              Reordering successful, solve generalized Sylvester
00663 *              equation for R and L,
00664 *                         A22 * R - L * A11 = A12
00665 *                         B22 * R - L * B11 = B12,
00666 *              and compute estimate of Difl((A11,B11), (A22, B22)).
00667 *
00668                N1 = 1
00669                IF( WORK( 2 ).NE.ZERO )
00670      $            N1 = 2
00671                N2 = N - N1
00672                IF( N2.EQ.0 ) THEN
00673                   DIF( KS ) = COND
00674                ELSE
00675                   I = N*N + 1
00676                   IZ = 2*N*N + 1
00677                   CALL STGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ),
00678      $                         N, WORK, N, WORK( N1+1 ), N,
00679      $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
00680      $                         WORK( N1+I ), N, SCALE, DIF( KS ),
00681      $                         WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR )
00682 *
00683                   IF( PAIR )
00684      $               DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ),
00685      $                           COND )
00686                END IF
00687             END IF
00688             IF( PAIR )
00689      $         DIF( KS+1 ) = DIF( KS )
00690          END IF
00691          IF( PAIR )
00692      $      KS = KS + 1
00693 *
00694    20 CONTINUE
00695       WORK( 1 ) = LWMIN
00696       RETURN
00697 *
00698 *     End of STGSNA
00699 *
00700       END
 All Files Functions