LAPACK 3.3.0

ctgsna.f

Go to the documentation of this file.
00001       SUBROUTINE CTGSNA( 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       REAL               DIF( * ), S( * )
00018       COMPLEX            A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
00019      $                   VR( LDVR, * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  CTGSNA estimates reciprocal condition numbers for specified
00026 *  eigenvalues and/or eigenvectors of a matrix pair (A, B).
00027 *
00028 *  (A, B) must be in generalized Schur canonical form, that is, A and
00029 *  B are both upper triangular.
00030 *
00031 *  Arguments
00032 *  =========
00033 *
00034 *  JOB     (input) CHARACTER*1
00035 *          Specifies whether condition numbers are required for
00036 *          eigenvalues (S) or eigenvectors (DIF):
00037 *          = 'E': for eigenvalues only (S);
00038 *          = 'V': for eigenvectors only (DIF);
00039 *          = 'B': for both eigenvalues and eigenvectors (S and DIF).
00040 *
00041 *  HOWMNY  (input) CHARACTER*1
00042 *          = 'A': compute condition numbers for all eigenpairs;
00043 *          = 'S': compute condition numbers for selected eigenpairs
00044 *                 specified by the array SELECT.
00045 *
00046 *  SELECT  (input) LOGICAL array, dimension (N)
00047 *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00048 *          condition numbers are required. To select condition numbers
00049 *          for the corresponding j-th eigenvalue and/or eigenvector,
00050 *          SELECT(j) must be set to .TRUE..
00051 *          If HOWMNY = 'A', SELECT is not referenced.
00052 *
00053 *  N       (input) INTEGER
00054 *          The order of the square matrix pair (A, B). N >= 0.
00055 *
00056 *  A       (input) COMPLEX array, dimension (LDA,N)
00057 *          The upper triangular matrix A in the pair (A,B).
00058 *
00059 *  LDA     (input) INTEGER
00060 *          The leading dimension of the array A. LDA >= max(1,N).
00061 *
00062 *  B       (input) COMPLEX array, dimension (LDB,N)
00063 *          The upper triangular matrix B in the pair (A, B).
00064 *
00065 *  LDB     (input) INTEGER
00066 *          The leading dimension of the array B. LDB >= max(1,N).
00067 *
00068 *  VL      (input) COMPLEX array, dimension (LDVL,M)
00069 *          IF JOB = 'E' or 'B', VL must contain left eigenvectors of
00070 *          (A, B), corresponding to the eigenpairs specified by HOWMNY
00071 *          and SELECT.  The eigenvectors must be stored in consecutive
00072 *          columns of VL, as returned by CTGEVC.
00073 *          If JOB = 'V', VL is not referenced.
00074 *
00075 *  LDVL    (input) INTEGER
00076 *          The leading dimension of the array VL. LDVL >= 1; and
00077 *          If JOB = 'E' or 'B', LDVL >= N.
00078 *
00079 *  VR      (input) COMPLEX array, dimension (LDVR,M)
00080 *          IF JOB = 'E' or 'B', VR must contain right eigenvectors of
00081 *          (A, B), corresponding to the eigenpairs specified by HOWMNY
00082 *          and SELECT.  The eigenvectors must be stored in consecutive
00083 *          columns of VR, as returned by CTGEVC.
00084 *          If JOB = 'V', VR is not referenced.
00085 *
00086 *  LDVR    (input) INTEGER
00087 *          The leading dimension of the array VR. LDVR >= 1;
00088 *          If JOB = 'E' or 'B', LDVR >= N.
00089 *
00090 *  S       (output) REAL array, dimension (MM)
00091 *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00092 *          selected eigenvalues, stored in consecutive elements of the
00093 *          array.
00094 *          If JOB = 'V', S is not referenced.
00095 *
00096 *  DIF     (output) REAL array, dimension (MM)
00097 *          If JOB = 'V' or 'B', the estimated reciprocal condition
00098 *          numbers of the selected eigenvectors, stored in consecutive
00099 *          elements of the array.
00100 *          If the eigenvalues cannot be reordered to compute DIF(j),
00101 *          DIF(j) is set to 0; this can only occur when the true value
00102 *          would be very small anyway.
00103 *          For each eigenvalue/vector specified by SELECT, DIF stores
00104 *          a Frobenius norm-based estimate of Difl.
00105 *          If JOB = 'E', DIF is not referenced.
00106 *
00107 *  MM      (input) INTEGER
00108 *          The number of elements in the arrays S and DIF. MM >= M.
00109 *
00110 *  M       (output) INTEGER
00111 *          The number of elements of the arrays S and DIF used to store
00112 *          the specified condition numbers; for each selected eigenvalue
00113 *          one element is used. If HOWMNY = 'A', M is set to N.
00114 *
00115 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00116 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00117 *
00118 *  LWORK  (input) INTEGER
00119 *          The dimension of the array WORK. LWORK >= max(1,N).
00120 *          If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
00121 *
00122 *  IWORK   (workspace) INTEGER array, dimension (N+2)
00123 *          If JOB = 'E', IWORK is not referenced.
00124 *
00125 *  INFO    (output) INTEGER
00126 *          = 0: Successful exit
00127 *          < 0: If INFO = -i, the i-th argument had an illegal value
00128 *
00129 *  Further Details
00130 *  ===============
00131 *
00132 *  The reciprocal of the condition number of the i-th generalized
00133 *  eigenvalue w = (a, b) is defined as
00134 *
00135 *          S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v))
00136 *
00137 *  where u and v are the right and left eigenvectors of (A, B)
00138 *  corresponding to w; |z| denotes the absolute value of the complex
00139 *  number, and norm(u) denotes the 2-norm of the vector u. The pair
00140 *  (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the
00141 *  matrix pair (A, B). If both a and b equal zero, then (A,B) is
00142 *  singular and S(I) = -1 is returned.
00143 *
00144 *  An approximate error bound on the chordal distance between the i-th
00145 *  computed generalized eigenvalue w and the corresponding exact
00146 *  eigenvalue lambda is
00147 *
00148 *          chord(w, lambda) <=   EPS * norm(A, B) / S(I),
00149 *
00150 *  where EPS is the machine precision.
00151 *
00152 *  The reciprocal of the condition number of the right eigenvector u
00153 *  and left eigenvector v corresponding to the generalized eigenvalue w
00154 *  is defined as follows. Suppose
00155 *
00156 *                   (A, B) = ( a   *  ) ( b  *  )  1
00157 *                            ( 0  A22 ),( 0 B22 )  n-1
00158 *                              1  n-1     1 n-1
00159 *
00160 *  Then the reciprocal condition number DIF(I) is
00161 *
00162 *          Difl[(a, b), (A22, B22)]  = sigma-min( Zl )
00163 *
00164 *  where sigma-min(Zl) denotes the smallest singular value of
00165 *
00166 *         Zl = [ kron(a, In-1) -kron(1, A22) ]
00167 *              [ kron(b, In-1) -kron(1, B22) ].
00168 *
00169 *  Here In-1 is the identity matrix of size n-1 and X' is the conjugate
00170 *  transpose of X. kron(X, Y) is the Kronecker product between the
00171 *  matrices X and Y.
00172 *
00173 *  We approximate the smallest singular value of Zl with an upper
00174 *  bound. This is done by CLATDF.
00175 *
00176 *  An approximate error bound for a computed eigenvector VL(i) or
00177 *  VR(i) is given by
00178 *
00179 *                      EPS * norm(A, B) / DIF(i).
00180 *
00181 *  See ref. [2-3] for more details and further references.
00182 *
00183 *  Based on contributions by
00184 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00185 *     Umea University, S-901 87 Umea, Sweden.
00186 *
00187 *  References
00188 *  ==========
00189 *
00190 *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00191 *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00192 *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00193 *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00194 *
00195 *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
00196 *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
00197 *      Estimation: Theory, Algorithms and Software, Report
00198 *      UMINF - 94.04, Department of Computing Science, Umea University,
00199 *      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
00200 *      To appear in Numerical Algorithms, 1996.
00201 *
00202 *  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
00203 *      for Solving the Generalized Sylvester Equation and Estimating the
00204 *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
00205 *      Department of Computing Science, Umea University, S-901 87 Umea,
00206 *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
00207 *      Note 75.
00208 *      To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
00209 *
00210 *  =====================================================================
00211 *
00212 *     .. Parameters ..
00213       REAL               ZERO, ONE
00214       INTEGER            IDIFJB
00215       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, IDIFJB = 3 )
00216 *     ..
00217 *     .. Local Scalars ..
00218       LOGICAL            LQUERY, SOMCON, WANTBH, WANTDF, WANTS
00219       INTEGER            I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
00220       REAL               BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
00221       COMPLEX            YHAX, YHBX
00222 *     ..
00223 *     .. Local Arrays ..
00224       COMPLEX            DUMMY( 1 ), DUMMY1( 1 )
00225 *     ..
00226 *     .. External Functions ..
00227       LOGICAL            LSAME
00228       REAL               SCNRM2, SLAMCH, SLAPY2
00229       COMPLEX            CDOTC
00230       EXTERNAL           LSAME, SCNRM2, SLAMCH, SLAPY2, CDOTC
00231 *     ..
00232 *     .. External Subroutines ..
00233       EXTERNAL           CGEMV, CLACPY, CTGEXC, CTGSYL, SLABAD, XERBLA
00234 *     ..
00235 *     .. Intrinsic Functions ..
00236       INTRINSIC          ABS, CMPLX, MAX
00237 *     ..
00238 *     .. Executable Statements ..
00239 *
00240 *     Decode and test the input parameters
00241 *
00242       WANTBH = LSAME( JOB, 'B' )
00243       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00244       WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
00245 *
00246       SOMCON = LSAME( HOWMNY, 'S' )
00247 *
00248       INFO = 0
00249       LQUERY = ( LWORK.EQ.-1 )
00250 *
00251       IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
00252          INFO = -1
00253       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00254          INFO = -2
00255       ELSE IF( N.LT.0 ) THEN
00256          INFO = -4
00257       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00258          INFO = -6
00259       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00260          INFO = -8
00261       ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
00262          INFO = -10
00263       ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
00264          INFO = -12
00265       ELSE
00266 *
00267 *        Set M to the number of eigenpairs for which condition numbers
00268 *        are required, and test MM.
00269 *
00270          IF( SOMCON ) THEN
00271             M = 0
00272             DO 10 K = 1, N
00273                IF( SELECT( K ) )
00274      $            M = M + 1
00275    10       CONTINUE
00276          ELSE
00277             M = N
00278          END IF
00279 *
00280          IF( N.EQ.0 ) THEN
00281             LWMIN = 1
00282          ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
00283             LWMIN = 2*N*N
00284          ELSE
00285             LWMIN = N
00286          END IF
00287          WORK( 1 ) = LWMIN
00288 *
00289          IF( MM.LT.M ) THEN
00290             INFO = -15
00291          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00292             INFO = -18
00293          END IF
00294       END IF
00295 *
00296       IF( INFO.NE.0 ) THEN
00297          CALL XERBLA( 'CTGSNA', -INFO )
00298          RETURN
00299       ELSE IF( LQUERY ) THEN
00300          RETURN
00301       END IF
00302 *
00303 *     Quick return if possible
00304 *
00305       IF( N.EQ.0 )
00306      $   RETURN
00307 *
00308 *     Get machine constants
00309 *
00310       EPS = SLAMCH( 'P' )
00311       SMLNUM = SLAMCH( 'S' ) / EPS
00312       BIGNUM = ONE / SMLNUM
00313       CALL SLABAD( SMLNUM, BIGNUM )
00314       KS = 0
00315       DO 20 K = 1, N
00316 *
00317 *        Determine whether condition numbers are required for the k-th
00318 *        eigenpair.
00319 *
00320          IF( SOMCON ) THEN
00321             IF( .NOT.SELECT( K ) )
00322      $         GO TO 20
00323          END IF
00324 *
00325          KS = KS + 1
00326 *
00327          IF( WANTS ) THEN
00328 *
00329 *           Compute the reciprocal condition number of the k-th
00330 *           eigenvalue.
00331 *
00332             RNRM = SCNRM2( N, VR( 1, KS ), 1 )
00333             LNRM = SCNRM2( N, VL( 1, KS ), 1 )
00334             CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), A, LDA,
00335      $                  VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 )
00336             YHAX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 )
00337             CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), B, LDB,
00338      $                  VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 )
00339             YHBX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 )
00340             COND = SLAPY2( ABS( YHAX ), ABS( YHBX ) )
00341             IF( COND.EQ.ZERO ) THEN
00342                S( KS ) = -ONE
00343             ELSE
00344                S( KS ) = COND / ( RNRM*LNRM )
00345             END IF
00346          END IF
00347 *
00348          IF( WANTDF ) THEN
00349             IF( N.EQ.1 ) THEN
00350                DIF( KS ) = SLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
00351             ELSE
00352 *
00353 *              Estimate the reciprocal condition number of the k-th
00354 *              eigenvectors.
00355 *
00356 *              Copy the matrix (A, B) to the array WORK and move the
00357 *              (k,k)th pair to the (1,1) position.
00358 *
00359                CALL CLACPY( 'Full', N, N, A, LDA, WORK, N )
00360                CALL CLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
00361                IFST = K
00362                ILST = 1
00363 *
00364                CALL CTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
00365      $                      N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
00366 *
00367                IF( IERR.GT.0 ) THEN
00368 *
00369 *                 Ill-conditioned problem - swap rejected.
00370 *
00371                   DIF( KS ) = ZERO
00372                ELSE
00373 *
00374 *                 Reordering successful, solve generalized Sylvester
00375 *                 equation for R and L,
00376 *                            A22 * R - L * A11 = A12
00377 *                            B22 * R - L * B11 = B12,
00378 *                 and compute estimate of Difl[(A11,B11), (A22, B22)].
00379 *
00380                   N1 = 1
00381                   N2 = N - N1
00382                   I = N*N + 1
00383                   CALL CTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
00384      $                         N, WORK, N, WORK( N1+1 ), N,
00385      $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
00386      $                         WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
00387      $                         1, IWORK, IERR )
00388                END IF
00389             END IF
00390          END IF
00391 *
00392    20 CONTINUE
00393       WORK( 1 ) = LWMIN
00394       RETURN
00395 *
00396 *     End of CTGSNA
00397 *
00398       END
 All Files Functions