LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
ztrsen.f
Go to the documentation of this file.
1 *> \brief \b ZTRSEN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZTRSEN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsen.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsen.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsen.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
22 * SEP, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER COMPQ, JOB
26 * INTEGER INFO, LDQ, LDT, LWORK, M, N
27 * DOUBLE PRECISION S, SEP
28 * ..
29 * .. Array Arguments ..
30 * LOGICAL SELECT( * )
31 * COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZTRSEN reorders the Schur factorization of a complex matrix
41 *> A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
42 *> the leading positions on the diagonal of the upper triangular matrix
43 *> T, and the leading columns of Q form an orthonormal basis of the
44 *> corresponding right invariant subspace.
45 *>
46 *> Optionally the routine computes the reciprocal condition numbers of
47 *> the cluster of eigenvalues and/or the invariant subspace.
48 *> \endverbatim
49 *
50 * Arguments:
51 * ==========
52 *
53 *> \param[in] JOB
54 *> \verbatim
55 *> JOB is CHARACTER*1
56 *> Specifies whether condition numbers are required for the
57 *> cluster of eigenvalues (S) or the invariant subspace (SEP):
58 *> = 'N': none;
59 *> = 'E': for eigenvalues only (S);
60 *> = 'V': for invariant subspace only (SEP);
61 *> = 'B': for both eigenvalues and invariant subspace (S and
62 *> SEP).
63 *> \endverbatim
64 *>
65 *> \param[in] COMPQ
66 *> \verbatim
67 *> COMPQ is CHARACTER*1
68 *> = 'V': update the matrix Q of Schur vectors;
69 *> = 'N': do not update Q.
70 *> \endverbatim
71 *>
72 *> \param[in] SELECT
73 *> \verbatim
74 *> SELECT is LOGICAL array, dimension (N)
75 *> SELECT specifies the eigenvalues in the selected cluster. To
76 *> select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
77 *> \endverbatim
78 *>
79 *> \param[in] N
80 *> \verbatim
81 *> N is INTEGER
82 *> The order of the matrix T. N >= 0.
83 *> \endverbatim
84 *>
85 *> \param[in,out] T
86 *> \verbatim
87 *> T is COMPLEX*16 array, dimension (LDT,N)
88 *> On entry, the upper triangular matrix T.
89 *> On exit, T is overwritten by the reordered matrix T, with the
90 *> selected eigenvalues as the leading diagonal elements.
91 *> \endverbatim
92 *>
93 *> \param[in] LDT
94 *> \verbatim
95 *> LDT is INTEGER
96 *> The leading dimension of the array T. LDT >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[in,out] Q
100 *> \verbatim
101 *> Q is COMPLEX*16 array, dimension (LDQ,N)
102 *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
103 *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
104 *> unitary transformation matrix which reorders T; the leading M
105 *> columns of Q form an orthonormal basis for the specified
106 *> invariant subspace.
107 *> If COMPQ = 'N', Q is not referenced.
108 *> \endverbatim
109 *>
110 *> \param[in] LDQ
111 *> \verbatim
112 *> LDQ is INTEGER
113 *> The leading dimension of the array Q.
114 *> LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
115 *> \endverbatim
116 *>
117 *> \param[out] W
118 *> \verbatim
119 *> W is COMPLEX*16 array, dimension (N)
120 *> The reordered eigenvalues of T, in the same order as they
121 *> appear on the diagonal of T.
122 *> \endverbatim
123 *>
124 *> \param[out] M
125 *> \verbatim
126 *> M is INTEGER
127 *> The dimension of the specified invariant subspace.
128 *> 0 <= M <= N.
129 *> \endverbatim
130 *>
131 *> \param[out] S
132 *> \verbatim
133 *> S is DOUBLE PRECISION
134 *> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
135 *> condition number for the selected cluster of eigenvalues.
136 *> S cannot underestimate the true reciprocal condition number
137 *> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
138 *> If JOB = 'N' or 'V', S is not referenced.
139 *> \endverbatim
140 *>
141 *> \param[out] SEP
142 *> \verbatim
143 *> SEP is DOUBLE PRECISION
144 *> If JOB = 'V' or 'B', SEP is the estimated reciprocal
145 *> condition number of the specified invariant subspace. If
146 *> M = 0 or N, SEP = norm(T).
147 *> If JOB = 'N' or 'E', SEP is not referenced.
148 *> \endverbatim
149 *>
150 *> \param[out] WORK
151 *> \verbatim
152 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
153 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
154 *> \endverbatim
155 *>
156 *> \param[in] LWORK
157 *> \verbatim
158 *> LWORK is INTEGER
159 *> The dimension of the array WORK.
160 *> If JOB = 'N', LWORK >= 1;
161 *> if JOB = 'E', LWORK = max(1,M*(N-M));
162 *> if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
163 *>
164 *> If LWORK = -1, then a workspace query is assumed; the routine
165 *> only calculates the optimal size of the WORK array, returns
166 *> this value as the first entry of the WORK array, and no error
167 *> message related to LWORK is issued by XERBLA.
168 *> \endverbatim
169 *>
170 *> \param[out] INFO
171 *> \verbatim
172 *> INFO is INTEGER
173 *> = 0: successful exit
174 *> < 0: if INFO = -i, the i-th argument had an illegal value
175 *> \endverbatim
176 *
177 * Authors:
178 * ========
179 *
180 *> \author Univ. of Tennessee
181 *> \author Univ. of California Berkeley
182 *> \author Univ. of Colorado Denver
183 *> \author NAG Ltd.
184 *
185 *> \ingroup complex16OTHERcomputational
186 *
187 *> \par Further Details:
188 * =====================
189 *>
190 *> \verbatim
191 *>
192 *> ZTRSEN first collects the selected eigenvalues by computing a unitary
193 *> transformation Z to move them to the top left corner of T. In other
194 *> words, the selected eigenvalues are the eigenvalues of T11 in:
195 *>
196 *> Z**H * T * Z = ( T11 T12 ) n1
197 *> ( 0 T22 ) n2
198 *> n1 n2
199 *>
200 *> where N = n1+n2. The first
201 *> n1 columns of Z span the specified invariant subspace of T.
202 *>
203 *> If T has been obtained from the Schur factorization of a matrix
204 *> A = Q*T*Q**H, then the reordered Schur factorization of A is given by
205 *> A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
206 *> corresponding invariant subspace of A.
207 *>
208 *> The reciprocal condition number of the average of the eigenvalues of
209 *> T11 may be returned in S. S lies between 0 (very badly conditioned)
210 *> and 1 (very well conditioned). It is computed as follows. First we
211 *> compute R so that
212 *>
213 *> P = ( I R ) n1
214 *> ( 0 0 ) n2
215 *> n1 n2
216 *>
217 *> is the projector on the invariant subspace associated with T11.
218 *> R is the solution of the Sylvester equation:
219 *>
220 *> T11*R - R*T22 = T12.
221 *>
222 *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
223 *> the two-norm of M. Then S is computed as the lower bound
224 *>
225 *> (1 + F-norm(R)**2)**(-1/2)
226 *>
227 *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
228 *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
229 *> sqrt(N).
230 *>
231 *> An approximate error bound for the computed average of the
232 *> eigenvalues of T11 is
233 *>
234 *> EPS * norm(T) / S
235 *>
236 *> where EPS is the machine precision.
237 *>
238 *> The reciprocal condition number of the right invariant subspace
239 *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
240 *> SEP is defined as the separation of T11 and T22:
241 *>
242 *> sep( T11, T22 ) = sigma-min( C )
243 *>
244 *> where sigma-min(C) is the smallest singular value of the
245 *> n1*n2-by-n1*n2 matrix
246 *>
247 *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
248 *>
249 *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
250 *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
251 *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
252 *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
253 *>
254 *> When SEP is small, small changes in T can cause large changes in
255 *> the invariant subspace. An approximate bound on the maximum angular
256 *> error in the computed right invariant subspace is
257 *>
258 *> EPS * norm(T) / SEP
259 *> \endverbatim
260 *>
261 * =====================================================================
262  SUBROUTINE ztrsen( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
263  $ SEP, WORK, LWORK, INFO )
264 *
265 * -- LAPACK computational routine --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 *
269 * .. Scalar Arguments ..
270  CHARACTER COMPQ, JOB
271  INTEGER INFO, LDQ, LDT, LWORK, M, N
272  DOUBLE PRECISION S, SEP
273 * ..
274 * .. Array Arguments ..
275  LOGICAL SELECT( * )
276  COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  DOUBLE PRECISION ZERO, ONE
283  parameter( zero = 0.0d+0, one = 1.0d+0 )
284 * ..
285 * .. Local Scalars ..
286  LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
287  INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
288  DOUBLE PRECISION EST, RNORM, SCALE
289 * ..
290 * .. Local Arrays ..
291  INTEGER ISAVE( 3 )
292  DOUBLE PRECISION RWORK( 1 )
293 * ..
294 * .. External Functions ..
295  LOGICAL LSAME
296  DOUBLE PRECISION ZLANGE
297  EXTERNAL lsame, zlange
298 * ..
299 * .. External Subroutines ..
300  EXTERNAL xerbla, zlacn2, zlacpy, ztrexc, ztrsyl
301 * ..
302 * .. Intrinsic Functions ..
303  INTRINSIC max, sqrt
304 * ..
305 * .. Executable Statements ..
306 *
307 * Decode and test the input parameters.
308 *
309  wantbh = lsame( job, 'B' )
310  wants = lsame( job, 'E' ) .OR. wantbh
311  wantsp = lsame( job, 'V' ) .OR. wantbh
312  wantq = lsame( compq, 'V' )
313 *
314 * Set M to the number of selected eigenvalues.
315 *
316  m = 0
317  DO 10 k = 1, n
318  IF( SELECT( k ) )
319  $ m = m + 1
320  10 CONTINUE
321 *
322  n1 = m
323  n2 = n - m
324  nn = n1*n2
325 *
326  info = 0
327  lquery = ( lwork.EQ.-1 )
328 *
329  IF( wantsp ) THEN
330  lwmin = max( 1, 2*nn )
331  ELSE IF( lsame( job, 'N' ) ) THEN
332  lwmin = 1
333  ELSE IF( lsame( job, 'E' ) ) THEN
334  lwmin = max( 1, nn )
335  END IF
336 *
337  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
338  $ THEN
339  info = -1
340  ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
341  info = -2
342  ELSE IF( n.LT.0 ) THEN
343  info = -4
344  ELSE IF( ldt.LT.max( 1, n ) ) THEN
345  info = -6
346  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
347  info = -8
348  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
349  info = -14
350  END IF
351 *
352  IF( info.EQ.0 ) THEN
353  work( 1 ) = lwmin
354  END IF
355 *
356  IF( info.NE.0 ) THEN
357  CALL xerbla( 'ZTRSEN', -info )
358  RETURN
359  ELSE IF( lquery ) THEN
360  RETURN
361  END IF
362 *
363 * Quick return if possible
364 *
365  IF( m.EQ.n .OR. m.EQ.0 ) THEN
366  IF( wants )
367  $ s = one
368  IF( wantsp )
369  $ sep = zlange( '1', n, n, t, ldt, rwork )
370  GO TO 40
371  END IF
372 *
373 * Collect the selected eigenvalues at the top left corner of T.
374 *
375  ks = 0
376  DO 20 k = 1, n
377  IF( SELECT( k ) ) THEN
378  ks = ks + 1
379 *
380 * Swap the K-th eigenvalue to position KS.
381 *
382  IF( k.NE.ks )
383  $ CALL ztrexc( compq, n, t, ldt, q, ldq, k, ks, ierr )
384  END IF
385  20 CONTINUE
386 *
387  IF( wants ) THEN
388 *
389 * Solve the Sylvester equation for R:
390 *
391 * T11*R - R*T22 = scale*T12
392 *
393  CALL zlacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
394  CALL ztrsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
395  $ ldt, work, n1, scale, ierr )
396 *
397 * Estimate the reciprocal of the condition number of the cluster
398 * of eigenvalues.
399 *
400  rnorm = zlange( 'F', n1, n2, work, n1, rwork )
401  IF( rnorm.EQ.zero ) THEN
402  s = one
403  ELSE
404  s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
405  $ sqrt( rnorm ) )
406  END IF
407  END IF
408 *
409  IF( wantsp ) THEN
410 *
411 * Estimate sep(T11,T22).
412 *
413  est = zero
414  kase = 0
415  30 CONTINUE
416  CALL zlacn2( nn, work( nn+1 ), work, est, kase, isave )
417  IF( kase.NE.0 ) THEN
418  IF( kase.EQ.1 ) THEN
419 *
420 * Solve T11*R - R*T22 = scale*X.
421 *
422  CALL ztrsyl( 'N', 'N', -1, n1, n2, t, ldt,
423  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
424  $ ierr )
425  ELSE
426 *
427 * Solve T11**H*R - R*T22**H = scale*X.
428 *
429  CALL ztrsyl( 'C', 'C', -1, n1, n2, t, ldt,
430  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
431  $ ierr )
432  END IF
433  GO TO 30
434  END IF
435 *
436  sep = scale / est
437  END IF
438 *
439  40 CONTINUE
440 *
441 * Copy reordered eigenvalues to W.
442 *
443  DO 50 k = 1, n
444  w( k ) = t( k, k )
445  50 CONTINUE
446 *
447  work( 1 ) = lwmin
448 *
449  RETURN
450 *
451 * End of ZTRSEN
452 *
453  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:133
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
Definition: ztrexc.f:126
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:264
subroutine ztrsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
ZTRSYL
Definition: ztrsyl.f:157