LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
186 *
187 *> \ingroup complex16OTHERcomputational
188 *
189 *> \par Further Details:
190 * =====================
191 *>
192 *> \verbatim
193 *>
194 *> ZTRSEN first collects the selected eigenvalues by computing a unitary
195 *> transformation Z to move them to the top left corner of T. In other
196 *> words, the selected eigenvalues are the eigenvalues of T11 in:
197 *>
198 *> Z**H * T * Z = ( T11 T12 ) n1
199 *> ( 0 T22 ) n2
200 *> n1 n2
201 *>
202 *> where N = n1+n2. The first
203 *> n1 columns of Z span the specified invariant subspace of T.
204 *>
205 *> If T has been obtained from the Schur factorization of a matrix
206 *> A = Q*T*Q**H, then the reordered Schur factorization of A is given by
207 *> A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
208 *> corresponding invariant subspace of A.
209 *>
210 *> The reciprocal condition number of the average of the eigenvalues of
211 *> T11 may be returned in S. S lies between 0 (very badly conditioned)
212 *> and 1 (very well conditioned). It is computed as follows. First we
213 *> compute R so that
214 *>
215 *> P = ( I R ) n1
216 *> ( 0 0 ) n2
217 *> n1 n2
218 *>
219 *> is the projector on the invariant subspace associated with T11.
220 *> R is the solution of the Sylvester equation:
221 *>
222 *> T11*R - R*T22 = T12.
223 *>
224 *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
225 *> the two-norm of M. Then S is computed as the lower bound
226 *>
227 *> (1 + F-norm(R)**2)**(-1/2)
228 *>
229 *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
230 *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
231 *> sqrt(N).
232 *>
233 *> An approximate error bound for the computed average of the
234 *> eigenvalues of T11 is
235 *>
236 *> EPS * norm(T) / S
237 *>
238 *> where EPS is the machine precision.
239 *>
240 *> The reciprocal condition number of the right invariant subspace
241 *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
242 *> SEP is defined as the separation of T11 and T22:
243 *>
244 *> sep( T11, T22 ) = sigma-min( C )
245 *>
246 *> where sigma-min(C) is the smallest singular value of the
247 *> n1*n2-by-n1*n2 matrix
248 *>
249 *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
250 *>
251 *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
252 *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
253 *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
254 *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
255 *>
256 *> When SEP is small, small changes in T can cause large changes in
257 *> the invariant subspace. An approximate bound on the maximum angular
258 *> error in the computed right invariant subspace is
259 *>
260 *> EPS * norm(T) / SEP
261 *> \endverbatim
262 *>
263 * =====================================================================
264  SUBROUTINE ztrsen( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
265  $ sep, work, lwork, info )
266 *
267 * -- LAPACK computational routine (version 3.4.0) --
268 * -- LAPACK is a software package provided by Univ. of Tennessee, --
269 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
270 * November 2011
271 *
272 * .. Scalar Arguments ..
273  CHARACTER compq, job
274  INTEGER info, ldq, ldt, lwork, m, n
275  DOUBLE PRECISION s, sep
276 * ..
277 * .. Array Arguments ..
278  LOGICAL select( * )
279  COMPLEX*16 q( ldq, * ), t( ldt, * ), w( * ), work( * )
280 * ..
281 *
282 * =====================================================================
283 *
284 * .. Parameters ..
285  DOUBLE PRECISION zero, one
286  parameter( zero = 0.0d+0, one = 1.0d+0 )
287 * ..
288 * .. Local Scalars ..
289  LOGICAL lquery, wantbh, wantq, wants, wantsp
290  INTEGER ierr, k, kase, ks, lwmin, n1, n2, nn
291  DOUBLE PRECISION est, rnorm, scale
292 * ..
293 * .. Local Arrays ..
294  INTEGER isave( 3 )
295  DOUBLE PRECISION rwork( 1 )
296 * ..
297 * .. External Functions ..
298  LOGICAL lsame
299  DOUBLE PRECISION zlange
300  EXTERNAL lsame, zlange
301 * ..
302 * .. External Subroutines ..
303  EXTERNAL xerbla, zlacn2, zlacpy, ztrexc, ztrsyl
304 * ..
305 * .. Intrinsic Functions ..
306  INTRINSIC max, sqrt
307 * ..
308 * .. Executable Statements ..
309 *
310 * Decode and test the input parameters.
311 *
312  wantbh = lsame( job, 'B' )
313  wants = lsame( job, 'E' ) .OR. wantbh
314  wantsp = lsame( job, 'V' ) .OR. wantbh
315  wantq = lsame( compq, 'V' )
316 *
317 * Set M to the number of selected eigenvalues.
318 *
319  m = 0
320  DO 10 k = 1, n
321  IF( SELECT( k ) )
322  $ m = m + 1
323  10 continue
324 *
325  n1 = m
326  n2 = n - m
327  nn = n1*n2
328 *
329  info = 0
330  lquery = ( lwork.EQ.-1 )
331 *
332  IF( wantsp ) THEN
333  lwmin = max( 1, 2*nn )
334  ELSE IF( lsame( job, 'N' ) ) THEN
335  lwmin = 1
336  ELSE IF( lsame( job, 'E' ) ) THEN
337  lwmin = max( 1, nn )
338  END IF
339 *
340  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
341  $ THEN
342  info = -1
343  ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
344  info = -2
345  ELSE IF( n.LT.0 ) THEN
346  info = -4
347  ELSE IF( ldt.LT.max( 1, n ) ) THEN
348  info = -6
349  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
350  info = -8
351  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
352  info = -14
353  END IF
354 *
355  IF( info.EQ.0 ) THEN
356  work( 1 ) = lwmin
357  END IF
358 *
359  IF( info.NE.0 ) THEN
360  CALL xerbla( 'ZTRSEN', -info )
361  return
362  ELSE IF( lquery ) THEN
363  return
364  END IF
365 *
366 * Quick return if possible
367 *
368  IF( m.EQ.n .OR. m.EQ.0 ) THEN
369  IF( wants )
370  $ s = one
371  IF( wantsp )
372  $ sep = zlange( '1', n, n, t, ldt, rwork )
373  go to 40
374  END IF
375 *
376 * Collect the selected eigenvalues at the top left corner of T.
377 *
378  ks = 0
379  DO 20 k = 1, n
380  IF( SELECT( k ) ) THEN
381  ks = ks + 1
382 *
383 * Swap the K-th eigenvalue to position KS.
384 *
385  IF( k.NE.ks )
386  $ CALL ztrexc( compq, n, t, ldt, q, ldq, k, ks, ierr )
387  END IF
388  20 continue
389 *
390  IF( wants ) THEN
391 *
392 * Solve the Sylvester equation for R:
393 *
394 * T11*R - R*T22 = scale*T12
395 *
396  CALL zlacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
397  CALL ztrsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
398  $ ldt, work, n1, scale, ierr )
399 *
400 * Estimate the reciprocal of the condition number of the cluster
401 * of eigenvalues.
402 *
403  rnorm = zlange( 'F', n1, n2, work, n1, rwork )
404  IF( rnorm.EQ.zero ) THEN
405  s = one
406  ELSE
407  s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
408  $ sqrt( rnorm ) )
409  END IF
410  END IF
411 *
412  IF( wantsp ) THEN
413 *
414 * Estimate sep(T11,T22).
415 *
416  est = zero
417  kase = 0
418  30 continue
419  CALL zlacn2( nn, work( nn+1 ), work, est, kase, isave )
420  IF( kase.NE.0 ) THEN
421  IF( kase.EQ.1 ) THEN
422 *
423 * Solve T11*R - R*T22 = scale*X.
424 *
425  CALL ztrsyl( 'N', 'N', -1, n1, n2, t, ldt,
426  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
427  $ ierr )
428  ELSE
429 *
430 * Solve T11**H*R - R*T22**H = scale*X.
431 *
432  CALL ztrsyl( 'C', 'C', -1, n1, n2, t, ldt,
433  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
434  $ ierr )
435  END IF
436  go to 30
437  END IF
438 *
439  sep = scale / est
440  END IF
441 *
442  40 continue
443 *
444 * Copy reordered eigenvalues to W.
445 *
446  DO 50 k = 1, n
447  w( k ) = t( k, k )
448  50 continue
449 *
450  work( 1 ) = lwmin
451 *
452  return
453 *
454 * End of ZTRSEN
455 *
456  END