LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup doubleOTHERcomputational
234 *
235 *> \par Further Details:
236 * =====================
237 *>
238 *> \verbatim
239 *>
240 *> DTRSEN first collects the selected eigenvalues by computing an
241 *> orthogonal transformation Z to move them to the top left corner of T.
242 *> In other words, the selected eigenvalues are the eigenvalues of T11
243 *> in:
244 *>
245 *> Z**T * T * Z = ( T11 T12 ) n1
246 *> ( 0 T22 ) n2
247 *> n1 n2
248 *>
249 *> where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
250 *> of Z span the specified invariant subspace of T.
251 *>
252 *> If T has been obtained from the real Schur factorization of a matrix
253 *> A = Q*T*Q**T, then the reordered real Schur factorization of A is given
254 *> by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span
255 *> the corresponding invariant subspace of A.
256 *>
257 *> The reciprocal condition number of the average of the eigenvalues of
258 *> T11 may be returned in S. S lies between 0 (very badly conditioned)
259 *> and 1 (very well conditioned). It is computed as follows. First we
260 *> compute R so that
261 *>
262 *> P = ( I R ) n1
263 *> ( 0 0 ) n2
264 *> n1 n2
265 *>
266 *> is the projector on the invariant subspace associated with T11.
267 *> R is the solution of the Sylvester equation:
268 *>
269 *> T11*R - R*T22 = T12.
270 *>
271 *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
272 *> the two-norm of M. Then S is computed as the lower bound
273 *>
274 *> (1 + F-norm(R)**2)**(-1/2)
275 *>
276 *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
277 *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
278 *> sqrt(N).
279 *>
280 *> An approximate error bound for the computed average of the
281 *> eigenvalues of T11 is
282 *>
283 *> EPS * norm(T) / S
284 *>
285 *> where EPS is the machine precision.
286 *>
287 *> The reciprocal condition number of the right invariant subspace
288 *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
289 *> SEP is defined as the separation of T11 and T22:
290 *>
291 *> sep( T11, T22 ) = sigma-min( C )
292 *>
293 *> where sigma-min(C) is the smallest singular value of the
294 *> n1*n2-by-n1*n2 matrix
295 *>
296 *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
297 *>
298 *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
299 *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
300 *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
301 *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
302 *>
303 *> When SEP is small, small changes in T can cause large changes in
304 *> the invariant subspace. An approximate bound on the maximum angular
305 *> error in the computed right invariant subspace is
306 *>
307 *> EPS * norm(T) / SEP
308 *> \endverbatim
309 *>
310 * =====================================================================
311  SUBROUTINE dtrsen( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
312  $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
313 *
314 * -- LAPACK computational routine --
315 * -- LAPACK is a software package provided by Univ. of Tennessee, --
316 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
317 *
318 * .. Scalar Arguments ..
319  CHARACTER COMPQ, JOB
320  INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
321  DOUBLE PRECISION S, SEP
322 * ..
323 * .. Array Arguments ..
324  LOGICAL SELECT( * )
325  INTEGER IWORK( * )
326  DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
327  $ wr( * )
328 * ..
329 *
330 * =====================================================================
331 *
332 * .. Parameters ..
333  DOUBLE PRECISION ZERO, ONE
334  parameter( zero = 0.0d+0, one = 1.0d+0 )
335 * ..
336 * .. Local Scalars ..
337  LOGICAL LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS,
338  $ wantsp
339  INTEGER IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2,
340  $ nn
341  DOUBLE PRECISION EST, RNORM, SCALE
342 * ..
343 * .. Local Arrays ..
344  INTEGER ISAVE( 3 )
345 * ..
346 * .. External Functions ..
347  LOGICAL LSAME
348  DOUBLE PRECISION DLANGE
349  EXTERNAL lsame, dlange
350 * ..
351 * .. External Subroutines ..
352  EXTERNAL dlacn2, dlacpy, dtrexc, dtrsyl, xerbla
353 * ..
354 * .. Intrinsic Functions ..
355  INTRINSIC abs, max, sqrt
356 * ..
357 * .. Executable Statements ..
358 *
359 * Decode and test the input parameters
360 *
361  wantbh = lsame( job, 'B' )
362  wants = lsame( job, 'E' ) .OR. wantbh
363  wantsp = lsame( job, 'V' ) .OR. wantbh
364  wantq = lsame( compq, 'V' )
365 *
366  info = 0
367  lquery = ( lwork.EQ.-1 )
368  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
369  $ THEN
370  info = -1
371  ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
372  info = -2
373  ELSE IF( n.LT.0 ) THEN
374  info = -4
375  ELSE IF( ldt.LT.max( 1, n ) ) THEN
376  info = -6
377  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
378  info = -8
379  ELSE
380 *
381 * Set M to the dimension of the specified invariant subspace,
382 * and test LWORK and LIWORK.
383 *
384  m = 0
385  pair = .false.
386  DO 10 k = 1, n
387  IF( pair ) THEN
388  pair = .false.
389  ELSE
390  IF( k.LT.n ) THEN
391  IF( t( k+1, k ).EQ.zero ) THEN
392  IF( SELECT( k ) )
393  $ m = m + 1
394  ELSE
395  pair = .true.
396  IF( SELECT( k ) .OR. SELECT( k+1 ) )
397  $ m = m + 2
398  END IF
399  ELSE
400  IF( SELECT( n ) )
401  $ m = m + 1
402  END IF
403  END IF
404  10 CONTINUE
405 *
406  n1 = m
407  n2 = n - m
408  nn = n1*n2
409 *
410  IF( wantsp ) THEN
411  lwmin = max( 1, 2*nn )
412  liwmin = max( 1, nn )
413  ELSE IF( lsame( job, 'N' ) ) THEN
414  lwmin = max( 1, n )
415  liwmin = 1
416  ELSE IF( lsame( job, 'E' ) ) THEN
417  lwmin = max( 1, nn )
418  liwmin = 1
419  END IF
420 *
421  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
422  info = -15
423  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
424  info = -17
425  END IF
426  END IF
427 *
428  IF( info.EQ.0 ) THEN
429  work( 1 ) = lwmin
430  iwork( 1 ) = liwmin
431  END IF
432 *
433  IF( info.NE.0 ) THEN
434  CALL xerbla( 'DTRSEN', -info )
435  RETURN
436  ELSE IF( lquery ) THEN
437  RETURN
438  END IF
439 *
440 * Quick return if possible.
441 *
442  IF( m.EQ.n .OR. m.EQ.0 ) THEN
443  IF( wants )
444  $ s = one
445  IF( wantsp )
446  $ sep = dlange( '1', n, n, t, ldt, work )
447  GO TO 40
448  END IF
449 *
450 * Collect the selected blocks at the top-left corner of T.
451 *
452  ks = 0
453  pair = .false.
454  DO 20 k = 1, n
455  IF( pair ) THEN
456  pair = .false.
457  ELSE
458  swap = SELECT( k )
459  IF( k.LT.n ) THEN
460  IF( t( k+1, k ).NE.zero ) THEN
461  pair = .true.
462  swap = swap .OR. SELECT( k+1 )
463  END IF
464  END IF
465  IF( swap ) THEN
466  ks = ks + 1
467 *
468 * Swap the K-th block to position KS.
469 *
470  ierr = 0
471  kk = k
472  IF( k.NE.ks )
473  $ CALL dtrexc( compq, n, t, ldt, q, ldq, kk, ks, work,
474  $ ierr )
475  IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
476 *
477 * Blocks too close to swap: exit.
478 *
479  info = 1
480  IF( wants )
481  $ s = zero
482  IF( wantsp )
483  $ sep = zero
484  GO TO 40
485  END IF
486  IF( pair )
487  $ ks = ks + 1
488  END IF
489  END IF
490  20 CONTINUE
491 *
492  IF( wants ) THEN
493 *
494 * Solve Sylvester equation for R:
495 *
496 * T11*R - R*T22 = scale*T12
497 *
498  CALL dlacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
499  CALL dtrsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
500  $ ldt, work, n1, scale, ierr )
501 *
502 * Estimate the reciprocal of the condition number of the cluster
503 * of eigenvalues.
504 *
505  rnorm = dlange( 'F', n1, n2, work, n1, work )
506  IF( rnorm.EQ.zero ) THEN
507  s = one
508  ELSE
509  s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
510  $ sqrt( rnorm ) )
511  END IF
512  END IF
513 *
514  IF( wantsp ) THEN
515 *
516 * Estimate sep(T11,T22).
517 *
518  est = zero
519  kase = 0
520  30 CONTINUE
521  CALL dlacn2( nn, work( nn+1 ), work, iwork, est, kase, isave )
522  IF( kase.NE.0 ) THEN
523  IF( kase.EQ.1 ) THEN
524 *
525 * Solve T11*R - R*T22 = scale*X.
526 *
527  CALL dtrsyl( 'N', 'N', -1, n1, n2, t, ldt,
528  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
529  $ ierr )
530  ELSE
531 *
532 * Solve T11**T*R - R*T22**T = scale*X.
533 *
534  CALL dtrsyl( 'T', 'T', -1, n1, n2, t, ldt,
535  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
536  $ ierr )
537  END IF
538  GO TO 30
539  END IF
540 *
541  sep = scale / est
542  END IF
543 *
544  40 CONTINUE
545 *
546 * Store the output eigenvalues in WR and WI.
547 *
548  DO 50 k = 1, n
549  wr( k ) = t( k, k )
550  wi( k ) = zero
551  50 CONTINUE
552  DO 60 k = 1, n - 1
553  IF( t( k+1, k ).NE.zero ) THEN
554  wi( k ) = sqrt( abs( t( k, k+1 ) ) )*
555  $ sqrt( abs( t( k+1, k ) ) )
556  wi( k+1 ) = -wi( k )
557  END IF
558  60 CONTINUE
559 *
560  work( 1 ) = lwmin
561  iwork( 1 ) = liwmin
562 *
563  RETURN
564 *
565 * End of DTRSEN
566 *
567  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:136
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
Definition: dtrexc.f:148
subroutine dtrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
DTRSEN
Definition: dtrsen.f:313
subroutine dtrsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
DTRSYL
Definition: dtrsyl.f:164