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