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: *  -- 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 ZLACN2 in place of ZLACON, 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:       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
019:       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
020:      $                   WORK( LDWORK, * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  ZTRSNA 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*16 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*16 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: *          ZHSEIN or ZTREVC.
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*16 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: *          ZHSEIN or ZTREVC.
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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*16 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) DOUBLE PRECISION 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:       DOUBLE PRECISION   ZERO, ONE
165:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
166: *     ..
167: *     .. Local Scalars ..
168:       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
169:       CHARACTER          NORMIN
170:       INTEGER            I, IERR, IX, J, K, KASE, KS
171:       DOUBLE PRECISION   BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
172:      $                   XNORM
173:       COMPLEX*16         CDUM, PROD
174: *     ..
175: *     .. Local Arrays ..
176:       INTEGER            ISAVE( 3 )
177:       COMPLEX*16         DUMMY( 1 )
178: *     ..
179: *     .. External Functions ..
180:       LOGICAL            LSAME
181:       INTEGER            IZAMAX
182:       DOUBLE PRECISION   DLAMCH, DZNRM2
183:       COMPLEX*16         ZDOTC
184:       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
185: *     ..
186: *     .. External Subroutines ..
187:       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
188: *     ..
189: *     .. Intrinsic Functions ..
190:       INTRINSIC          ABS, DBLE, DIMAG, MAX
191: *     ..
192: *     .. Statement Functions ..
193:       DOUBLE PRECISION   CABS1
194: *     ..
195: *     .. Statement Function definitions ..
196:       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( 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( 'ZTRSNA', -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 = DLAMCH( 'P' )
264:       SMLNUM = DLAMCH( 'S' ) / EPS
265:       BIGNUM = ONE / SMLNUM
266:       CALL DLABAD( 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 = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
282:             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
283:             LNRM = DZNRM2( 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 ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
297:             CALL ZTREXC( '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 ZLACN2( 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 ZLATRS( '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 ZLATRS( '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 = IZAMAX( 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 ZDRSCL( 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 ZTRSNA
355: *
356:       END
357: