LAPACK 3.3.0

dtgsna.f

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