LAPACK 3.3.0

strsna.f

Go to the documentation of this file.
00001       SUBROUTINE STRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00002      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
00003      $                   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 *     Modified to call SLACN2 in place of SLACON, 7 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       INTEGER            IWORK( * )
00019       REAL               S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
00020      $                   VR( LDVR, * ), WORK( LDWORK, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  STRSNA estimates reciprocal condition numbers for specified
00027 *  eigenvalues and/or right eigenvectors of a real upper
00028 *  quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
00029 *  orthogonal).
00030 *
00031 *  T must be in Schur canonical form (as returned by SHSEQR), that is,
00032 *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
00033 *  2-by-2 diagonal block has its diagonal elements equal and its
00034 *  off-diagonal elements of opposite sign.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  JOB     (input) CHARACTER*1
00040 *          Specifies whether condition numbers are required for
00041 *          eigenvalues (S) or eigenvectors (SEP):
00042 *          = 'E': for eigenvalues only (S);
00043 *          = 'V': for eigenvectors only (SEP);
00044 *          = 'B': for both eigenvalues and eigenvectors (S and SEP).
00045 *
00046 *  HOWMNY  (input) CHARACTER*1
00047 *          = 'A': compute condition numbers for all eigenpairs;
00048 *          = 'S': compute condition numbers for selected eigenpairs
00049 *                 specified by the array SELECT.
00050 *
00051 *  SELECT  (input) LOGICAL array, dimension (N)
00052 *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00053 *          condition numbers are required. To select condition numbers
00054 *          for the eigenpair corresponding to a real eigenvalue w(j),
00055 *          SELECT(j) must be set to .TRUE.. To select condition numbers
00056 *          corresponding to a complex conjugate pair of eigenvalues w(j)
00057 *          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
00058 *          set to .TRUE..
00059 *          If HOWMNY = 'A', SELECT is not referenced.
00060 *
00061 *  N       (input) INTEGER
00062 *          The order of the matrix T. N >= 0.
00063 *
00064 *  T       (input) REAL array, dimension (LDT,N)
00065 *          The upper quasi-triangular matrix T, in Schur canonical form.
00066 *
00067 *  LDT     (input) INTEGER
00068 *          The leading dimension of the array T. LDT >= max(1,N).
00069 *
00070 *  VL      (input) REAL array, dimension (LDVL,M)
00071 *          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
00072 *          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
00073 *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00074 *          must be stored in consecutive columns of VL, as returned by
00075 *          SHSEIN or STREVC.
00076 *          If JOB = 'V', VL is not referenced.
00077 *
00078 *  LDVL    (input) INTEGER
00079 *          The leading dimension of the array VL.
00080 *          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
00081 *
00082 *  VR      (input) REAL array, dimension (LDVR,M)
00083 *          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
00084 *          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
00085 *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00086 *          must be stored in consecutive columns of VR, as returned by
00087 *          SHSEIN or STREVC.
00088 *          If JOB = 'V', VR is not referenced.
00089 *
00090 *  LDVR    (input) INTEGER
00091 *          The leading dimension of the array VR.
00092 *          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
00093 *
00094 *  S       (output) REAL array, dimension (MM)
00095 *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00096 *          selected eigenvalues, stored in consecutive elements of the
00097 *          array. For a complex conjugate pair of eigenvalues two
00098 *          consecutive elements of S are set to the same value. Thus
00099 *          S(j), SEP(j), and the j-th columns of VL and VR all
00100 *          correspond to the same eigenpair (but not in general the
00101 *          j-th eigenpair, unless all eigenpairs are selected).
00102 *          If JOB = 'V', S is not referenced.
00103 *
00104 *  SEP     (output) REAL array, dimension (MM)
00105 *          If JOB = 'V' or 'B', the estimated reciprocal condition
00106 *          numbers of the selected eigenvectors, stored in consecutive
00107 *          elements of the array. For a complex eigenvector two
00108 *          consecutive elements of SEP are set to the same value. If
00109 *          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
00110 *          is set to 0; this can only occur when the true value would be
00111 *          very small anyway.
00112 *          If JOB = 'E', SEP is not referenced.
00113 *
00114 *  MM      (input) INTEGER
00115 *          The number of elements in the arrays S (if JOB = 'E' or 'B')
00116 *           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
00117 *
00118 *  M       (output) INTEGER
00119 *          The number of elements of the arrays S and/or SEP actually
00120 *          used to store the estimated condition numbers.
00121 *          If HOWMNY = 'A', M is set to N.
00122 *
00123 *  WORK    (workspace) REAL array, dimension (LDWORK,N+6)
00124 *          If JOB = 'E', WORK is not referenced.
00125 *
00126 *  LDWORK  (input) INTEGER
00127 *          The leading dimension of the array WORK.
00128 *          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
00129 *
00130 *  IWORK   (workspace) INTEGER array, dimension (2*(N-1))
00131 *          If JOB = 'E', IWORK is not referenced.
00132 *
00133 *  INFO    (output) INTEGER
00134 *          = 0: successful exit
00135 *          < 0: if INFO = -i, the i-th argument had an illegal value
00136 *
00137 *  Further Details
00138 *  ===============
00139 *
00140 *  The reciprocal of the condition number of an eigenvalue lambda is
00141 *  defined as
00142 *
00143 *          S(lambda) = |v'*u| / (norm(u)*norm(v))
00144 *
00145 *  where u and v are the right and left eigenvectors of T corresponding
00146 *  to lambda; v' denotes the conjugate-transpose of v, and norm(u)
00147 *  denotes the Euclidean norm. These reciprocal condition numbers always
00148 *  lie between zero (very badly conditioned) and one (very well
00149 *  conditioned). If n = 1, S(lambda) is defined to be 1.
00150 *
00151 *  An approximate error bound for a computed eigenvalue W(i) is given by
00152 *
00153 *                      EPS * norm(T) / S(i)
00154 *
00155 *  where EPS is the machine precision.
00156 *
00157 *  The reciprocal of the condition number of the right eigenvector u
00158 *  corresponding to lambda is defined as follows. Suppose
00159 *
00160 *              T = ( lambda  c  )
00161 *                  (   0    T22 )
00162 *
00163 *  Then the reciprocal condition number is
00164 *
00165 *          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
00166 *
00167 *  where sigma-min denotes the smallest singular value. We approximate
00168 *  the smallest singular value by the reciprocal of an estimate of the
00169 *  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
00170 *  defined to be abs(T(1,1)).
00171 *
00172 *  An approximate error bound for a computed right eigenvector VR(i)
00173 *  is given by
00174 *
00175 *                      EPS * norm(T) / SEP(i)
00176 *
00177 *  =====================================================================
00178 *
00179 *     .. Parameters ..
00180       REAL               ZERO, ONE, TWO
00181       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 )
00182 *     ..
00183 *     .. Local Scalars ..
00184       LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
00185       INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
00186       REAL               BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
00187      $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
00188 *     ..
00189 *     .. Local Arrays ..
00190       INTEGER            ISAVE( 3 )
00191       REAL               DUMMY( 1 )
00192 *     ..
00193 *     .. External Functions ..
00194       LOGICAL            LSAME
00195       REAL               SDOT, SLAMCH, SLAPY2, SNRM2
00196       EXTERNAL           LSAME, SDOT, SLAMCH, SLAPY2, SNRM2
00197 *     ..
00198 *     .. External Subroutines ..
00199       EXTERNAL           SLABAD, SLACN2, SLACPY, SLAQTR, STREXC, XERBLA
00200 *     ..
00201 *     .. Intrinsic Functions ..
00202       INTRINSIC          ABS, MAX, SQRT
00203 *     ..
00204 *     .. Executable Statements ..
00205 *
00206 *     Decode and test the input parameters
00207 *
00208       WANTBH = LSAME( JOB, 'B' )
00209       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00210       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00211 *
00212       SOMCON = LSAME( HOWMNY, 'S' )
00213 *
00214       INFO = 0
00215       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
00216          INFO = -1
00217       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00218          INFO = -2
00219       ELSE IF( N.LT.0 ) THEN
00220          INFO = -4
00221       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00222          INFO = -6
00223       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
00224          INFO = -8
00225       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
00226          INFO = -10
00227       ELSE
00228 *
00229 *        Set M to the number of eigenpairs for which condition numbers
00230 *        are required, and test MM.
00231 *
00232          IF( SOMCON ) THEN
00233             M = 0
00234             PAIR = .FALSE.
00235             DO 10 K = 1, N
00236                IF( PAIR ) THEN
00237                   PAIR = .FALSE.
00238                ELSE
00239                   IF( K.LT.N ) THEN
00240                      IF( T( K+1, K ).EQ.ZERO ) THEN
00241                         IF( SELECT( K ) )
00242      $                     M = M + 1
00243                      ELSE
00244                         PAIR = .TRUE.
00245                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
00246      $                     M = M + 2
00247                      END IF
00248                   ELSE
00249                      IF( SELECT( N ) )
00250      $                  M = M + 1
00251                   END IF
00252                END IF
00253    10       CONTINUE
00254          ELSE
00255             M = N
00256          END IF
00257 *
00258          IF( MM.LT.M ) THEN
00259             INFO = -13
00260          ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
00261             INFO = -16
00262          END IF
00263       END IF
00264       IF( INFO.NE.0 ) THEN
00265          CALL XERBLA( 'STRSNA', -INFO )
00266          RETURN
00267       END IF
00268 *
00269 *     Quick return if possible
00270 *
00271       IF( N.EQ.0 )
00272      $   RETURN
00273 *
00274       IF( N.EQ.1 ) THEN
00275          IF( SOMCON ) THEN
00276             IF( .NOT.SELECT( 1 ) )
00277      $         RETURN
00278          END IF
00279          IF( WANTS )
00280      $      S( 1 ) = ONE
00281          IF( WANTSP )
00282      $      SEP( 1 ) = ABS( T( 1, 1 ) )
00283          RETURN
00284       END IF
00285 *
00286 *     Get machine constants
00287 *
00288       EPS = SLAMCH( 'P' )
00289       SMLNUM = SLAMCH( 'S' ) / EPS
00290       BIGNUM = ONE / SMLNUM
00291       CALL SLABAD( SMLNUM, BIGNUM )
00292 *
00293       KS = 0
00294       PAIR = .FALSE.
00295       DO 60 K = 1, N
00296 *
00297 *        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
00298 *
00299          IF( PAIR ) THEN
00300             PAIR = .FALSE.
00301             GO TO 60
00302          ELSE
00303             IF( K.LT.N )
00304      $         PAIR = T( K+1, K ).NE.ZERO
00305          END IF
00306 *
00307 *        Determine whether condition numbers are required for the k-th
00308 *        eigenpair.
00309 *
00310          IF( SOMCON ) THEN
00311             IF( PAIR ) THEN
00312                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
00313      $            GO TO 60
00314             ELSE
00315                IF( .NOT.SELECT( K ) )
00316      $            GO TO 60
00317             END IF
00318          END IF
00319 *
00320          KS = KS + 1
00321 *
00322          IF( WANTS ) THEN
00323 *
00324 *           Compute the reciprocal condition number of the k-th
00325 *           eigenvalue.
00326 *
00327             IF( .NOT.PAIR ) THEN
00328 *
00329 *              Real eigenvalue.
00330 *
00331                PROD = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00332                RNRM = SNRM2( N, VR( 1, KS ), 1 )
00333                LNRM = SNRM2( N, VL( 1, KS ), 1 )
00334                S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
00335             ELSE
00336 *
00337 *              Complex eigenvalue.
00338 *
00339                PROD1 = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00340                PROD1 = PROD1 + SDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
00341      $                 1 )
00342                PROD2 = SDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
00343                PROD2 = PROD2 - SDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
00344      $                 1 )
00345                RNRM = SLAPY2( SNRM2( N, VR( 1, KS ), 1 ),
00346      $                SNRM2( N, VR( 1, KS+1 ), 1 ) )
00347                LNRM = SLAPY2( SNRM2( N, VL( 1, KS ), 1 ),
00348      $                SNRM2( N, VL( 1, KS+1 ), 1 ) )
00349                COND = SLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
00350                S( KS ) = COND
00351                S( KS+1 ) = COND
00352             END IF
00353          END IF
00354 *
00355          IF( WANTSP ) THEN
00356 *
00357 *           Estimate the reciprocal condition number of the k-th
00358 *           eigenvector.
00359 *
00360 *           Copy the matrix T to the array WORK and swap the diagonal
00361 *           block beginning at T(k,k) to the (1,1) position.
00362 *
00363             CALL SLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
00364             IFST = K
00365             ILST = 1
00366             CALL STREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
00367      $                   WORK( 1, N+1 ), IERR )
00368 *
00369             IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
00370 *
00371 *              Could not swap because blocks not well separated
00372 *
00373                SCALE = ONE
00374                EST = BIGNUM
00375             ELSE
00376 *
00377 *              Reordering successful
00378 *
00379                IF( WORK( 2, 1 ).EQ.ZERO ) THEN
00380 *
00381 *                 Form C = T22 - lambda*I in WORK(2:N,2:N).
00382 *
00383                   DO 20 I = 2, N
00384                      WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
00385    20             CONTINUE
00386                   N2 = 1
00387                   NN = N - 1
00388                ELSE
00389 *
00390 *                 Triangularize the 2 by 2 block by unitary
00391 *                 transformation U = [  cs   i*ss ]
00392 *                                    [ i*ss   cs  ].
00393 *                 such that the (1,1) position of WORK is complex
00394 *                 eigenvalue lambda with positive imaginary part. (2,2)
00395 *                 position of WORK is the complex eigenvalue lambda
00396 *                 with negative imaginary  part.
00397 *
00398                   MU = SQRT( ABS( WORK( 1, 2 ) ) )*
00399      $                 SQRT( ABS( WORK( 2, 1 ) ) )
00400                   DELTA = SLAPY2( MU, WORK( 2, 1 ) )
00401                   CS = MU / DELTA
00402                   SN = -WORK( 2, 1 ) / DELTA
00403 *
00404 *                 Form
00405 *
00406 *                 C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
00407 *                                        [   mu                     ]
00408 *                                        [         ..               ]
00409 *                                        [             ..           ]
00410 *                                        [                  mu      ]
00411 *                 where C' is conjugate transpose of complex matrix C,
00412 *                 and RWORK is stored starting in the N+1-st column of
00413 *                 WORK.
00414 *
00415                   DO 30 J = 3, N
00416                      WORK( 2, J ) = CS*WORK( 2, J )
00417                      WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
00418    30             CONTINUE
00419                   WORK( 2, 2 ) = ZERO
00420 *
00421                   WORK( 1, N+1 ) = TWO*MU
00422                   DO 40 I = 2, N - 1
00423                      WORK( I, N+1 ) = SN*WORK( 1, I+1 )
00424    40             CONTINUE
00425                   N2 = 2
00426                   NN = 2*( N-1 )
00427                END IF
00428 *
00429 *              Estimate norm(inv(C'))
00430 *
00431                EST = ZERO
00432                KASE = 0
00433    50          CONTINUE
00434                CALL SLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
00435      $                      EST, KASE, ISAVE )
00436                IF( KASE.NE.0 ) THEN
00437                   IF( KASE.EQ.1 ) THEN
00438                      IF( N2.EQ.1 ) THEN
00439 *
00440 *                       Real eigenvalue: solve C'*x = scale*c.
00441 *
00442                         CALL SLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
00443      $                               LDWORK, DUMMY, DUMM, SCALE,
00444      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00445      $                               IERR )
00446                      ELSE
00447 *
00448 *                       Complex eigenvalue: solve
00449 *                       C'*(p+iq) = scale*(c+id) in real arithmetic.
00450 *
00451                         CALL SLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
00452      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
00453      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00454      $                               IERR )
00455                      END IF
00456                   ELSE
00457                      IF( N2.EQ.1 ) THEN
00458 *
00459 *                       Real eigenvalue: solve C*x = scale*c.
00460 *
00461                         CALL SLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
00462      $                               LDWORK, DUMMY, DUMM, SCALE,
00463      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00464      $                               IERR )
00465                      ELSE
00466 *
00467 *                       Complex eigenvalue: solve
00468 *                       C*(p+iq) = scale*(c+id) in real arithmetic.
00469 *
00470                         CALL SLAQTR( .FALSE., .FALSE., N-1,
00471      $                               WORK( 2, 2 ), LDWORK,
00472      $                               WORK( 1, N+1 ), MU, SCALE,
00473      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
00474      $                               IERR )
00475 *
00476                      END IF
00477                   END IF
00478 *
00479                   GO TO 50
00480                END IF
00481             END IF
00482 *
00483             SEP( KS ) = SCALE / MAX( EST, SMLNUM )
00484             IF( PAIR )
00485      $         SEP( KS+1 ) = SEP( KS )
00486          END IF
00487 *
00488          IF( PAIR )
00489      $      KS = KS + 1
00490 *
00491    60 CONTINUE
00492       RETURN
00493 *
00494 *     End of STRSNA
00495 *
00496       END
 All Files Functions