001:       SUBROUTINE ZTRSNA( 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 ZLACN2 in place of ZLACON, 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:       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
018:       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
019:      $                   WORK( LDWORK, * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  ZTRSNA 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*16 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*16 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: *          ZHSEIN or ZTREVC.
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*16 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: *          ZHSEIN or ZTREVC.
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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*16 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) DOUBLE PRECISION 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:       DOUBLE PRECISION   ZERO, ONE
164:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
165: *     ..
166: *     .. Local Scalars ..
167:       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
168:       CHARACTER          NORMIN
169:       INTEGER            I, IERR, IX, J, K, KASE, KS
170:       DOUBLE PRECISION   BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
171:      $                   XNORM
172:       COMPLEX*16         CDUM, PROD
173: *     ..
174: *     .. Local Arrays ..
175:       INTEGER            ISAVE( 3 )
176:       COMPLEX*16         DUMMY( 1 )
177: *     ..
178: *     .. External Functions ..
179:       LOGICAL            LSAME
180:       INTEGER            IZAMAX
181:       DOUBLE PRECISION   DLAMCH, DZNRM2
182:       COMPLEX*16         ZDOTC
183:       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
184: *     ..
185: *     .. External Subroutines ..
186:       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
187: *     ..
188: *     .. Intrinsic Functions ..
189:       INTRINSIC          ABS, DBLE, DIMAG, MAX
190: *     ..
191: *     .. Statement Functions ..
192:       DOUBLE PRECISION   CABS1
193: *     ..
194: *     .. Statement Function definitions ..
195:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
196: *     ..
197: *     .. Executable Statements ..
198: *
199: *     Decode and test the input parameters
200: *
201:       WANTBH = LSAME( JOB, 'B' )
202:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
203:       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
204: *
205:       SOMCON = LSAME( HOWMNY, 'S' )
206: *
207: *     Set M to the number of eigenpairs for which condition numbers are
208: *     to be computed.
209: *
210:       IF( SOMCON ) THEN
211:          M = 0
212:          DO 10 J = 1, N
213:             IF( SELECT( J ) )
214:      $         M = M + 1
215:    10    CONTINUE
216:       ELSE
217:          M = N
218:       END IF
219: *
220:       INFO = 0
221:       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
222:          INFO = -1
223:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
224:          INFO = -2
225:       ELSE IF( N.LT.0 ) THEN
226:          INFO = -4
227:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
228:          INFO = -6
229:       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
230:          INFO = -8
231:       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
232:          INFO = -10
233:       ELSE IF( MM.LT.M ) THEN
234:          INFO = -13
235:       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
236:          INFO = -16
237:       END IF
238:       IF( INFO.NE.0 ) THEN
239:          CALL XERBLA( 'ZTRSNA', -INFO )
240:          RETURN
241:       END IF
242: *
243: *     Quick return if possible
244: *
245:       IF( N.EQ.0 )
246:      $   RETURN
247: *
248:       IF( N.EQ.1 ) THEN
249:          IF( SOMCON ) THEN
250:             IF( .NOT.SELECT( 1 ) )
251:      $         RETURN
252:          END IF
253:          IF( WANTS )
254:      $      S( 1 ) = ONE
255:          IF( WANTSP )
256:      $      SEP( 1 ) = ABS( T( 1, 1 ) )
257:          RETURN
258:       END IF
259: *
260: *     Get machine constants
261: *
262:       EPS = DLAMCH( 'P' )
263:       SMLNUM = DLAMCH( 'S' ) / EPS
264:       BIGNUM = ONE / SMLNUM
265:       CALL DLABAD( SMLNUM, BIGNUM )
266: *
267:       KS = 1
268:       DO 50 K = 1, N
269: *
270:          IF( SOMCON ) THEN
271:             IF( .NOT.SELECT( K ) )
272:      $         GO TO 50
273:          END IF
274: *
275:          IF( WANTS ) THEN
276: *
277: *           Compute the reciprocal condition number of the k-th
278: *           eigenvalue.
279: *
280:             PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
281:             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
282:             LNRM = DZNRM2( N, VL( 1, KS ), 1 )
283:             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
284: *
285:          END IF
286: *
287:          IF( WANTSP ) THEN
288: *
289: *           Estimate the reciprocal condition number of the k-th
290: *           eigenvector.
291: *
292: *           Copy the matrix T to the array WORK and swap the k-th
293: *           diagonal element to the (1,1) position.
294: *
295:             CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
296:             CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
297: *
298: *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
299: *
300:             DO 20 I = 2, N
301:                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
302:    20       CONTINUE
303: *
304: *           Estimate a lower bound for the 1-norm of inv(C'). The 1st
305: *           and (N+1)th columns of WORK are used to store work vectors.
306: *
307:             SEP( KS ) = ZERO
308:             EST = ZERO
309:             KASE = 0
310:             NORMIN = 'N'
311:    30       CONTINUE
312:             CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
313: *
314:             IF( KASE.NE.0 ) THEN
315:                IF( KASE.EQ.1 ) THEN
316: *
317: *                 Solve C'*x = scale*b
318: *
319:                   CALL ZLATRS( 'Upper', 'Conjugate transpose',
320:      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
321:      $                         LDWORK, WORK, SCALE, RWORK, IERR )
322:                ELSE
323: *
324: *                 Solve C*x = scale*b
325: *
326:                   CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
327:      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
328:      $                         SCALE, RWORK, IERR )
329:                END IF
330:                NORMIN = 'Y'
331:                IF( SCALE.NE.ONE ) THEN
332: *
333: *                 Multiply by 1/SCALE if doing so will not cause
334: *                 overflow.
335: *
336:                   IX = IZAMAX( N-1, WORK, 1 )
337:                   XNORM = CABS1( WORK( IX, 1 ) )
338:                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
339:      $               GO TO 40
340:                   CALL ZDRSCL( N, SCALE, WORK, 1 )
341:                END IF
342:                GO TO 30
343:             END IF
344: *
345:             SEP( KS ) = ONE / MAX( EST, SMLNUM )
346:          END IF
347: *
348:    40    CONTINUE
349:          KS = KS + 1
350:    50 CONTINUE
351:       RETURN
352: *
353: *     End of ZTRSNA
354: *
355:       END
356: