LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dtrsen.f
Go to the documentation of this file.
1 *> \brief \b DTRSEN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DTRSEN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsen.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsen.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsen.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
22 * M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER COMPQ, JOB
26 * INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
27 * DOUBLE PRECISION S, SEP
28 * ..
29 * .. Array Arguments ..
30 * LOGICAL SELECT( * )
31 * INTEGER IWORK( * )
32 * DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
33 * $ WR( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> DTRSEN reorders the real Schur factorization of a real matrix
43 *> A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in
44 *> the leading diagonal blocks of the upper quasi-triangular matrix T,
45 *> and the leading columns of Q form an orthonormal basis of the
46 *> corresponding right invariant subspace.
47 *>
48 *> Optionally the routine computes the reciprocal condition numbers of
49 *> the cluster of eigenvalues and/or the invariant subspace.
50 *>
51 *> T must be in Schur canonical form (as returned by DHSEQR), that is,
52 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
53 *> 2-by-2 diagonal block has its diagonal elements equal and its
54 *> off-diagonal elements of opposite sign.
55 *> \endverbatim
56 *
57 * Arguments:
58 * ==========
59 *
60 *> \param[in] JOB
61 *> \verbatim
62 *> JOB is CHARACTER*1
63 *> Specifies whether condition numbers are required for the
64 *> cluster of eigenvalues (S) or the invariant subspace (SEP):
65 *> = 'N': none;
66 *> = 'E': for eigenvalues only (S);
67 *> = 'V': for invariant subspace only (SEP);
68 *> = 'B': for both eigenvalues and invariant subspace (S and
69 *> SEP).
70 *> \endverbatim
71 *>
72 *> \param[in] COMPQ
73 *> \verbatim
74 *> COMPQ is CHARACTER*1
75 *> = 'V': update the matrix Q of Schur vectors;
76 *> = 'N': do not update Q.
77 *> \endverbatim
78 *>
79 *> \param[in] SELECT
80 *> \verbatim
81 *> SELECT is LOGICAL array, dimension (N)
82 *> SELECT specifies the eigenvalues in the selected cluster. To
83 *> select a real eigenvalue w(j), SELECT(j) must be set to
84 *> .TRUE.. To select a complex conjugate pair of eigenvalues
85 *> w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
86 *> either SELECT(j) or SELECT(j+1) or both must be set to
87 *> .TRUE.; a complex conjugate pair of eigenvalues must be
88 *> either both included in the cluster or both excluded.
89 *> \endverbatim
90 *>
91 *> \param[in] N
92 *> \verbatim
93 *> N is INTEGER
94 *> The order of the matrix T. N >= 0.
95 *> \endverbatim
96 *>
97 *> \param[in,out] T
98 *> \verbatim
99 *> T is DOUBLE PRECISION array, dimension (LDT,N)
100 *> On entry, the upper quasi-triangular matrix T, in Schur
101 *> canonical form.
102 *> On exit, T is overwritten by the reordered matrix T, again in
103 *> Schur canonical form, with the selected eigenvalues in the
104 *> leading diagonal blocks.
105 *> \endverbatim
106 *>
107 *> \param[in] LDT
108 *> \verbatim
109 *> LDT is INTEGER
110 *> The leading dimension of the array T. LDT >= max(1,N).
111 *> \endverbatim
112 *>
113 *> \param[in,out] Q
114 *> \verbatim
115 *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
116 *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
117 *> On exit, if COMPQ = 'V', Q has been postmultiplied by the
118 *> orthogonal transformation matrix which reorders T; the
119 *> leading M columns of Q form an orthonormal basis for the
120 *> specified invariant subspace.
121 *> If COMPQ = 'N', Q is not referenced.
122 *> \endverbatim
123 *>
124 *> \param[in] LDQ
125 *> \verbatim
126 *> LDQ is INTEGER
127 *> The leading dimension of the array Q.
128 *> LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
129 *> \endverbatim
130 *>
131 *> \param[out] WR
132 *> \verbatim
133 *> WR is DOUBLE PRECISION array, dimension (N)
134 *> \endverbatim
135 *> \param[out] WI
136 *> \verbatim
137 *> WI is DOUBLE PRECISION array, dimension (N)
138 *>
139 *> The real and imaginary parts, respectively, of the reordered
140 *> eigenvalues of T. The eigenvalues are stored in the same
141 *> order as on the diagonal of T, with WR(i) = T(i,i) and, if
142 *> T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
143 *> WI(i+1) = -WI(i). Note that if a complex eigenvalue is
144 *> sufficiently ill-conditioned, then its value may differ
145 *> significantly from its value before reordering.
146 *> \endverbatim
147 *>
148 *> \param[out] M
149 *> \verbatim
150 *> M is INTEGER
151 *> The dimension of the specified invariant subspace.
152 *> 0 < = M <= N.
153 *> \endverbatim
154 *>
155 *> \param[out] S
156 *> \verbatim
157 *> S is DOUBLE PRECISION
158 *> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
159 *> condition number for the selected cluster of eigenvalues.
160 *> S cannot underestimate the true reciprocal condition number
161 *> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
162 *> If JOB = 'N' or 'V', S is not referenced.
163 *> \endverbatim
164 *>
165 *> \param[out] SEP
166 *> \verbatim
167 *> SEP is DOUBLE PRECISION
168 *> If JOB = 'V' or 'B', SEP is the estimated reciprocal
169 *> condition number of the specified invariant subspace. If
170 *> M = 0 or N, SEP = norm(T).
171 *> If JOB = 'N' or 'E', SEP is not referenced.
172 *> \endverbatim
173 *>
174 *> \param[out] WORK
175 *> \verbatim
176 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
177 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
178 *> \endverbatim
179 *>
180 *> \param[in] LWORK
181 *> \verbatim
182 *> LWORK is INTEGER
183 *> The dimension of the array WORK.
184 *> If JOB = 'N', LWORK >= max(1,N);
185 *> if JOB = 'E', LWORK >= max(1,M*(N-M));
186 *> if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
187 *>
188 *> If LWORK = -1, then a workspace query is assumed; the routine
189 *> only calculates the optimal size of the WORK array, returns
190 *> this value as the first entry of the WORK array, and no error
191 *> message related to LWORK is issued by XERBLA.
192 *> \endverbatim
193 *>
194 *> \param[out] IWORK
195 *> \verbatim
196 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
197 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
198 *> \endverbatim
199 *>
200 *> \param[in] LIWORK
201 *> \verbatim
202 *> LIWORK is INTEGER
203 *> The dimension of the array IWORK.
204 *> If JOB = 'N' or 'E', LIWORK >= 1;
205 *> if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
206 *>
207 *> If LIWORK = -1, then a workspace query is assumed; the
208 *> routine only calculates the optimal size of the IWORK array,
209 *> returns this value as the first entry of the IWORK array, and
210 *> no error message related to LIWORK is issued by XERBLA.
211 *> \endverbatim
212 *>
213 *> \param[out] INFO
214 *> \verbatim
215 *> INFO is INTEGER
216 *> = 0: successful exit
217 *> < 0: if INFO = -i, the i-th argument had an illegal value
218 *> = 1: reordering of T failed because some eigenvalues are too
219 *> close to separate (the problem is very ill-conditioned);
220 *> T may have been partially reordered, and WR and WI
221 *> contain the eigenvalues in the same order as in T; S and
222 *> SEP (if requested) are set to zero.
223 *> \endverbatim
224 *
225 * Authors:
226 * ========
227 *
228 *> \author Univ. of Tennessee
229 *> \author Univ. of California Berkeley
230 *> \author Univ. of Colorado Denver
231 *> \author NAG Ltd.
232 *
233 *> \date April 2012
234 *
235 *> \ingroup doubleOTHERcomputational
236 *
237 *> \par Further Details:
238 * =====================
239 *>
240 *> \verbatim
241 *>
242 *> DTRSEN first collects the selected eigenvalues by computing an
243 *> orthogonal transformation Z to move them to the top left corner of T.
244 *> In other words, the selected eigenvalues are the eigenvalues of T11
245 *> in:
246 *>
247 *> Z**T * T * Z = ( T11 T12 ) n1
248 *> ( 0 T22 ) n2
249 *> n1 n2
250 *>
251 *> where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
252 *> of Z span the specified invariant subspace of T.
253 *>
254 *> If T has been obtained from the real Schur factorization of a matrix
255 *> A = Q*T*Q**T, then the reordered real Schur factorization of A is given
256 *> by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span
257 *> the corresponding invariant subspace of A.
258 *>
259 *> The reciprocal condition number of the average of the eigenvalues of
260 *> T11 may be returned in S. S lies between 0 (very badly conditioned)
261 *> and 1 (very well conditioned). It is computed as follows. First we
262 *> compute R so that
263 *>
264 *> P = ( I R ) n1
265 *> ( 0 0 ) n2
266 *> n1 n2
267 *>
268 *> is the projector on the invariant subspace associated with T11.
269 *> R is the solution of the Sylvester equation:
270 *>
271 *> T11*R - R*T22 = T12.
272 *>
273 *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
274 *> the two-norm of M. Then S is computed as the lower bound
275 *>
276 *> (1 + F-norm(R)**2)**(-1/2)
277 *>
278 *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
279 *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
280 *> sqrt(N).
281 *>
282 *> An approximate error bound for the computed average of the
283 *> eigenvalues of T11 is
284 *>
285 *> EPS * norm(T) / S
286 *>
287 *> where EPS is the machine precision.
288 *>
289 *> The reciprocal condition number of the right invariant subspace
290 *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
291 *> SEP is defined as the separation of T11 and T22:
292 *>
293 *> sep( T11, T22 ) = sigma-min( C )
294 *>
295 *> where sigma-min(C) is the smallest singular value of the
296 *> n1*n2-by-n1*n2 matrix
297 *>
298 *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
299 *>
300 *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
301 *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
302 *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
303 *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
304 *>
305 *> When SEP is small, small changes in T can cause large changes in
306 *> the invariant subspace. An approximate bound on the maximum angular
307 *> error in the computed right invariant subspace is
308 *>
309 *> EPS * norm(T) / SEP
310 *> \endverbatim
311 *>
312 * =====================================================================
313  SUBROUTINE dtrsen( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
314  $ m, s, sep, work, lwork, iwork, liwork, info )
315 *
316 * -- LAPACK computational routine (version 3.4.1) --
317 * -- LAPACK is a software package provided by Univ. of Tennessee, --
318 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
319 * April 2012
320 *
321 * .. Scalar Arguments ..
322  CHARACTER compq, job
323  INTEGER info, ldq, ldt, liwork, lwork, m, n
324  DOUBLE PRECISION s, sep
325 * ..
326 * .. Array Arguments ..
327  LOGICAL select( * )
328  INTEGER iwork( * )
329  DOUBLE PRECISION q( ldq, * ), t( ldt, * ), wi( * ), work( * ),
330  $ wr( * )
331 * ..
332 *
333 * =====================================================================
334 *
335 * .. Parameters ..
336  DOUBLE PRECISION zero, one
337  parameter( zero = 0.0d+0, one = 1.0d+0 )
338 * ..
339 * .. Local Scalars ..
340  LOGICAL lquery, pair, swap, wantbh, wantq, wants,
341  $ wantsp
342  INTEGER ierr, k, kase, kk, ks, liwmin, lwmin, n1, n2,
343  $ nn
344  DOUBLE PRECISION est, rnorm, scale
345 * ..
346 * .. Local Arrays ..
347  INTEGER isave( 3 )
348 * ..
349 * .. External Functions ..
350  LOGICAL lsame
351  DOUBLE PRECISION dlange
352  EXTERNAL lsame, dlange
353 * ..
354 * .. External Subroutines ..
355  EXTERNAL dlacn2, dlacpy, dtrexc, dtrsyl, xerbla
356 * ..
357 * .. Intrinsic Functions ..
358  INTRINSIC abs, max, sqrt
359 * ..
360 * .. Executable Statements ..
361 *
362 * Decode and test the input parameters
363 *
364  wantbh = lsame( job, 'B' )
365  wants = lsame( job, 'E' ) .OR. wantbh
366  wantsp = lsame( job, 'V' ) .OR. wantbh
367  wantq = lsame( compq, 'V' )
368 *
369  info = 0
370  lquery = ( lwork.EQ.-1 )
371  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
372  $ THEN
373  info = -1
374  ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
375  info = -2
376  ELSE IF( n.LT.0 ) THEN
377  info = -4
378  ELSE IF( ldt.LT.max( 1, n ) ) THEN
379  info = -6
380  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
381  info = -8
382  ELSE
383 *
384 * Set M to the dimension of the specified invariant subspace,
385 * and test LWORK and LIWORK.
386 *
387  m = 0
388  pair = .false.
389  DO 10 k = 1, n
390  IF( pair ) THEN
391  pair = .false.
392  ELSE
393  IF( k.LT.n ) THEN
394  IF( t( k+1, k ).EQ.zero ) THEN
395  IF( SELECT( k ) )
396  $ m = m + 1
397  ELSE
398  pair = .true.
399  IF( SELECT( k ) .OR. SELECT( k+1 ) )
400  $ m = m + 2
401  END IF
402  ELSE
403  IF( SELECT( n ) )
404  $ m = m + 1
405  END IF
406  END IF
407  10 continue
408 *
409  n1 = m
410  n2 = n - m
411  nn = n1*n2
412 *
413  IF( wantsp ) THEN
414  lwmin = max( 1, 2*nn )
415  liwmin = max( 1, nn )
416  ELSE IF( lsame( job, 'N' ) ) THEN
417  lwmin = max( 1, n )
418  liwmin = 1
419  ELSE IF( lsame( job, 'E' ) ) THEN
420  lwmin = max( 1, nn )
421  liwmin = 1
422  END IF
423 *
424  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
425  info = -15
426  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
427  info = -17
428  END IF
429  END IF
430 *
431  IF( info.EQ.0 ) THEN
432  work( 1 ) = lwmin
433  iwork( 1 ) = liwmin
434  END IF
435 *
436  IF( info.NE.0 ) THEN
437  CALL xerbla( 'DTRSEN', -info )
438  return
439  ELSE IF( lquery ) THEN
440  return
441  END IF
442 *
443 * Quick return if possible.
444 *
445  IF( m.EQ.n .OR. m.EQ.0 ) THEN
446  IF( wants )
447  $ s = one
448  IF( wantsp )
449  $ sep = dlange( '1', n, n, t, ldt, work )
450  go to 40
451  END IF
452 *
453 * Collect the selected blocks at the top-left corner of T.
454 *
455  ks = 0
456  pair = .false.
457  DO 20 k = 1, n
458  IF( pair ) THEN
459  pair = .false.
460  ELSE
461  swap = SELECT( k )
462  IF( k.LT.n ) THEN
463  IF( t( k+1, k ).NE.zero ) THEN
464  pair = .true.
465  swap = swap .OR. SELECT( k+1 )
466  END IF
467  END IF
468  IF( swap ) THEN
469  ks = ks + 1
470 *
471 * Swap the K-th block to position KS.
472 *
473  ierr = 0
474  kk = k
475  IF( k.NE.ks )
476  $ CALL dtrexc( compq, n, t, ldt, q, ldq, kk, ks, work,
477  $ ierr )
478  IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
479 *
480 * Blocks too close to swap: exit.
481 *
482  info = 1
483  IF( wants )
484  $ s = zero
485  IF( wantsp )
486  $ sep = zero
487  go to 40
488  END IF
489  IF( pair )
490  $ ks = ks + 1
491  END IF
492  END IF
493  20 continue
494 *
495  IF( wants ) THEN
496 *
497 * Solve Sylvester equation for R:
498 *
499 * T11*R - R*T22 = scale*T12
500 *
501  CALL dlacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
502  CALL dtrsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
503  $ ldt, work, n1, scale, ierr )
504 *
505 * Estimate the reciprocal of the condition number of the cluster
506 * of eigenvalues.
507 *
508  rnorm = dlange( 'F', n1, n2, work, n1, work )
509  IF( rnorm.EQ.zero ) THEN
510  s = one
511  ELSE
512  s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
513  $ sqrt( rnorm ) )
514  END IF
515  END IF
516 *
517  IF( wantsp ) THEN
518 *
519 * Estimate sep(T11,T22).
520 *
521  est = zero
522  kase = 0
523  30 continue
524  CALL dlacn2( nn, work( nn+1 ), work, iwork, est, kase, isave )
525  IF( kase.NE.0 ) THEN
526  IF( kase.EQ.1 ) THEN
527 *
528 * Solve T11*R - R*T22 = scale*X.
529 *
530  CALL dtrsyl( 'N', 'N', -1, n1, n2, t, ldt,
531  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
532  $ ierr )
533  ELSE
534 *
535 * Solve T11**T*R - R*T22**T = scale*X.
536 *
537  CALL dtrsyl( 'T', 'T', -1, n1, n2, t, ldt,
538  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
539  $ ierr )
540  END IF
541  go to 30
542  END IF
543 *
544  sep = scale / est
545  END IF
546 *
547  40 continue
548 *
549 * Store the output eigenvalues in WR and WI.
550 *
551  DO 50 k = 1, n
552  wr( k ) = t( k, k )
553  wi( k ) = zero
554  50 continue
555  DO 60 k = 1, n - 1
556  IF( t( k+1, k ).NE.zero ) THEN
557  wi( k ) = sqrt( abs( t( k, k+1 ) ) )*
558  $ sqrt( abs( t( k+1, k ) ) )
559  wi( k+1 ) = -wi( k )
560  END IF
561  60 continue
562 *
563  work( 1 ) = lwmin
564  iwork( 1 ) = liwmin
565 *
566  return
567 *
568 * End of DTRSEN
569 *
570  END