001:       SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
002:      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
003:      $                   INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH.
011: *
012: *     .. Scalar Arguments ..
013:       CHARACTER          HOWMNY, JOB
014:       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
015: *     ..
016: *     .. Array Arguments ..
017:       LOGICAL            SELECT( * )
018:       REAL               RWORK( * ), S( * ), SEP( * )
019:       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
020:      $                   WORK( LDWORK, * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  CTRSNA estimates reciprocal condition numbers for specified
027: *  eigenvalues and/or right eigenvectors of a complex upper triangular
028: *  matrix T (or of any matrix Q*T*Q**H with Q unitary).
029: *
030: *  Arguments
031: *  =========
032: *
033: *  JOB     (input) CHARACTER*1
034: *          Specifies whether condition numbers are required for
035: *          eigenvalues (S) or eigenvectors (SEP):
036: *          = 'E': for eigenvalues only (S);
037: *          = 'V': for eigenvectors only (SEP);
038: *          = 'B': for both eigenvalues and eigenvectors (S and SEP).
039: *
040: *  HOWMNY  (input) CHARACTER*1
041: *          = 'A': compute condition numbers for all eigenpairs;
042: *          = 'S': compute condition numbers for selected eigenpairs
043: *                 specified by the array SELECT.
044: *
045: *  SELECT  (input) LOGICAL array, dimension (N)
046: *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
047: *          condition numbers are required. To select condition numbers
048: *          for the j-th eigenpair, SELECT(j) must be set to .TRUE..
049: *          If HOWMNY = 'A', SELECT is not referenced.
050: *
051: *  N       (input) INTEGER
052: *          The order of the matrix T. N >= 0.
053: *
054: *  T       (input) COMPLEX array, dimension (LDT,N)
055: *          The upper triangular matrix T.
056: *
057: *  LDT     (input) INTEGER
058: *          The leading dimension of the array T. LDT >= max(1,N).
059: *
060: *  VL      (input) COMPLEX array, dimension (LDVL,M)
061: *          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
062: *          (or of any Q*T*Q**H with Q unitary), corresponding to the
063: *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
064: *          must be stored in consecutive columns of VL, as returned by
065: *          CHSEIN or CTREVC.
066: *          If JOB = 'V', VL is not referenced.
067: *
068: *  LDVL    (input) INTEGER
069: *          The leading dimension of the array VL.
070: *          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
071: *
072: *  VR      (input) COMPLEX array, dimension (LDVR,M)
073: *          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
074: *          (or of any Q*T*Q**H with Q unitary), corresponding to the
075: *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
076: *          must be stored in consecutive columns of VR, as returned by
077: *          CHSEIN or CTREVC.
078: *          If JOB = 'V', VR is not referenced.
079: *
080: *  LDVR    (input) INTEGER
081: *          The leading dimension of the array VR.
082: *          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
083: *
084: *  S       (output) REAL array, dimension (MM)
085: *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
086: *          selected eigenvalues, stored in consecutive elements of the
087: *          array. Thus S(j), SEP(j), and the j-th columns of VL and VR
088: *          all correspond to the same eigenpair (but not in general the
089: *          j-th eigenpair, unless all eigenpairs are selected).
090: *          If JOB = 'V', S is not referenced.
091: *
092: *  SEP     (output) REAL array, dimension (MM)
093: *          If JOB = 'V' or 'B', the estimated reciprocal condition
094: *          numbers of the selected eigenvectors, stored in consecutive
095: *          elements of the array.
096: *          If JOB = 'E', SEP is not referenced.
097: *
098: *  MM      (input) INTEGER
099: *          The number of elements in the arrays S (if JOB = 'E' or 'B')
100: *           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
101: *
102: *  M       (output) INTEGER
103: *          The number of elements of the arrays S and/or SEP actually
104: *          used to store the estimated condition numbers.
105: *          If HOWMNY = 'A', M is set to N.
106: *
107: *  WORK    (workspace) COMPLEX array, dimension (LDWORK,N+6)
108: *          If JOB = 'E', WORK is not referenced.
109: *
110: *  LDWORK  (input) INTEGER
111: *          The leading dimension of the array WORK.
112: *          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
113: *
114: *  RWORK   (workspace) REAL array, dimension (N)
115: *          If JOB = 'E', RWORK is not referenced.
116: *
117: *  INFO    (output) INTEGER
118: *          = 0: successful exit
119: *          < 0: if INFO = -i, the i-th argument had an illegal value
120: *
121: *  Further Details
122: *  ===============
123: *
124: *  The reciprocal of the condition number of an eigenvalue lambda is
125: *  defined as
126: *
127: *          S(lambda) = |v'*u| / (norm(u)*norm(v))
128: *
129: *  where u and v are the right and left eigenvectors of T corresponding
130: *  to lambda; v' denotes the conjugate transpose of v, and norm(u)
131: *  denotes the Euclidean norm. These reciprocal condition numbers always
132: *  lie between zero (very badly conditioned) and one (very well
133: *  conditioned). If n = 1, S(lambda) is defined to be 1.
134: *
135: *  An approximate error bound for a computed eigenvalue W(i) is given by
136: *
137: *                      EPS * norm(T) / S(i)
138: *
139: *  where EPS is the machine precision.
140: *
141: *  The reciprocal of the condition number of the right eigenvector u
142: *  corresponding to lambda is defined as follows. Suppose
143: *
144: *              T = ( lambda  c  )
145: *                  (   0    T22 )
146: *
147: *  Then the reciprocal condition number is
148: *
149: *          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
150: *
151: *  where sigma-min denotes the smallest singular value. We approximate
152: *  the smallest singular value by the reciprocal of an estimate of the
153: *  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
154: *  defined to be abs(T(1,1)).
155: *
156: *  An approximate error bound for a computed right eigenvector VR(i)
157: *  is given by
158: *
159: *                      EPS * norm(T) / SEP(i)
160: *
161: *  =====================================================================
162: *
163: *     .. Parameters ..
164:       REAL               ZERO, ONE
165:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0+0 )
166: *     ..
167: *     .. Local Scalars ..
168:       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
169:       CHARACTER          NORMIN
170:       INTEGER            I, IERR, IX, J, K, KASE, KS
171:       REAL               BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
172:      $                   XNORM
173:       COMPLEX            CDUM, PROD
174: *     ..
175: *     .. Local Arrays ..
176:       INTEGER            ISAVE( 3 )
177:       COMPLEX            DUMMY( 1 )
178: *     ..
179: *     .. External Functions ..
180:       LOGICAL            LSAME
181:       INTEGER            ICAMAX
182:       REAL               SCNRM2, SLAMCH
183:       COMPLEX            CDOTC
184:       EXTERNAL           LSAME, ICAMAX, SCNRM2, SLAMCH, CDOTC
185: *     ..
186: *     .. External Subroutines ..
187:       EXTERNAL           CLACN2, CLACPY, CLATRS, CSRSCL, CTREXC, SLABAD,
188:      $                   XERBLA
189: *     ..
190: *     .. Intrinsic Functions ..
191:       INTRINSIC          ABS, AIMAG, MAX, REAL
192: *     ..
193: *     .. Statement Functions ..
194:       REAL               CABS1
195: *     ..
196: *     .. Statement Function definitions ..
197:       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
198: *     ..
199: *     .. Executable Statements ..
200: *
201: *     Decode and test the input parameters
202: *
203:       WANTBH = LSAME( JOB, 'B' )
204:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
205:       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
206: *
207:       SOMCON = LSAME( HOWMNY, 'S' )
208: *
209: *     Set M to the number of eigenpairs for which condition numbers are
210: *     to be computed.
211: *
212:       IF( SOMCON ) THEN
213:          M = 0
214:          DO 10 J = 1, N
215:             IF( SELECT( J ) )
216:      $         M = M + 1
217:    10    CONTINUE
218:       ELSE
219:          M = N
220:       END IF
221: *
222:       INFO = 0
223:       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
224:          INFO = -1
225:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
226:          INFO = -2
227:       ELSE IF( N.LT.0 ) THEN
228:          INFO = -4
229:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
230:          INFO = -6
231:       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
232:          INFO = -8
233:       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
234:          INFO = -10
235:       ELSE IF( MM.LT.M ) THEN
236:          INFO = -13
237:       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
238:          INFO = -16
239:       END IF
240:       IF( INFO.NE.0 ) THEN
241:          CALL XERBLA( 'CTRSNA', -INFO )
242:          RETURN
243:       END IF
244: *
245: *     Quick return if possible
246: *
247:       IF( N.EQ.0 )
248:      $   RETURN
249: *
250:       IF( N.EQ.1 ) THEN
251:          IF( SOMCON ) THEN
252:             IF( .NOT.SELECT( 1 ) )
253:      $         RETURN
254:          END IF
255:          IF( WANTS )
256:      $      S( 1 ) = ONE
257:          IF( WANTSP )
258:      $      SEP( 1 ) = ABS( T( 1, 1 ) )
259:          RETURN
260:       END IF
261: *
262: *     Get machine constants
263: *
264:       EPS = SLAMCH( 'P' )
265:       SMLNUM = SLAMCH( 'S' ) / EPS
266:       BIGNUM = ONE / SMLNUM
267:       CALL SLABAD( SMLNUM, BIGNUM )
268: *
269:       KS = 1
270:       DO 50 K = 1, N
271: *
272:          IF( SOMCON ) THEN
273:             IF( .NOT.SELECT( K ) )
274:      $         GO TO 50
275:          END IF
276: *
277:          IF( WANTS ) THEN
278: *
279: *           Compute the reciprocal condition number of the k-th
280: *           eigenvalue.
281: *
282:             PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
283:             RNRM = SCNRM2( N, VR( 1, KS ), 1 )
284:             LNRM = SCNRM2( N, VL( 1, KS ), 1 )
285:             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
286: *
287:          END IF
288: *
289:          IF( WANTSP ) THEN
290: *
291: *           Estimate the reciprocal condition number of the k-th
292: *           eigenvector.
293: *
294: *           Copy the matrix T to the array WORK and swap the k-th
295: *           diagonal element to the (1,1) position.
296: *
297:             CALL CLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
298:             CALL CTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
299: *
300: *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
301: *
302:             DO 20 I = 2, N
303:                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
304:    20       CONTINUE
305: *
306: *           Estimate a lower bound for the 1-norm of inv(C'). The 1st
307: *           and (N+1)th columns of WORK are used to store work vectors.
308: *
309:             SEP( KS ) = ZERO
310:             EST = ZERO
311:             KASE = 0
312:             NORMIN = 'N'
313:    30       CONTINUE
314:             CALL CLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
315: *
316:             IF( KASE.NE.0 ) THEN
317:                IF( KASE.EQ.1 ) THEN
318: *
319: *                 Solve C'*x = scale*b
320: *
321:                   CALL CLATRS( 'Upper', 'Conjugate transpose',
322:      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
323:      $                         LDWORK, WORK, SCALE, RWORK, IERR )
324:                ELSE
325: *
326: *                 Solve C*x = scale*b
327: *
328:                   CALL CLATRS( 'Upper', 'No transpose', 'Nonunit',
329:      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
330:      $                         SCALE, RWORK, IERR )
331:                END IF
332:                NORMIN = 'Y'
333:                IF( SCALE.NE.ONE ) THEN
334: *
335: *                 Multiply by 1/SCALE if doing so will not cause
336: *                 overflow.
337: *
338:                   IX = ICAMAX( N-1, WORK, 1 )
339:                   XNORM = CABS1( WORK( IX, 1 ) )
340:                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
341:      $               GO TO 40
342:                   CALL CSRSCL( N, SCALE, WORK, 1 )
343:                END IF
344:                GO TO 30
345:             END IF
346: *
347:             SEP( KS ) = ONE / MAX( EST, SMLNUM )
348:          END IF
349: *
350:    40    CONTINUE
351:          KS = KS + 1
352:    50 CONTINUE
353:       RETURN
354: *
355: *     End of CTRSNA
356: *
357:       END
358: