--- rpl/lapack/lapack/dtgsna.f 2011/07/22 07:38:12 1.8
+++ rpl/lapack/lapack/dtgsna.f 2011/11/21 20:43:06 1.9
@@ -1,11 +1,390 @@
+*> \brief \b DTGSNA
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DTGSNA + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
+* LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
+* IWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER HOWMNY, JOB
+* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+* LOGICAL SELECT( * )
+* INTEGER IWORK( * )
+* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
+* $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DTGSNA estimates reciprocal condition numbers for specified
+*> eigenvalues and/or eigenvectors of a matrix pair (A, B) in
+*> generalized real Schur canonical form (or of any matrix pair
+*> (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
+*> Z**T denotes the transpose of Z.
+*>
+*> (A, B) must be in generalized real Schur form (as returned by DGGES),
+*> i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
+*> blocks. B is upper triangular.
+*>
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] JOB
+*> \verbatim
+*> JOB is CHARACTER*1
+*> Specifies whether condition numbers are required for
+*> eigenvalues (S) or eigenvectors (DIF):
+*> = 'E': for eigenvalues only (S);
+*> = 'V': for eigenvectors only (DIF);
+*> = 'B': for both eigenvalues and eigenvectors (S and DIF).
+*> \endverbatim
+*>
+*> \param[in] HOWMNY
+*> \verbatim
+*> HOWMNY is CHARACTER*1
+*> = 'A': compute condition numbers for all eigenpairs;
+*> = 'S': compute condition numbers for selected eigenpairs
+*> specified by the array SELECT.
+*> \endverbatim
+*>
+*> \param[in] SELECT
+*> \verbatim
+*> SELECT is LOGICAL array, dimension (N)
+*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
+*> condition numbers are required. To select condition numbers
+*> for the eigenpair corresponding to a real eigenvalue w(j),
+*> SELECT(j) must be set to .TRUE.. To select condition numbers
+*> corresponding to a complex conjugate pair of eigenvalues w(j)
+*> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
+*> set to .TRUE..
+*> If HOWMNY = 'A', SELECT is not referenced.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the square matrix pair (A, B). N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is DOUBLE PRECISION array, dimension (LDA,N)
+*> The upper quasi-triangular matrix A in the pair (A,B).
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] B
+*> \verbatim
+*> B is DOUBLE PRECISION array, dimension (LDB,N)
+*> The upper triangular matrix B in the pair (A,B).
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] VL
+*> \verbatim
+*> VL is DOUBLE PRECISION array, dimension (LDVL,M)
+*> If JOB = 'E' or 'B', VL must contain left eigenvectors of
+*> (A, B), corresponding to the eigenpairs specified by HOWMNY
+*> and SELECT. The eigenvectors must be stored in consecutive
+*> columns of VL, as returned by DTGEVC.
+*> If JOB = 'V', VL is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*> LDVL is INTEGER
+*> The leading dimension of the array VL. LDVL >= 1.
+*> If JOB = 'E' or 'B', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[in] VR
+*> \verbatim
+*> VR is DOUBLE PRECISION array, dimension (LDVR,M)
+*> If JOB = 'E' or 'B', VR must contain right eigenvectors of
+*> (A, B), corresponding to the eigenpairs specified by HOWMNY
+*> and SELECT. The eigenvectors must be stored in consecutive
+*> columns ov VR, as returned by DTGEVC.
+*> If JOB = 'V', VR is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*> LDVR is INTEGER
+*> The leading dimension of the array VR. LDVR >= 1.
+*> If JOB = 'E' or 'B', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[out] S
+*> \verbatim
+*> S is DOUBLE PRECISION array, dimension (MM)
+*> If JOB = 'E' or 'B', the reciprocal condition numbers of the
+*> selected eigenvalues, stored in consecutive elements of the
+*> array. For a complex conjugate pair of eigenvalues two
+*> consecutive elements of S are set to the same value. Thus
+*> S(j), DIF(j), and the j-th columns of VL and VR all
+*> correspond to the same eigenpair (but not in general the
+*> j-th eigenpair, unless all eigenpairs are selected).
+*> If JOB = 'V', S is not referenced.
+*> \endverbatim
+*>
+*> \param[out] DIF
+*> \verbatim
+*> DIF is DOUBLE PRECISION array, dimension (MM)
+*> If JOB = 'V' or 'B', the estimated reciprocal condition
+*> numbers of the selected eigenvectors, stored in consecutive
+*> elements of the array. For a complex eigenvector two
+*> consecutive elements of DIF are set to the same value. If
+*> the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
+*> is set to 0; this can only occur when the true value would be
+*> very small anyway.
+*> If JOB = 'E', DIF is not referenced.
+*> \endverbatim
+*>
+*> \param[in] MM
+*> \verbatim
+*> MM is INTEGER
+*> The number of elements in the arrays S and DIF. MM >= M.
+*> \endverbatim
+*>
+*> \param[out] M
+*> \verbatim
+*> M is INTEGER
+*> The number of elements of the arrays S and DIF used to store
+*> the specified condition numbers; for each selected real
+*> eigenvalue one element is used, and for each selected complex
+*> conjugate pair of eigenvalues, two elements are used.
+*> If HOWMNY = 'A', M is set to N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK. LWORK >= max(1,N).
+*> If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
+*>
+*> If LWORK = -1, then a workspace query is assumed; the routine
+*> only calculates the optimal size of the WORK array, returns
+*> this value as the first entry of the WORK array, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*> IWORK is INTEGER array, dimension (N + 6)
+*> If JOB = 'E', IWORK is not referenced.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> =0: Successful exit
+*> <0: If INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup doubleOTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> The reciprocal of the condition number of a generalized eigenvalue
+*> w = (a, b) is defined as
+*>
+*> S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
+*>
+*> where u and v are the left and right eigenvectors of (A, B)
+*> corresponding to w; |z| denotes the absolute value of the complex
+*> number, and norm(u) denotes the 2-norm of the vector u.
+*> The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
+*> of the matrix pair (A, B). If both a and b equal zero, then (A B) is
+*> singular and S(I) = -1 is returned.
+*>
+*> An approximate error bound on the chordal distance between the i-th
+*> computed generalized eigenvalue w and the corresponding exact
+*> eigenvalue lambda is
+*>
+*> chord(w, lambda) <= EPS * norm(A, B) / S(I)
+*>
+*> where EPS is the machine precision.
+*>
+*> The reciprocal of the condition number DIF(i) of right eigenvector u
+*> and left eigenvector v corresponding to the generalized eigenvalue w
+*> is defined as follows:
+*>
+*> a) If the i-th eigenvalue w = (a,b) is real
+*>
+*> Suppose U and V are orthogonal transformations such that
+*>
+*> U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
+*> ( 0 S22 ),( 0 T22 ) n-1
+*> 1 n-1 1 n-1
+*>
+*> Then the reciprocal condition number DIF(i) is
+*>
+*> Difl((a, b), (S22, T22)) = sigma-min( Zl ),
+*>
+*> where sigma-min(Zl) denotes the smallest singular value of the
+*> 2(n-1)-by-2(n-1) matrix
+*>
+*> Zl = [ kron(a, In-1) -kron(1, S22) ]
+*> [ kron(b, In-1) -kron(1, T22) ] .
+*>
+*> Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
+*> Kronecker product between the matrices X and Y.
+*>
+*> Note that if the default method for computing DIF(i) is wanted
+*> (see DLATDF), then the parameter DIFDRI (see below) should be
+*> changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
+*> See DTGSYL for more details.
+*>
+*> b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
+*>
+*> Suppose U and V are orthogonal transformations such that
+*>
+*> U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
+*> ( 0 S22 ),( 0 T22) n-2
+*> 2 n-2 2 n-2
+*>
+*> and (S11, T11) corresponds to the complex conjugate eigenvalue
+*> pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
+*> that
+*>
+*> U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
+*> ( 0 s22 ) ( 0 t22 )
+*>
+*> where the generalized eigenvalues w = s11/t11 and
+*> conjg(w) = s22/t22.
+*>
+*> Then the reciprocal condition number DIF(i) is bounded by
+*>
+*> min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
+*>
+*> where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
+*> Z1 is the complex 2-by-2 matrix
+*>
+*> Z1 = [ s11 -s22 ]
+*> [ t11 -t22 ],
+*>
+*> This is done by computing (using real arithmetic) the
+*> roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
+*> where Z1**T denotes the transpose of Z1 and det(X) denotes
+*> the determinant of X.
+*>
+*> and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
+*> upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
+*>
+*> Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ]
+*> [ kron(T11**T, In-2) -kron(I2, T22) ]
+*>
+*> Note that if the default method for computing DIF is wanted (see
+*> DLATDF), then the parameter DIFDRI (see below) should be changed
+*> from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
+*> for more details.
+*>
+*> For each eigenvalue/vector specified by SELECT, DIF stores a
+*> Frobenius norm-based estimate of Difl.
+*>
+*> An approximate error bound for the i-th computed eigenvector VL(i) or
+*> VR(i) is given by
+*>
+*> EPS * norm(A, B) / DIF(i).
+*>
+*> See ref. [2-3] for more details and further references.
+*> \endverbatim
+*
+*> \par Contributors:
+* ==================
+*>
+*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
+*> Umea University, S-901 87 Umea, Sweden.
+*
+*> \par References:
+* ================
+*>
+*> \verbatim
+*>
+*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
+*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
+*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
+*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
+*>
+*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
+*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
+*> Estimation: Theory, Algorithms and Software,
+*> Report UMINF - 94.04, Department of Computing Science, Umea
+*> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
+*> Note 87. To appear in Numerical Algorithms, 1996.
+*>
+*> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
+*> for Solving the Generalized Sylvester Equation and Estimating the
+*> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
+*> Department of Computing Science, Umea University, S-901 87 Umea,
+*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
+*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
+*> No 1, 1996.
+*> \endverbatim
+*>
+* =====================================================================
SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
$ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
$ IWORK, INFO )
*
-* -- LAPACK routine (version 3.3.1) --
+* -- LAPACK computational routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* -- April 2011 --
+* November 2011
*
* .. Scalar Arguments ..
CHARACTER HOWMNY, JOB
@@ -18,266 +397,6 @@
$ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
* ..
*
-* Purpose
-* =======
-*
-* DTGSNA estimates reciprocal condition numbers for specified
-* eigenvalues and/or eigenvectors of a matrix pair (A, B) in
-* generalized real Schur canonical form (or of any matrix pair
-* (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
-* Z**T denotes the transpose of Z.
-*
-* (A, B) must be in generalized real Schur form (as returned by DGGES),
-* i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
-* blocks. B is upper triangular.
-*
-*
-* Arguments
-* =========
-*
-* JOB (input) CHARACTER*1
-* Specifies whether condition numbers are required for
-* eigenvalues (S) or eigenvectors (DIF):
-* = 'E': for eigenvalues only (S);
-* = 'V': for eigenvectors only (DIF);
-* = 'B': for both eigenvalues and eigenvectors (S and DIF).
-*
-* HOWMNY (input) CHARACTER*1
-* = 'A': compute condition numbers for all eigenpairs;
-* = 'S': compute condition numbers for selected eigenpairs
-* specified by the array SELECT.
-*
-* SELECT (input) LOGICAL array, dimension (N)
-* If HOWMNY = 'S', SELECT specifies the eigenpairs for which
-* condition numbers are required. To select condition numbers
-* for the eigenpair corresponding to a real eigenvalue w(j),
-* SELECT(j) must be set to .TRUE.. To select condition numbers
-* corresponding to a complex conjugate pair of eigenvalues w(j)
-* and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
-* set to .TRUE..
-* If HOWMNY = 'A', SELECT is not referenced.
-*
-* N (input) INTEGER
-* The order of the square matrix pair (A, B). N >= 0.
-*
-* A (input) DOUBLE PRECISION array, dimension (LDA,N)
-* The upper quasi-triangular matrix A in the pair (A,B).
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,N).
-*
-* B (input) DOUBLE PRECISION array, dimension (LDB,N)
-* The upper triangular matrix B in the pair (A,B).
-*
-* LDB (input) INTEGER
-* The leading dimension of the array B. LDB >= max(1,N).
-*
-* VL (input) DOUBLE PRECISION array, dimension (LDVL,M)
-* If JOB = 'E' or 'B', VL must contain left eigenvectors of
-* (A, B), corresponding to the eigenpairs specified by HOWMNY
-* and SELECT. The eigenvectors must be stored in consecutive
-* columns of VL, as returned by DTGEVC.
-* If JOB = 'V', VL is not referenced.
-*
-* LDVL (input) INTEGER
-* The leading dimension of the array VL. LDVL >= 1.
-* If JOB = 'E' or 'B', LDVL >= N.
-*
-* VR (input) DOUBLE PRECISION array, dimension (LDVR,M)
-* If JOB = 'E' or 'B', VR must contain right eigenvectors of
-* (A, B), corresponding to the eigenpairs specified by HOWMNY
-* and SELECT. The eigenvectors must be stored in consecutive
-* columns ov VR, as returned by DTGEVC.
-* If JOB = 'V', VR is not referenced.
-*
-* LDVR (input) INTEGER
-* The leading dimension of the array VR. LDVR >= 1.
-* If JOB = 'E' or 'B', LDVR >= N.
-*
-* S (output) DOUBLE PRECISION array, dimension (MM)
-* If JOB = 'E' or 'B', the reciprocal condition numbers of the
-* selected eigenvalues, stored in consecutive elements of the
-* array. For a complex conjugate pair of eigenvalues two
-* consecutive elements of S are set to the same value. Thus
-* S(j), DIF(j), and the j-th columns of VL and VR all
-* correspond to the same eigenpair (but not in general the
-* j-th eigenpair, unless all eigenpairs are selected).
-* If JOB = 'V', S is not referenced.
-*
-* DIF (output) DOUBLE PRECISION array, dimension (MM)
-* If JOB = 'V' or 'B', the estimated reciprocal condition
-* numbers of the selected eigenvectors, stored in consecutive
-* elements of the array. For a complex eigenvector two
-* consecutive elements of DIF are set to the same value. If
-* the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
-* is set to 0; this can only occur when the true value would be
-* very small anyway.
-* If JOB = 'E', DIF is not referenced.
-*
-* MM (input) INTEGER
-* The number of elements in the arrays S and DIF. MM >= M.
-*
-* M (output) INTEGER
-* The number of elements of the arrays S and DIF used to store
-* the specified condition numbers; for each selected real
-* eigenvalue one element is used, and for each selected complex
-* conjugate pair of eigenvalues, two elements are used.
-* If HOWMNY = 'A', M is set to N.
-*
-* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
-* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
-*
-* LWORK (input) INTEGER
-* The dimension of the array WORK. LWORK >= max(1,N).
-* If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
-*
-* If LWORK = -1, then a workspace query is assumed; the routine
-* only calculates the optimal size of the WORK array, returns
-* this value as the first entry of the WORK array, and no error
-* message related to LWORK is issued by XERBLA.
-*
-* IWORK (workspace) INTEGER array, dimension (N + 6)
-* If JOB = 'E', IWORK is not referenced.
-*
-* INFO (output) INTEGER
-* =0: Successful exit
-* <0: If INFO = -i, the i-th argument had an illegal value
-*
-*
-* Further Details
-* ===============
-*
-* The reciprocal of the condition number of a generalized eigenvalue
-* w = (a, b) is defined as
-*
-* S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
-*
-* where u and v are the left and right eigenvectors of (A, B)
-* corresponding to w; |z| denotes the absolute value of the complex
-* number, and norm(u) denotes the 2-norm of the vector u.
-* The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
-* of the matrix pair (A, B). If both a and b equal zero, then (A B) is
-* singular and S(I) = -1 is returned.
-*
-* An approximate error bound on the chordal distance between the i-th
-* computed generalized eigenvalue w and the corresponding exact
-* eigenvalue lambda is
-*
-* chord(w, lambda) <= EPS * norm(A, B) / S(I)
-*
-* where EPS is the machine precision.
-*
-* The reciprocal of the condition number DIF(i) of right eigenvector u
-* and left eigenvector v corresponding to the generalized eigenvalue w
-* is defined as follows:
-*
-* a) If the i-th eigenvalue w = (a,b) is real
-*
-* Suppose U and V are orthogonal transformations such that
-*
-* U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
-* ( 0 S22 ),( 0 T22 ) n-1
-* 1 n-1 1 n-1
-*
-* Then the reciprocal condition number DIF(i) is
-*
-* Difl((a, b), (S22, T22)) = sigma-min( Zl ),
-*
-* where sigma-min(Zl) denotes the smallest singular value of the
-* 2(n-1)-by-2(n-1) matrix
-*
-* Zl = [ kron(a, In-1) -kron(1, S22) ]
-* [ kron(b, In-1) -kron(1, T22) ] .
-*
-* Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
-* Kronecker product between the matrices X and Y.
-*
-* Note that if the default method for computing DIF(i) is wanted
-* (see DLATDF), then the parameter DIFDRI (see below) should be
-* changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
-* See DTGSYL for more details.
-*
-* b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
-*
-* Suppose U and V are orthogonal transformations such that
-*
-* U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
-* ( 0 S22 ),( 0 T22) n-2
-* 2 n-2 2 n-2
-*
-* and (S11, T11) corresponds to the complex conjugate eigenvalue
-* pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
-* that
-*
-* U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
-* ( 0 s22 ) ( 0 t22 )
-*
-* where the generalized eigenvalues w = s11/t11 and
-* conjg(w) = s22/t22.
-*
-* Then the reciprocal condition number DIF(i) is bounded by
-*
-* min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
-*
-* where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
-* Z1 is the complex 2-by-2 matrix
-*
-* Z1 = [ s11 -s22 ]
-* [ t11 -t22 ],
-*
-* This is done by computing (using real arithmetic) the
-* roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
-* where Z1**T denotes the transpose of Z1 and det(X) denotes
-* the determinant of X.
-*
-* and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
-* upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
-*
-* Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ]
-* [ kron(T11**T, In-2) -kron(I2, T22) ]
-*
-* Note that if the default method for computing DIF is wanted (see
-* DLATDF), then the parameter DIFDRI (see below) should be changed
-* from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
-* for more details.
-*
-* For each eigenvalue/vector specified by SELECT, DIF stores a
-* Frobenius norm-based estimate of Difl.
-*
-* An approximate error bound for the i-th computed eigenvector VL(i) or
-* VR(i) is given by
-*
-* EPS * norm(A, B) / DIF(i).
-*
-* See ref. [2-3] for more details and further references.
-*
-* Based on contributions by
-* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
-* Umea University, S-901 87 Umea, Sweden.
-*
-* References
-* ==========
-*
-* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
-* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
-* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
-* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
-*
-* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
-* Eigenvalues of a Regular Matrix Pair (A, B) and Condition
-* Estimation: Theory, Algorithms and Software,
-* Report UMINF - 94.04, Department of Computing Science, Umea
-* University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
-* Note 87. To appear in Numerical Algorithms, 1996.
-*
-* [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
-* for Solving the Generalized Sylvester Equation and Estimating the
-* Separation between Regular Matrix Pairs, Report UMINF - 93.23,
-* Department of Computing Science, Umea University, S-901 87 Umea,
-* Sweden, December 1993, Revised April 1994, Also as LAPACK Working
-* Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
-* No 1, 1996.
-*
* =====================================================================
*
* .. Parameters ..