LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ctgsna.f
Go to the documentation of this file.
00001 *> \brief \b CTGSNA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CTGSNA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsna.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsna.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsna.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CTGSNA( 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               DIF( * ), S( * )
00033 *       COMPLEX            A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
00034 *      $                   VR( LDVR, * ), WORK( * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> CTGSNA estimates reciprocal condition numbers for specified
00044 *> eigenvalues and/or eigenvectors of a matrix pair (A, B).
00045 *>
00046 *> (A, B) must be in generalized Schur canonical form, that is, A and
00047 *> B are both upper triangular.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] JOB
00054 *> \verbatim
00055 *>          JOB is CHARACTER*1
00056 *>          Specifies whether condition numbers are required for
00057 *>          eigenvalues (S) or eigenvectors (DIF):
00058 *>          = 'E': for eigenvalues only (S);
00059 *>          = 'V': for eigenvectors only (DIF);
00060 *>          = 'B': for both eigenvalues and eigenvectors (S and DIF).
00061 *> \endverbatim
00062 *>
00063 *> \param[in] HOWMNY
00064 *> \verbatim
00065 *>          HOWMNY is CHARACTER*1
00066 *>          = 'A': compute condition numbers for all eigenpairs;
00067 *>          = 'S': compute condition numbers for selected eigenpairs
00068 *>                 specified by the array SELECT.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] SELECT
00072 *> \verbatim
00073 *>          SELECT is LOGICAL array, dimension (N)
00074 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00075 *>          condition numbers are required. To select condition numbers
00076 *>          for the corresponding j-th eigenvalue and/or eigenvector,
00077 *>          SELECT(j) must be set to .TRUE..
00078 *>          If HOWMNY = 'A', SELECT is not referenced.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] N
00082 *> \verbatim
00083 *>          N is INTEGER
00084 *>          The order of the square matrix pair (A, B). N >= 0.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] A
00088 *> \verbatim
00089 *>          A is COMPLEX array, dimension (LDA,N)
00090 *>          The upper triangular matrix A in the pair (A,B).
00091 *> \endverbatim
00092 *>
00093 *> \param[in] LDA
00094 *> \verbatim
00095 *>          LDA is INTEGER
00096 *>          The leading dimension of the array A. LDA >= max(1,N).
00097 *> \endverbatim
00098 *>
00099 *> \param[in] B
00100 *> \verbatim
00101 *>          B is COMPLEX array, dimension (LDB,N)
00102 *>          The upper triangular matrix B in the pair (A, B).
00103 *> \endverbatim
00104 *>
00105 *> \param[in] LDB
00106 *> \verbatim
00107 *>          LDB is INTEGER
00108 *>          The leading dimension of the array B. LDB >= max(1,N).
00109 *> \endverbatim
00110 *>
00111 *> \param[in] VL
00112 *> \verbatim
00113 *>          VL is COMPLEX array, dimension (LDVL,M)
00114 *>          IF JOB = 'E' or 'B', VL must contain left eigenvectors of
00115 *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
00116 *>          and SELECT.  The eigenvectors must be stored in consecutive
00117 *>          columns of VL, as returned by CTGEVC.
00118 *>          If JOB = 'V', VL is not referenced.
00119 *> \endverbatim
00120 *>
00121 *> \param[in] LDVL
00122 *> \verbatim
00123 *>          LDVL is INTEGER
00124 *>          The leading dimension of the array VL. LDVL >= 1; and
00125 *>          If JOB = 'E' or 'B', LDVL >= N.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] VR
00129 *> \verbatim
00130 *>          VR is COMPLEX array, dimension (LDVR,M)
00131 *>          IF JOB = 'E' or 'B', VR must contain right eigenvectors of
00132 *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
00133 *>          and SELECT.  The eigenvectors must be stored in consecutive
00134 *>          columns of VR, as returned by CTGEVC.
00135 *>          If JOB = 'V', VR is not referenced.
00136 *> \endverbatim
00137 *>
00138 *> \param[in] LDVR
00139 *> \verbatim
00140 *>          LDVR is INTEGER
00141 *>          The leading dimension of the array VR. LDVR >= 1;
00142 *>          If JOB = 'E' or 'B', LDVR >= N.
00143 *> \endverbatim
00144 *>
00145 *> \param[out] S
00146 *> \verbatim
00147 *>          S is REAL array, dimension (MM)
00148 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00149 *>          selected eigenvalues, stored in consecutive elements of the
00150 *>          array.
00151 *>          If JOB = 'V', S is not referenced.
00152 *> \endverbatim
00153 *>
00154 *> \param[out] DIF
00155 *> \verbatim
00156 *>          DIF is REAL array, dimension (MM)
00157 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
00158 *>          numbers of the selected eigenvectors, stored in consecutive
00159 *>          elements of the array.
00160 *>          If the eigenvalues cannot be reordered to compute DIF(j),
00161 *>          DIF(j) is set to 0; this can only occur when the true value
00162 *>          would be very small anyway.
00163 *>          For each eigenvalue/vector specified by SELECT, DIF stores
00164 *>          a Frobenius norm-based estimate of Difl.
00165 *>          If JOB = 'E', DIF is not referenced.
00166 *> \endverbatim
00167 *>
00168 *> \param[in] MM
00169 *> \verbatim
00170 *>          MM is INTEGER
00171 *>          The number of elements in the arrays S and DIF. MM >= M.
00172 *> \endverbatim
00173 *>
00174 *> \param[out] M
00175 *> \verbatim
00176 *>          M is INTEGER
00177 *>          The number of elements of the arrays S and DIF used to store
00178 *>          the specified condition numbers; for each selected eigenvalue
00179 *>          one element is used. If HOWMNY = 'A', M is set to N.
00180 *> \endverbatim
00181 *>
00182 *> \param[out] WORK
00183 *> \verbatim
00184 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
00185 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00186 *> \endverbatim
00187 *>
00188 *> \param[in] LWORK
00189 *> \verbatim
00190 *>          LWORK is INTEGER
00191 *>          The dimension of the array WORK. LWORK >= max(1,N).
00192 *>          If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
00193 *> \endverbatim
00194 *>
00195 *> \param[out] IWORK
00196 *> \verbatim
00197 *>          IWORK is INTEGER array, dimension (N+2)
00198 *>          If JOB = 'E', IWORK is not referenced.
00199 *> \endverbatim
00200 *>
00201 *> \param[out] INFO
00202 *> \verbatim
00203 *>          INFO is INTEGER
00204 *>          = 0: Successful exit
00205 *>          < 0: If INFO = -i, the i-th argument had an illegal value
00206 *> \endverbatim
00207 *
00208 *  Authors:
00209 *  ========
00210 *
00211 *> \author Univ. of Tennessee 
00212 *> \author Univ. of California Berkeley 
00213 *> \author Univ. of Colorado Denver 
00214 *> \author NAG Ltd. 
00215 *
00216 *> \date November 2011
00217 *
00218 *> \ingroup complexOTHERcomputational
00219 *
00220 *> \par Further Details:
00221 *  =====================
00222 *>
00223 *> \verbatim
00224 *>
00225 *>  The reciprocal of the condition number of the i-th generalized
00226 *>  eigenvalue w = (a, b) is defined as
00227 *>
00228 *>          S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
00229 *>
00230 *>  where u and v are the right and left eigenvectors of (A, B)
00231 *>  corresponding to w; |z| denotes the absolute value of the complex
00232 *>  number, and norm(u) denotes the 2-norm of the vector u. The pair
00233 *>  (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
00234 *>  matrix pair (A, B). If both a and b equal zero, then (A,B) is
00235 *>  singular and S(I) = -1 is returned.
00236 *>
00237 *>  An approximate error bound on the chordal distance between the i-th
00238 *>  computed generalized eigenvalue w and the corresponding exact
00239 *>  eigenvalue lambda is
00240 *>
00241 *>          chord(w, lambda) <=   EPS * norm(A, B) / S(I),
00242 *>
00243 *>  where EPS is the machine precision.
00244 *>
00245 *>  The reciprocal of the condition number of the right eigenvector u
00246 *>  and left eigenvector v corresponding to the generalized eigenvalue w
00247 *>  is defined as follows. Suppose
00248 *>
00249 *>                   (A, B) = ( a   *  ) ( b  *  )  1
00250 *>                            ( 0  A22 ),( 0 B22 )  n-1
00251 *>                              1  n-1     1 n-1
00252 *>
00253 *>  Then the reciprocal condition number DIF(I) is
00254 *>
00255 *>          Difl[(a, b), (A22, B22)]  = sigma-min( Zl )
00256 *>
00257 *>  where sigma-min(Zl) denotes the smallest singular value of
00258 *>
00259 *>         Zl = [ kron(a, In-1) -kron(1, A22) ]
00260 *>              [ kron(b, In-1) -kron(1, B22) ].
00261 *>
00262 *>  Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
00263 *>  transpose of X. kron(X, Y) is the Kronecker product between the
00264 *>  matrices X and Y.
00265 *>
00266 *>  We approximate the smallest singular value of Zl with an upper
00267 *>  bound. This is done by CLATDF.
00268 *>
00269 *>  An approximate error bound for a computed eigenvector VL(i) or
00270 *>  VR(i) is given by
00271 *>
00272 *>                      EPS * norm(A, B) / DIF(i).
00273 *>
00274 *>  See ref. [2-3] for more details and further references.
00275 *> \endverbatim
00276 *
00277 *> \par Contributors:
00278 *  ==================
00279 *>
00280 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00281 *>     Umea University, S-901 87 Umea, Sweden.
00282 *
00283 *> \par References:
00284 *  ================
00285 *>
00286 *> \verbatim
00287 *>
00288 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00289 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00290 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00291 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00292 *>
00293 *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
00294 *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
00295 *>      Estimation: Theory, Algorithms and Software, Report
00296 *>      UMINF - 94.04, Department of Computing Science, Umea University,
00297 *>      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
00298 *>      To appear in Numerical Algorithms, 1996.
00299 *>
00300 *>  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
00301 *>      for Solving the Generalized Sylvester Equation and Estimating the
00302 *>      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
00303 *>      Department of Computing Science, Umea University, S-901 87 Umea,
00304 *>      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
00305 *>      Note 75.
00306 *>      To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
00307 *> \endverbatim
00308 *>
00309 *  =====================================================================
00310       SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
00311      $                   LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
00312      $                   IWORK, INFO )
00313 *
00314 *  -- LAPACK computational routine (version 3.4.0) --
00315 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00316 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00317 *     November 2011
00318 *
00319 *     .. Scalar Arguments ..
00320       CHARACTER          HOWMNY, JOB
00321       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
00322 *     ..
00323 *     .. Array Arguments ..
00324       LOGICAL            SELECT( * )
00325       INTEGER            IWORK( * )
00326       REAL               DIF( * ), S( * )
00327       COMPLEX            A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
00328      $                   VR( LDVR, * ), WORK( * )
00329 *     ..
00330 *
00331 *  =====================================================================
00332 *
00333 *     .. Parameters ..
00334       REAL               ZERO, ONE
00335       INTEGER            IDIFJB
00336       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, IDIFJB = 3 )
00337 *     ..
00338 *     .. Local Scalars ..
00339       LOGICAL            LQUERY, SOMCON, WANTBH, WANTDF, WANTS
00340       INTEGER            I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
00341       REAL               BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
00342       COMPLEX            YHAX, YHBX
00343 *     ..
00344 *     .. Local Arrays ..
00345       COMPLEX            DUMMY( 1 ), DUMMY1( 1 )
00346 *     ..
00347 *     .. External Functions ..
00348       LOGICAL            LSAME
00349       REAL               SCNRM2, SLAMCH, SLAPY2
00350       COMPLEX            CDOTC
00351       EXTERNAL           LSAME, SCNRM2, SLAMCH, SLAPY2, CDOTC
00352 *     ..
00353 *     .. External Subroutines ..
00354       EXTERNAL           CGEMV, CLACPY, CTGEXC, CTGSYL, SLABAD, XERBLA
00355 *     ..
00356 *     .. Intrinsic Functions ..
00357       INTRINSIC          ABS, CMPLX, MAX
00358 *     ..
00359 *     .. Executable Statements ..
00360 *
00361 *     Decode and test the input parameters
00362 *
00363       WANTBH = LSAME( JOB, 'B' )
00364       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00365       WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
00366 *
00367       SOMCON = LSAME( HOWMNY, 'S' )
00368 *
00369       INFO = 0
00370       LQUERY = ( LWORK.EQ.-1 )
00371 *
00372       IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
00373          INFO = -1
00374       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00375          INFO = -2
00376       ELSE IF( N.LT.0 ) THEN
00377          INFO = -4
00378       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00379          INFO = -6
00380       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00381          INFO = -8
00382       ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
00383          INFO = -10
00384       ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
00385          INFO = -12
00386       ELSE
00387 *
00388 *        Set M to the number of eigenpairs for which condition numbers
00389 *        are required, and test MM.
00390 *
00391          IF( SOMCON ) THEN
00392             M = 0
00393             DO 10 K = 1, N
00394                IF( SELECT( K ) )
00395      $            M = M + 1
00396    10       CONTINUE
00397          ELSE
00398             M = N
00399          END IF
00400 *
00401          IF( N.EQ.0 ) THEN
00402             LWMIN = 1
00403          ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
00404             LWMIN = 2*N*N
00405          ELSE
00406             LWMIN = N
00407          END IF
00408          WORK( 1 ) = LWMIN
00409 *
00410          IF( MM.LT.M ) THEN
00411             INFO = -15
00412          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00413             INFO = -18
00414          END IF
00415       END IF
00416 *
00417       IF( INFO.NE.0 ) THEN
00418          CALL XERBLA( 'CTGSNA', -INFO )
00419          RETURN
00420       ELSE IF( LQUERY ) THEN
00421          RETURN
00422       END IF
00423 *
00424 *     Quick return if possible
00425 *
00426       IF( N.EQ.0 )
00427      $   RETURN
00428 *
00429 *     Get machine constants
00430 *
00431       EPS = SLAMCH( 'P' )
00432       SMLNUM = SLAMCH( 'S' ) / EPS
00433       BIGNUM = ONE / SMLNUM
00434       CALL SLABAD( SMLNUM, BIGNUM )
00435       KS = 0
00436       DO 20 K = 1, N
00437 *
00438 *        Determine whether condition numbers are required for the k-th
00439 *        eigenpair.
00440 *
00441          IF( SOMCON ) THEN
00442             IF( .NOT.SELECT( K ) )
00443      $         GO TO 20
00444          END IF
00445 *
00446          KS = KS + 1
00447 *
00448          IF( WANTS ) THEN
00449 *
00450 *           Compute the reciprocal condition number of the k-th
00451 *           eigenvalue.
00452 *
00453             RNRM = SCNRM2( N, VR( 1, KS ), 1 )
00454             LNRM = SCNRM2( N, VL( 1, KS ), 1 )
00455             CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), A, LDA,
00456      $                  VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 )
00457             YHAX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 )
00458             CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), B, LDB,
00459      $                  VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 )
00460             YHBX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 )
00461             COND = SLAPY2( ABS( YHAX ), ABS( YHBX ) )
00462             IF( COND.EQ.ZERO ) THEN
00463                S( KS ) = -ONE
00464             ELSE
00465                S( KS ) = COND / ( RNRM*LNRM )
00466             END IF
00467          END IF
00468 *
00469          IF( WANTDF ) THEN
00470             IF( N.EQ.1 ) THEN
00471                DIF( KS ) = SLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
00472             ELSE
00473 *
00474 *              Estimate the reciprocal condition number of the k-th
00475 *              eigenvectors.
00476 *
00477 *              Copy the matrix (A, B) to the array WORK and move the
00478 *              (k,k)th pair to the (1,1) position.
00479 *
00480                CALL CLACPY( 'Full', N, N, A, LDA, WORK, N )
00481                CALL CLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
00482                IFST = K
00483                ILST = 1
00484 *
00485                CALL CTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
00486      $                      N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
00487 *
00488                IF( IERR.GT.0 ) THEN
00489 *
00490 *                 Ill-conditioned problem - swap rejected.
00491 *
00492                   DIF( KS ) = ZERO
00493                ELSE
00494 *
00495 *                 Reordering successful, solve generalized Sylvester
00496 *                 equation for R and L,
00497 *                            A22 * R - L * A11 = A12
00498 *                            B22 * R - L * B11 = B12,
00499 *                 and compute estimate of Difl[(A11,B11), (A22, B22)].
00500 *
00501                   N1 = 1
00502                   N2 = N - N1
00503                   I = N*N + 1
00504                   CALL CTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
00505      $                         N, WORK, N, WORK( N1+1 ), N,
00506      $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
00507      $                         WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
00508      $                         1, IWORK, IERR )
00509                END IF
00510             END IF
00511          END IF
00512 *
00513    20 CONTINUE
00514       WORK( 1 ) = LWMIN
00515       RETURN
00516 *
00517 *     End of CTGSNA
00518 *
00519       END
 All Files Functions