1 SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
2 $ SEP, WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
10 *
11 * .. Scalar Arguments ..
12 CHARACTER COMPQ, JOB
13 INTEGER INFO, LDQ, LDT, LWORK, M, N
14 DOUBLE PRECISION S, SEP
15 * ..
16 * .. Array Arguments ..
17 LOGICAL SELECT( * )
18 COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * ZTRSEN reorders the Schur factorization of a complex matrix
25 * A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
26 * the leading positions on the diagonal of the upper triangular matrix
27 * T, and the leading columns of Q form an orthonormal basis of the
28 * corresponding right invariant subspace.
29 *
30 * Optionally the routine computes the reciprocal condition numbers of
31 * the cluster of eigenvalues and/or the invariant subspace.
32 *
33 * Arguments
34 * =========
35 *
36 * JOB (input) CHARACTER*1
37 * Specifies whether condition numbers are required for the
38 * cluster of eigenvalues (S) or the invariant subspace (SEP):
39 * = 'N': none;
40 * = 'E': for eigenvalues only (S);
41 * = 'V': for invariant subspace only (SEP);
42 * = 'B': for both eigenvalues and invariant subspace (S and
43 * SEP).
44 *
45 * COMPQ (input) CHARACTER*1
46 * = 'V': update the matrix Q of Schur vectors;
47 * = 'N': do not update Q.
48 *
49 * SELECT (input) LOGICAL array, dimension (N)
50 * SELECT specifies the eigenvalues in the selected cluster. To
51 * select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
52 *
53 * N (input) INTEGER
54 * The order of the matrix T. N >= 0.
55 *
56 * T (input/output) COMPLEX*16 array, dimension (LDT,N)
57 * On entry, the upper triangular matrix T.
58 * On exit, T is overwritten by the reordered matrix T, with the
59 * selected eigenvalues as the leading diagonal elements.
60 *
61 * LDT (input) INTEGER
62 * The leading dimension of the array T. LDT >= max(1,N).
63 *
64 * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
65 * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
66 * On exit, if COMPQ = 'V', Q has been postmultiplied by the
67 * unitary transformation matrix which reorders T; the leading M
68 * columns of Q form an orthonormal basis for the specified
69 * invariant subspace.
70 * If COMPQ = 'N', Q is not referenced.
71 *
72 * LDQ (input) INTEGER
73 * The leading dimension of the array Q.
74 * LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
75 *
76 * W (output) COMPLEX*16 array, dimension (N)
77 * The reordered eigenvalues of T, in the same order as they
78 * appear on the diagonal of T.
79 *
80 * M (output) INTEGER
81 * The dimension of the specified invariant subspace.
82 * 0 <= M <= N.
83 *
84 * S (output) DOUBLE PRECISION
85 * If JOB = 'E' or 'B', S is a lower bound on the reciprocal
86 * condition number for the selected cluster of eigenvalues.
87 * S cannot underestimate the true reciprocal condition number
88 * by more than a factor of sqrt(N). If M = 0 or N, S = 1.
89 * If JOB = 'N' or 'V', S is not referenced.
90 *
91 * SEP (output) DOUBLE PRECISION
92 * If JOB = 'V' or 'B', SEP is the estimated reciprocal
93 * condition number of the specified invariant subspace. If
94 * M = 0 or N, SEP = norm(T).
95 * If JOB = 'N' or 'E', SEP is not referenced.
96 *
97 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
98 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
99 *
100 * LWORK (input) INTEGER
101 * The dimension of the array WORK.
102 * If JOB = 'N', LWORK >= 1;
103 * if JOB = 'E', LWORK = max(1,M*(N-M));
104 * if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
105 *
106 * If LWORK = -1, then a workspace query is assumed; the routine
107 * only calculates the optimal size of the WORK array, returns
108 * this value as the first entry of the WORK array, and no error
109 * message related to LWORK is issued by XERBLA.
110 *
111 * INFO (output) INTEGER
112 * = 0: successful exit
113 * < 0: if INFO = -i, the i-th argument had an illegal value
114 *
115 * Further Details
116 * ===============
117 *
118 * ZTRSEN first collects the selected eigenvalues by computing a unitary
119 * transformation Z to move them to the top left corner of T. In other
120 * words, the selected eigenvalues are the eigenvalues of T11 in:
121 *
122 * Z**H * T * Z = ( T11 T12 ) n1
123 * ( 0 T22 ) n2
124 * n1 n2
125 *
126 * where N = n1+n2. The first
127 * n1 columns of Z span the specified invariant subspace of T.
128 *
129 * If T has been obtained from the Schur factorization of a matrix
130 * A = Q*T*Q**H, then the reordered Schur factorization of A is given by
131 * A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
132 * corresponding invariant subspace of A.
133 *
134 * The reciprocal condition number of the average of the eigenvalues of
135 * T11 may be returned in S. S lies between 0 (very badly conditioned)
136 * and 1 (very well conditioned). It is computed as follows. First we
137 * compute R so that
138 *
139 * P = ( I R ) n1
140 * ( 0 0 ) n2
141 * n1 n2
142 *
143 * is the projector on the invariant subspace associated with T11.
144 * R is the solution of the Sylvester equation:
145 *
146 * T11*R - R*T22 = T12.
147 *
148 * Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
149 * the two-norm of M. Then S is computed as the lower bound
150 *
151 * (1 + F-norm(R)**2)**(-1/2)
152 *
153 * on the reciprocal of 2-norm(P), the true reciprocal condition number.
154 * S cannot underestimate 1 / 2-norm(P) by more than a factor of
155 * sqrt(N).
156 *
157 * An approximate error bound for the computed average of the
158 * eigenvalues of T11 is
159 *
160 * EPS * norm(T) / S
161 *
162 * where EPS is the machine precision.
163 *
164 * The reciprocal condition number of the right invariant subspace
165 * spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
166 * SEP is defined as the separation of T11 and T22:
167 *
168 * sep( T11, T22 ) = sigma-min( C )
169 *
170 * where sigma-min(C) is the smallest singular value of the
171 * n1*n2-by-n1*n2 matrix
172 *
173 * C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
174 *
175 * I(m) is an m by m identity matrix, and kprod denotes the Kronecker
176 * product. We estimate sigma-min(C) by the reciprocal of an estimate of
177 * the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
178 * cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
179 *
180 * When SEP is small, small changes in T can cause large changes in
181 * the invariant subspace. An approximate bound on the maximum angular
182 * error in the computed right invariant subspace is
183 *
184 * EPS * norm(T) / SEP
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189 DOUBLE PRECISION ZERO, ONE
190 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
191 * ..
192 * .. Local Scalars ..
193 LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
194 INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
195 DOUBLE PRECISION EST, RNORM, SCALE
196 * ..
197 * .. Local Arrays ..
198 INTEGER ISAVE( 3 )
199 DOUBLE PRECISION RWORK( 1 )
200 * ..
201 * .. External Functions ..
202 LOGICAL LSAME
203 DOUBLE PRECISION ZLANGE
204 EXTERNAL LSAME, ZLANGE
205 * ..
206 * .. External Subroutines ..
207 EXTERNAL XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
208 * ..
209 * .. Intrinsic Functions ..
210 INTRINSIC MAX, SQRT
211 * ..
212 * .. Executable Statements ..
213 *
214 * Decode and test the input parameters.
215 *
216 WANTBH = LSAME( JOB, 'B' )
217 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
218 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
219 WANTQ = LSAME( COMPQ, 'V' )
220 *
221 * Set M to the number of selected eigenvalues.
222 *
223 M = 0
224 DO 10 K = 1, N
225 IF( SELECT( K ) )
226 $ M = M + 1
227 10 CONTINUE
228 *
229 N1 = M
230 N2 = N - M
231 NN = N1*N2
232 *
233 INFO = 0
234 LQUERY = ( LWORK.EQ.-1 )
235 *
236 IF( WANTSP ) THEN
237 LWMIN = MAX( 1, 2*NN )
238 ELSE IF( LSAME( JOB, 'N' ) ) THEN
239 LWMIN = 1
240 ELSE IF( LSAME( JOB, 'E' ) ) THEN
241 LWMIN = MAX( 1, NN )
242 END IF
243 *
244 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
245 $ THEN
246 INFO = -1
247 ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
248 INFO = -2
249 ELSE IF( N.LT.0 ) THEN
250 INFO = -4
251 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
252 INFO = -6
253 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
254 INFO = -8
255 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
256 INFO = -14
257 END IF
258 *
259 IF( INFO.EQ.0 ) THEN
260 WORK( 1 ) = LWMIN
261 END IF
262 *
263 IF( INFO.NE.0 ) THEN
264 CALL XERBLA( 'ZTRSEN', -INFO )
265 RETURN
266 ELSE IF( LQUERY ) THEN
267 RETURN
268 END IF
269 *
270 * Quick return if possible
271 *
272 IF( M.EQ.N .OR. M.EQ.0 ) THEN
273 IF( WANTS )
274 $ S = ONE
275 IF( WANTSP )
276 $ SEP = ZLANGE( '1', N, N, T, LDT, RWORK )
277 GO TO 40
278 END IF
279 *
280 * Collect the selected eigenvalues at the top left corner of T.
281 *
282 KS = 0
283 DO 20 K = 1, N
284 IF( SELECT( K ) ) THEN
285 KS = KS + 1
286 *
287 * Swap the K-th eigenvalue to position KS.
288 *
289 IF( K.NE.KS )
290 $ CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
291 END IF
292 20 CONTINUE
293 *
294 IF( WANTS ) THEN
295 *
296 * Solve the Sylvester equation for R:
297 *
298 * T11*R - R*T22 = scale*T12
299 *
300 CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
301 CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
302 $ LDT, WORK, N1, SCALE, IERR )
303 *
304 * Estimate the reciprocal of the condition number of the cluster
305 * of eigenvalues.
306 *
307 RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK )
308 IF( RNORM.EQ.ZERO ) THEN
309 S = ONE
310 ELSE
311 S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
312 $ SQRT( RNORM ) )
313 END IF
314 END IF
315 *
316 IF( WANTSP ) THEN
317 *
318 * Estimate sep(T11,T22).
319 *
320 EST = ZERO
321 KASE = 0
322 30 CONTINUE
323 CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
324 IF( KASE.NE.0 ) THEN
325 IF( KASE.EQ.1 ) THEN
326 *
327 * Solve T11*R - R*T22 = scale*X.
328 *
329 CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
330 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
331 $ IERR )
332 ELSE
333 *
334 * Solve T11**H*R - R*T22**H = scale*X.
335 *
336 CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
337 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
338 $ IERR )
339 END IF
340 GO TO 30
341 END IF
342 *
343 SEP = SCALE / EST
344 END IF
345 *
346 40 CONTINUE
347 *
348 * Copy reordered eigenvalues to W.
349 *
350 DO 50 K = 1, N
351 W( K ) = T( K, K )
352 50 CONTINUE
353 *
354 WORK( 1 ) = LWMIN
355 *
356 RETURN
357 *
358 * End of ZTRSEN
359 *
360 END
2 $ SEP, WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
10 *
11 * .. Scalar Arguments ..
12 CHARACTER COMPQ, JOB
13 INTEGER INFO, LDQ, LDT, LWORK, M, N
14 DOUBLE PRECISION S, SEP
15 * ..
16 * .. Array Arguments ..
17 LOGICAL SELECT( * )
18 COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * ZTRSEN reorders the Schur factorization of a complex matrix
25 * A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
26 * the leading positions on the diagonal of the upper triangular matrix
27 * T, and the leading columns of Q form an orthonormal basis of the
28 * corresponding right invariant subspace.
29 *
30 * Optionally the routine computes the reciprocal condition numbers of
31 * the cluster of eigenvalues and/or the invariant subspace.
32 *
33 * Arguments
34 * =========
35 *
36 * JOB (input) CHARACTER*1
37 * Specifies whether condition numbers are required for the
38 * cluster of eigenvalues (S) or the invariant subspace (SEP):
39 * = 'N': none;
40 * = 'E': for eigenvalues only (S);
41 * = 'V': for invariant subspace only (SEP);
42 * = 'B': for both eigenvalues and invariant subspace (S and
43 * SEP).
44 *
45 * COMPQ (input) CHARACTER*1
46 * = 'V': update the matrix Q of Schur vectors;
47 * = 'N': do not update Q.
48 *
49 * SELECT (input) LOGICAL array, dimension (N)
50 * SELECT specifies the eigenvalues in the selected cluster. To
51 * select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
52 *
53 * N (input) INTEGER
54 * The order of the matrix T. N >= 0.
55 *
56 * T (input/output) COMPLEX*16 array, dimension (LDT,N)
57 * On entry, the upper triangular matrix T.
58 * On exit, T is overwritten by the reordered matrix T, with the
59 * selected eigenvalues as the leading diagonal elements.
60 *
61 * LDT (input) INTEGER
62 * The leading dimension of the array T. LDT >= max(1,N).
63 *
64 * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
65 * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
66 * On exit, if COMPQ = 'V', Q has been postmultiplied by the
67 * unitary transformation matrix which reorders T; the leading M
68 * columns of Q form an orthonormal basis for the specified
69 * invariant subspace.
70 * If COMPQ = 'N', Q is not referenced.
71 *
72 * LDQ (input) INTEGER
73 * The leading dimension of the array Q.
74 * LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
75 *
76 * W (output) COMPLEX*16 array, dimension (N)
77 * The reordered eigenvalues of T, in the same order as they
78 * appear on the diagonal of T.
79 *
80 * M (output) INTEGER
81 * The dimension of the specified invariant subspace.
82 * 0 <= M <= N.
83 *
84 * S (output) DOUBLE PRECISION
85 * If JOB = 'E' or 'B', S is a lower bound on the reciprocal
86 * condition number for the selected cluster of eigenvalues.
87 * S cannot underestimate the true reciprocal condition number
88 * by more than a factor of sqrt(N). If M = 0 or N, S = 1.
89 * If JOB = 'N' or 'V', S is not referenced.
90 *
91 * SEP (output) DOUBLE PRECISION
92 * If JOB = 'V' or 'B', SEP is the estimated reciprocal
93 * condition number of the specified invariant subspace. If
94 * M = 0 or N, SEP = norm(T).
95 * If JOB = 'N' or 'E', SEP is not referenced.
96 *
97 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
98 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
99 *
100 * LWORK (input) INTEGER
101 * The dimension of the array WORK.
102 * If JOB = 'N', LWORK >= 1;
103 * if JOB = 'E', LWORK = max(1,M*(N-M));
104 * if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
105 *
106 * If LWORK = -1, then a workspace query is assumed; the routine
107 * only calculates the optimal size of the WORK array, returns
108 * this value as the first entry of the WORK array, and no error
109 * message related to LWORK is issued by XERBLA.
110 *
111 * INFO (output) INTEGER
112 * = 0: successful exit
113 * < 0: if INFO = -i, the i-th argument had an illegal value
114 *
115 * Further Details
116 * ===============
117 *
118 * ZTRSEN first collects the selected eigenvalues by computing a unitary
119 * transformation Z to move them to the top left corner of T. In other
120 * words, the selected eigenvalues are the eigenvalues of T11 in:
121 *
122 * Z**H * T * Z = ( T11 T12 ) n1
123 * ( 0 T22 ) n2
124 * n1 n2
125 *
126 * where N = n1+n2. The first
127 * n1 columns of Z span the specified invariant subspace of T.
128 *
129 * If T has been obtained from the Schur factorization of a matrix
130 * A = Q*T*Q**H, then the reordered Schur factorization of A is given by
131 * A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
132 * corresponding invariant subspace of A.
133 *
134 * The reciprocal condition number of the average of the eigenvalues of
135 * T11 may be returned in S. S lies between 0 (very badly conditioned)
136 * and 1 (very well conditioned). It is computed as follows. First we
137 * compute R so that
138 *
139 * P = ( I R ) n1
140 * ( 0 0 ) n2
141 * n1 n2
142 *
143 * is the projector on the invariant subspace associated with T11.
144 * R is the solution of the Sylvester equation:
145 *
146 * T11*R - R*T22 = T12.
147 *
148 * Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
149 * the two-norm of M. Then S is computed as the lower bound
150 *
151 * (1 + F-norm(R)**2)**(-1/2)
152 *
153 * on the reciprocal of 2-norm(P), the true reciprocal condition number.
154 * S cannot underestimate 1 / 2-norm(P) by more than a factor of
155 * sqrt(N).
156 *
157 * An approximate error bound for the computed average of the
158 * eigenvalues of T11 is
159 *
160 * EPS * norm(T) / S
161 *
162 * where EPS is the machine precision.
163 *
164 * The reciprocal condition number of the right invariant subspace
165 * spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
166 * SEP is defined as the separation of T11 and T22:
167 *
168 * sep( T11, T22 ) = sigma-min( C )
169 *
170 * where sigma-min(C) is the smallest singular value of the
171 * n1*n2-by-n1*n2 matrix
172 *
173 * C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
174 *
175 * I(m) is an m by m identity matrix, and kprod denotes the Kronecker
176 * product. We estimate sigma-min(C) by the reciprocal of an estimate of
177 * the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
178 * cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
179 *
180 * When SEP is small, small changes in T can cause large changes in
181 * the invariant subspace. An approximate bound on the maximum angular
182 * error in the computed right invariant subspace is
183 *
184 * EPS * norm(T) / SEP
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189 DOUBLE PRECISION ZERO, ONE
190 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
191 * ..
192 * .. Local Scalars ..
193 LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
194 INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
195 DOUBLE PRECISION EST, RNORM, SCALE
196 * ..
197 * .. Local Arrays ..
198 INTEGER ISAVE( 3 )
199 DOUBLE PRECISION RWORK( 1 )
200 * ..
201 * .. External Functions ..
202 LOGICAL LSAME
203 DOUBLE PRECISION ZLANGE
204 EXTERNAL LSAME, ZLANGE
205 * ..
206 * .. External Subroutines ..
207 EXTERNAL XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
208 * ..
209 * .. Intrinsic Functions ..
210 INTRINSIC MAX, SQRT
211 * ..
212 * .. Executable Statements ..
213 *
214 * Decode and test the input parameters.
215 *
216 WANTBH = LSAME( JOB, 'B' )
217 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
218 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
219 WANTQ = LSAME( COMPQ, 'V' )
220 *
221 * Set M to the number of selected eigenvalues.
222 *
223 M = 0
224 DO 10 K = 1, N
225 IF( SELECT( K ) )
226 $ M = M + 1
227 10 CONTINUE
228 *
229 N1 = M
230 N2 = N - M
231 NN = N1*N2
232 *
233 INFO = 0
234 LQUERY = ( LWORK.EQ.-1 )
235 *
236 IF( WANTSP ) THEN
237 LWMIN = MAX( 1, 2*NN )
238 ELSE IF( LSAME( JOB, 'N' ) ) THEN
239 LWMIN = 1
240 ELSE IF( LSAME( JOB, 'E' ) ) THEN
241 LWMIN = MAX( 1, NN )
242 END IF
243 *
244 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
245 $ THEN
246 INFO = -1
247 ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
248 INFO = -2
249 ELSE IF( N.LT.0 ) THEN
250 INFO = -4
251 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
252 INFO = -6
253 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
254 INFO = -8
255 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
256 INFO = -14
257 END IF
258 *
259 IF( INFO.EQ.0 ) THEN
260 WORK( 1 ) = LWMIN
261 END IF
262 *
263 IF( INFO.NE.0 ) THEN
264 CALL XERBLA( 'ZTRSEN', -INFO )
265 RETURN
266 ELSE IF( LQUERY ) THEN
267 RETURN
268 END IF
269 *
270 * Quick return if possible
271 *
272 IF( M.EQ.N .OR. M.EQ.0 ) THEN
273 IF( WANTS )
274 $ S = ONE
275 IF( WANTSP )
276 $ SEP = ZLANGE( '1', N, N, T, LDT, RWORK )
277 GO TO 40
278 END IF
279 *
280 * Collect the selected eigenvalues at the top left corner of T.
281 *
282 KS = 0
283 DO 20 K = 1, N
284 IF( SELECT( K ) ) THEN
285 KS = KS + 1
286 *
287 * Swap the K-th eigenvalue to position KS.
288 *
289 IF( K.NE.KS )
290 $ CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
291 END IF
292 20 CONTINUE
293 *
294 IF( WANTS ) THEN
295 *
296 * Solve the Sylvester equation for R:
297 *
298 * T11*R - R*T22 = scale*T12
299 *
300 CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
301 CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
302 $ LDT, WORK, N1, SCALE, IERR )
303 *
304 * Estimate the reciprocal of the condition number of the cluster
305 * of eigenvalues.
306 *
307 RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK )
308 IF( RNORM.EQ.ZERO ) THEN
309 S = ONE
310 ELSE
311 S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
312 $ SQRT( RNORM ) )
313 END IF
314 END IF
315 *
316 IF( WANTSP ) THEN
317 *
318 * Estimate sep(T11,T22).
319 *
320 EST = ZERO
321 KASE = 0
322 30 CONTINUE
323 CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
324 IF( KASE.NE.0 ) THEN
325 IF( KASE.EQ.1 ) THEN
326 *
327 * Solve T11*R - R*T22 = scale*X.
328 *
329 CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
330 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
331 $ IERR )
332 ELSE
333 *
334 * Solve T11**H*R - R*T22**H = scale*X.
335 *
336 CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
337 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
338 $ IERR )
339 END IF
340 GO TO 30
341 END IF
342 *
343 SEP = SCALE / EST
344 END IF
345 *
346 40 CONTINUE
347 *
348 * Copy reordered eigenvalues to W.
349 *
350 DO 50 K = 1, N
351 W( K ) = T( K, K )
352 50 CONTINUE
353 *
354 WORK( 1 ) = LWMIN
355 *
356 RETURN
357 *
358 * End of ZTRSEN
359 *
360 END