1 SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
2 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
3 $ INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
11 *
12 * .. Scalar Arguments ..
13 CHARACTER HOWMNY, JOB
14 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
15 * ..
16 * .. Array Arguments ..
17 LOGICAL SELECT( * )
18 DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
19 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
20 $ WORK( LDWORK, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZTRSNA estimates reciprocal condition numbers for specified
27 * eigenvalues and/or right eigenvectors of a complex upper triangular
28 * matrix T (or of any matrix Q*T*Q**H with Q unitary).
29 *
30 * Arguments
31 * =========
32 *
33 * JOB (input) CHARACTER*1
34 * Specifies whether condition numbers are required for
35 * eigenvalues (S) or eigenvectors (SEP):
36 * = 'E': for eigenvalues only (S);
37 * = 'V': for eigenvectors only (SEP);
38 * = 'B': for both eigenvalues and eigenvectors (S and SEP).
39 *
40 * HOWMNY (input) CHARACTER*1
41 * = 'A': compute condition numbers for all eigenpairs;
42 * = 'S': compute condition numbers for selected eigenpairs
43 * specified by the array SELECT.
44 *
45 * SELECT (input) LOGICAL array, dimension (N)
46 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
47 * condition numbers are required. To select condition numbers
48 * for the j-th eigenpair, SELECT(j) must be set to .TRUE..
49 * If HOWMNY = 'A', SELECT is not referenced.
50 *
51 * N (input) INTEGER
52 * The order of the matrix T. N >= 0.
53 *
54 * T (input) COMPLEX*16 array, dimension (LDT,N)
55 * The upper triangular matrix T.
56 *
57 * LDT (input) INTEGER
58 * The leading dimension of the array T. LDT >= max(1,N).
59 *
60 * VL (input) COMPLEX*16 array, dimension (LDVL,M)
61 * If JOB = 'E' or 'B', VL must contain left eigenvectors of T
62 * (or of any Q*T*Q**H with Q unitary), corresponding to the
63 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
64 * must be stored in consecutive columns of VL, as returned by
65 * ZHSEIN or ZTREVC.
66 * If JOB = 'V', VL is not referenced.
67 *
68 * LDVL (input) INTEGER
69 * The leading dimension of the array VL.
70 * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
71 *
72 * VR (input) COMPLEX*16 array, dimension (LDVR,M)
73 * If JOB = 'E' or 'B', VR must contain right eigenvectors of T
74 * (or of any Q*T*Q**H with Q unitary), corresponding to the
75 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
76 * must be stored in consecutive columns of VR, as returned by
77 * ZHSEIN or ZTREVC.
78 * If JOB = 'V', VR is not referenced.
79 *
80 * LDVR (input) INTEGER
81 * The leading dimension of the array VR.
82 * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
83 *
84 * S (output) DOUBLE PRECISION array, dimension (MM)
85 * If JOB = 'E' or 'B', the reciprocal condition numbers of the
86 * selected eigenvalues, stored in consecutive elements of the
87 * array. Thus S(j), SEP(j), and the j-th columns of VL and VR
88 * all correspond to the same eigenpair (but not in general the
89 * j-th eigenpair, unless all eigenpairs are selected).
90 * If JOB = 'V', S is not referenced.
91 *
92 * SEP (output) DOUBLE PRECISION array, dimension (MM)
93 * If JOB = 'V' or 'B', the estimated reciprocal condition
94 * numbers of the selected eigenvectors, stored in consecutive
95 * elements of the array.
96 * If JOB = 'E', SEP is not referenced.
97 *
98 * MM (input) INTEGER
99 * 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**H*u| / (norm(u)*norm(v))
128 *
129 * where u and v are the right and left eigenvectors of T corresponding
130 * to lambda; v**H 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**H). 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**H*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
2 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
3 $ INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
11 *
12 * .. Scalar Arguments ..
13 CHARACTER HOWMNY, JOB
14 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
15 * ..
16 * .. Array Arguments ..
17 LOGICAL SELECT( * )
18 DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
19 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
20 $ WORK( LDWORK, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZTRSNA estimates reciprocal condition numbers for specified
27 * eigenvalues and/or right eigenvectors of a complex upper triangular
28 * matrix T (or of any matrix Q*T*Q**H with Q unitary).
29 *
30 * Arguments
31 * =========
32 *
33 * JOB (input) CHARACTER*1
34 * Specifies whether condition numbers are required for
35 * eigenvalues (S) or eigenvectors (SEP):
36 * = 'E': for eigenvalues only (S);
37 * = 'V': for eigenvectors only (SEP);
38 * = 'B': for both eigenvalues and eigenvectors (S and SEP).
39 *
40 * HOWMNY (input) CHARACTER*1
41 * = 'A': compute condition numbers for all eigenpairs;
42 * = 'S': compute condition numbers for selected eigenpairs
43 * specified by the array SELECT.
44 *
45 * SELECT (input) LOGICAL array, dimension (N)
46 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
47 * condition numbers are required. To select condition numbers
48 * for the j-th eigenpair, SELECT(j) must be set to .TRUE..
49 * If HOWMNY = 'A', SELECT is not referenced.
50 *
51 * N (input) INTEGER
52 * The order of the matrix T. N >= 0.
53 *
54 * T (input) COMPLEX*16 array, dimension (LDT,N)
55 * The upper triangular matrix T.
56 *
57 * LDT (input) INTEGER
58 * The leading dimension of the array T. LDT >= max(1,N).
59 *
60 * VL (input) COMPLEX*16 array, dimension (LDVL,M)
61 * If JOB = 'E' or 'B', VL must contain left eigenvectors of T
62 * (or of any Q*T*Q**H with Q unitary), corresponding to the
63 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
64 * must be stored in consecutive columns of VL, as returned by
65 * ZHSEIN or ZTREVC.
66 * If JOB = 'V', VL is not referenced.
67 *
68 * LDVL (input) INTEGER
69 * The leading dimension of the array VL.
70 * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
71 *
72 * VR (input) COMPLEX*16 array, dimension (LDVR,M)
73 * If JOB = 'E' or 'B', VR must contain right eigenvectors of T
74 * (or of any Q*T*Q**H with Q unitary), corresponding to the
75 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors
76 * must be stored in consecutive columns of VR, as returned by
77 * ZHSEIN or ZTREVC.
78 * If JOB = 'V', VR is not referenced.
79 *
80 * LDVR (input) INTEGER
81 * The leading dimension of the array VR.
82 * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
83 *
84 * S (output) DOUBLE PRECISION array, dimension (MM)
85 * If JOB = 'E' or 'B', the reciprocal condition numbers of the
86 * selected eigenvalues, stored in consecutive elements of the
87 * array. Thus S(j), SEP(j), and the j-th columns of VL and VR
88 * all correspond to the same eigenpair (but not in general the
89 * j-th eigenpair, unless all eigenpairs are selected).
90 * If JOB = 'V', S is not referenced.
91 *
92 * SEP (output) DOUBLE PRECISION array, dimension (MM)
93 * If JOB = 'V' or 'B', the estimated reciprocal condition
94 * numbers of the selected eigenvectors, stored in consecutive
95 * elements of the array.
96 * If JOB = 'E', SEP is not referenced.
97 *
98 * MM (input) INTEGER
99 * 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**H*u| / (norm(u)*norm(v))
128 *
129 * where u and v are the right and left eigenvectors of T corresponding
130 * to lambda; v**H 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**H). 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**H*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