LAPACK 3.3.1
Linear Algebra PACKage

ctrsna.f

Go to the documentation of this file.
00001       SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00002      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
00003      $                   INFO )
00004 *
00005 *  -- LAPACK routine (version 3.3.1) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *  -- April 2011                                                      --
00009 *
00010 *     Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH.
00011 *
00012 *     .. Scalar Arguments ..
00013       CHARACTER          HOWMNY, JOB
00014       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
00015 *     ..
00016 *     .. Array Arguments ..
00017       LOGICAL            SELECT( * )
00018       REAL               RWORK( * ), S( * ), SEP( * )
00019       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00020      $                   WORK( LDWORK, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CTRSNA estimates reciprocal condition numbers for specified
00027 *  eigenvalues and/or right eigenvectors of a complex upper triangular
00028 *  matrix T (or of any matrix Q*T*Q**H with Q unitary).
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  JOB     (input) CHARACTER*1
00034 *          Specifies whether condition numbers are required for
00035 *          eigenvalues (S) or eigenvectors (SEP):
00036 *          = 'E': for eigenvalues only (S);
00037 *          = 'V': for eigenvectors only (SEP);
00038 *          = 'B': for both eigenvalues and eigenvectors (S and SEP).
00039 *
00040 *  HOWMNY  (input) CHARACTER*1
00041 *          = 'A': compute condition numbers for all eigenpairs;
00042 *          = 'S': compute condition numbers for selected eigenpairs
00043 *                 specified by the array SELECT.
00044 *
00045 *  SELECT  (input) LOGICAL array, dimension (N)
00046 *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00047 *          condition numbers are required. To select condition numbers
00048 *          for the j-th eigenpair, SELECT(j) must be set to .TRUE..
00049 *          If HOWMNY = 'A', SELECT is not referenced.
00050 *
00051 *  N       (input) INTEGER
00052 *          The order of the matrix T. N >= 0.
00053 *
00054 *  T       (input) COMPLEX array, dimension (LDT,N)
00055 *          The upper triangular matrix T.
00056 *
00057 *  LDT     (input) INTEGER
00058 *          The leading dimension of the array T. LDT >= max(1,N).
00059 *
00060 *  VL      (input) COMPLEX array, dimension (LDVL,M)
00061 *          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
00062 *          (or of any Q*T*Q**H with Q unitary), corresponding to the
00063 *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00064 *          must be stored in consecutive columns of VL, as returned by
00065 *          CHSEIN or CTREVC.
00066 *          If JOB = 'V', VL is not referenced.
00067 *
00068 *  LDVL    (input) INTEGER
00069 *          The leading dimension of the array VL.
00070 *          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
00071 *
00072 *  VR      (input) COMPLEX array, dimension (LDVR,M)
00073 *          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
00074 *          (or of any Q*T*Q**H with Q unitary), corresponding to the
00075 *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00076 *          must be stored in consecutive columns of VR, as returned by
00077 *          CHSEIN or CTREVC.
00078 *          If JOB = 'V', VR is not referenced.
00079 *
00080 *  LDVR    (input) INTEGER
00081 *          The leading dimension of the array VR.
00082 *          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
00083 *
00084 *  S       (output) REAL array, dimension (MM)
00085 *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00086 *          selected eigenvalues, stored in consecutive elements of the
00087 *          array. Thus S(j), SEP(j), and the j-th columns of VL and VR
00088 *          all correspond to the same eigenpair (but not in general the
00089 *          j-th eigenpair, unless all eigenpairs are selected).
00090 *          If JOB = 'V', S is not referenced.
00091 *
00092 *  SEP     (output) REAL array, dimension (MM)
00093 *          If JOB = 'V' or 'B', the estimated reciprocal condition
00094 *          numbers of the selected eigenvectors, stored in consecutive
00095 *          elements of the array.
00096 *          If JOB = 'E', SEP is not referenced.
00097 *
00098 *  MM      (input) INTEGER
00099 *          The number of elements in the arrays S (if JOB = 'E' or 'B')
00100 *           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
00101 *
00102 *  M       (output) INTEGER
00103 *          The number of elements of the arrays S and/or SEP actually
00104 *          used to store the estimated condition numbers.
00105 *          If HOWMNY = 'A', M is set to N.
00106 *
00107 *  WORK    (workspace) COMPLEX array, dimension (LDWORK,N+6)
00108 *          If JOB = 'E', WORK is not referenced.
00109 *
00110 *  LDWORK  (input) INTEGER
00111 *          The leading dimension of the array WORK.
00112 *          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
00113 *
00114 *  RWORK   (workspace) REAL array, dimension (N)
00115 *          If JOB = 'E', RWORK is not referenced.
00116 *
00117 *  INFO    (output) INTEGER
00118 *          = 0: successful exit
00119 *          < 0: if INFO = -i, the i-th argument had an illegal value
00120 *
00121 *  Further Details
00122 *  ===============
00123 *
00124 *  The reciprocal of the condition number of an eigenvalue lambda is
00125 *  defined as
00126 *
00127 *          S(lambda) = |v**H*u| / (norm(u)*norm(v))
00128 *
00129 *  where u and v are the right and left eigenvectors of T corresponding
00130 *  to lambda; v**H denotes the conjugate transpose of v, and norm(u)
00131 *  denotes the Euclidean norm. These reciprocal condition numbers always
00132 *  lie between zero (very badly conditioned) and one (very well
00133 *  conditioned). If n = 1, S(lambda) is defined to be 1.
00134 *
00135 *  An approximate error bound for a computed eigenvalue W(i) is given by
00136 *
00137 *                      EPS * norm(T) / S(i)
00138 *
00139 *  where EPS is the machine precision.
00140 *
00141 *  The reciprocal of the condition number of the right eigenvector u
00142 *  corresponding to lambda is defined as follows. Suppose
00143 *
00144 *              T = ( lambda  c  )
00145 *                  (   0    T22 )
00146 *
00147 *  Then the reciprocal condition number is
00148 *
00149 *          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
00150 *
00151 *  where sigma-min denotes the smallest singular value. We approximate
00152 *  the smallest singular value by the reciprocal of an estimate of the
00153 *  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
00154 *  defined to be abs(T(1,1)).
00155 *
00156 *  An approximate error bound for a computed right eigenvector VR(i)
00157 *  is given by
00158 *
00159 *                      EPS * norm(T) / SEP(i)
00160 *
00161 *  =====================================================================
00162 *
00163 *     .. Parameters ..
00164       REAL               ZERO, ONE
00165       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0+0 )
00166 *     ..
00167 *     .. Local Scalars ..
00168       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
00169       CHARACTER          NORMIN
00170       INTEGER            I, IERR, IX, J, K, KASE, KS
00171       REAL               BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
00172      $                   XNORM
00173       COMPLEX            CDUM, PROD
00174 *     ..
00175 *     .. Local Arrays ..
00176       INTEGER            ISAVE( 3 )
00177       COMPLEX            DUMMY( 1 )
00178 *     ..
00179 *     .. External Functions ..
00180       LOGICAL            LSAME
00181       INTEGER            ICAMAX
00182       REAL               SCNRM2, SLAMCH
00183       COMPLEX            CDOTC
00184       EXTERNAL           LSAME, ICAMAX, SCNRM2, SLAMCH, CDOTC
00185 *     ..
00186 *     .. External Subroutines ..
00187       EXTERNAL           CLACN2, CLACPY, CLATRS, CSRSCL, CTREXC, SLABAD,
00188      $                   XERBLA
00189 *     ..
00190 *     .. Intrinsic Functions ..
00191       INTRINSIC          ABS, AIMAG, MAX, REAL
00192 *     ..
00193 *     .. Statement Functions ..
00194       REAL               CABS1
00195 *     ..
00196 *     .. Statement Function definitions ..
00197       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00198 *     ..
00199 *     .. Executable Statements ..
00200 *
00201 *     Decode and test the input parameters
00202 *
00203       WANTBH = LSAME( JOB, 'B' )
00204       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00205       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00206 *
00207       SOMCON = LSAME( HOWMNY, 'S' )
00208 *
00209 *     Set M to the number of eigenpairs for which condition numbers are
00210 *     to be computed.
00211 *
00212       IF( SOMCON ) THEN
00213          M = 0
00214          DO 10 J = 1, N
00215             IF( SELECT( J ) )
00216      $         M = M + 1
00217    10    CONTINUE
00218       ELSE
00219          M = N
00220       END IF
00221 *
00222       INFO = 0
00223       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
00224          INFO = -1
00225       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00226          INFO = -2
00227       ELSE IF( N.LT.0 ) THEN
00228          INFO = -4
00229       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00230          INFO = -6
00231       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
00232          INFO = -8
00233       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
00234          INFO = -10
00235       ELSE IF( MM.LT.M ) THEN
00236          INFO = -13
00237       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
00238          INFO = -16
00239       END IF
00240       IF( INFO.NE.0 ) THEN
00241          CALL XERBLA( 'CTRSNA', -INFO )
00242          RETURN
00243       END IF
00244 *
00245 *     Quick return if possible
00246 *
00247       IF( N.EQ.0 )
00248      $   RETURN
00249 *
00250       IF( N.EQ.1 ) THEN
00251          IF( SOMCON ) THEN
00252             IF( .NOT.SELECT( 1 ) )
00253      $         RETURN
00254          END IF
00255          IF( WANTS )
00256      $      S( 1 ) = ONE
00257          IF( WANTSP )
00258      $      SEP( 1 ) = ABS( T( 1, 1 ) )
00259          RETURN
00260       END IF
00261 *
00262 *     Get machine constants
00263 *
00264       EPS = SLAMCH( 'P' )
00265       SMLNUM = SLAMCH( 'S' ) / EPS
00266       BIGNUM = ONE / SMLNUM
00267       CALL SLABAD( SMLNUM, BIGNUM )
00268 *
00269       KS = 1
00270       DO 50 K = 1, N
00271 *
00272          IF( SOMCON ) THEN
00273             IF( .NOT.SELECT( K ) )
00274      $         GO TO 50
00275          END IF
00276 *
00277          IF( WANTS ) THEN
00278 *
00279 *           Compute the reciprocal condition number of the k-th
00280 *           eigenvalue.
00281 *
00282             PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00283             RNRM = SCNRM2( N, VR( 1, KS ), 1 )
00284             LNRM = SCNRM2( N, VL( 1, KS ), 1 )
00285             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
00286 *
00287          END IF
00288 *
00289          IF( WANTSP ) THEN
00290 *
00291 *           Estimate the reciprocal condition number of the k-th
00292 *           eigenvector.
00293 *
00294 *           Copy the matrix T to the array WORK and swap the k-th
00295 *           diagonal element to the (1,1) position.
00296 *
00297             CALL CLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
00298             CALL CTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
00299 *
00300 *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
00301 *
00302             DO 20 I = 2, N
00303                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
00304    20       CONTINUE
00305 *
00306 *           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
00307 *           and (N+1)th columns of WORK are used to store work vectors.
00308 *
00309             SEP( KS ) = ZERO
00310             EST = ZERO
00311             KASE = 0
00312             NORMIN = 'N'
00313    30       CONTINUE
00314             CALL CLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
00315 *
00316             IF( KASE.NE.0 ) THEN
00317                IF( KASE.EQ.1 ) THEN
00318 *
00319 *                 Solve C**H*x = scale*b
00320 *
00321                   CALL CLATRS( 'Upper', 'Conjugate transpose',
00322      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
00323      $                         LDWORK, WORK, SCALE, RWORK, IERR )
00324                ELSE
00325 *
00326 *                 Solve C*x = scale*b
00327 *
00328                   CALL CLATRS( 'Upper', 'No transpose', 'Nonunit',
00329      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
00330      $                         SCALE, RWORK, IERR )
00331                END IF
00332                NORMIN = 'Y'
00333                IF( SCALE.NE.ONE ) THEN
00334 *
00335 *                 Multiply by 1/SCALE if doing so will not cause
00336 *                 overflow.
00337 *
00338                   IX = ICAMAX( N-1, WORK, 1 )
00339                   XNORM = CABS1( WORK( IX, 1 ) )
00340                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
00341      $               GO TO 40
00342                   CALL CSRSCL( N, SCALE, WORK, 1 )
00343                END IF
00344                GO TO 30
00345             END IF
00346 *
00347             SEP( KS ) = ONE / MAX( EST, SMLNUM )
00348          END IF
00349 *
00350    40    CONTINUE
00351          KS = KS + 1
00352    50 CONTINUE
00353       RETURN
00354 *
00355 *     End of CTRSNA
00356 *
00357       END
 All Files Functions