LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
strsen.f
Go to the documentation of this file.
1 *> \brief \b STRSEN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STRSEN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsen.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsen.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsen.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STRSEN( 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 * REAL S, SEP
28 * ..
29 * .. Array Arguments ..
30 * LOGICAL SELECT( * )
31 * INTEGER IWORK( * )
32 * REAL Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
33 * $ WR( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> STRSEN 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 SHSEQR), 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 REAL 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 REAL 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 REAL array, dimension (N)
134 *> \endverbatim
135 *>
136 *> \param[out] WI
137 *> \verbatim
138 *> WI is REAL array, dimension (N)
139 *>
140 *> The real and imaginary parts, respectively, of the reordered
141 *> eigenvalues of T. The eigenvalues are stored in the same
142 *> order as on the diagonal of T, with WR(i) = T(i,i) and, if
143 *> T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
144 *> WI(i+1) = -WI(i). Note that if a complex eigenvalue is
145 *> sufficiently ill-conditioned, then its value may differ
146 *> significantly from its value before reordering.
147 *> \endverbatim
148 *>
149 *> \param[out] M
150 *> \verbatim
151 *> M is INTEGER
152 *> The dimension of the specified invariant subspace.
153 *> 0 < = M <= N.
154 *> \endverbatim
155 *>
156 *> \param[out] S
157 *> \verbatim
158 *> S is REAL
159 *> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
160 *> condition number for the selected cluster of eigenvalues.
161 *> S cannot underestimate the true reciprocal condition number
162 *> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
163 *> If JOB = 'N' or 'V', S is not referenced.
164 *> \endverbatim
165 *>
166 *> \param[out] SEP
167 *> \verbatim
168 *> SEP is REAL
169 *> If JOB = 'V' or 'B', SEP is the estimated reciprocal
170 *> condition number of the specified invariant subspace. If
171 *> M = 0 or N, SEP = norm(T).
172 *> If JOB = 'N' or 'E', SEP is not referenced.
173 *> \endverbatim
174 *>
175 *> \param[out] WORK
176 *> \verbatim
177 *> WORK is REAL array, dimension (MAX(1,LWORK))
178 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
179 *> \endverbatim
180 *>
181 *> \param[in] LWORK
182 *> \verbatim
183 *> LWORK is INTEGER
184 *> The dimension of the array WORK.
185 *> If JOB = 'N', LWORK >= max(1,N);
186 *> if JOB = 'E', LWORK >= max(1,M*(N-M));
187 *> if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
188 *>
189 *> If LWORK = -1, then a workspace query is assumed; the routine
190 *> only calculates the optimal size of the WORK array, returns
191 *> this value as the first entry of the WORK array, and no error
192 *> message related to LWORK is issued by XERBLA.
193 *> \endverbatim
194 *>
195 *> \param[out] IWORK
196 *> \verbatim
197 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
198 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
199 *> \endverbatim
200 *>
201 *> \param[in] LIWORK
202 *> \verbatim
203 *> LIWORK is INTEGER
204 *> The dimension of the array IWORK.
205 *> If JOB = 'N' or 'E', LIWORK >= 1;
206 *> if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
207 *>
208 *> If LIWORK = -1, then a workspace query is assumed; the
209 *> routine only calculates the optimal size of the IWORK array,
210 *> returns this value as the first entry of the IWORK array, and
211 *> no error message related to LIWORK is issued by XERBLA.
212 *> \endverbatim
213 *>
214 *> \param[out] INFO
215 *> \verbatim
216 *> INFO is INTEGER
217 *> = 0: successful exit
218 *> < 0: if INFO = -i, the i-th argument had an illegal value
219 *> = 1: reordering of T failed because some eigenvalues are too
220 *> close to separate (the problem is very ill-conditioned);
221 *> T may have been partially reordered, and WR and WI
222 *> contain the eigenvalues in the same order as in T; S and
223 *> SEP (if requested) are set to zero.
224 *> \endverbatim
225 *
226 * Authors:
227 * ========
228 *
229 *> \author Univ. of Tennessee
230 *> \author Univ. of California Berkeley
231 *> \author Univ. of Colorado Denver
232 *> \author NAG Ltd.
233 *
234 *> \ingroup realOTHERcomputational
235 *
236 *> \par Further Details:
237 * =====================
238 *>
239 *> \verbatim
240 *>
241 *> STRSEN first collects the selected eigenvalues by computing an
242 *> orthogonal transformation Z to move them to the top left corner of T.
243 *> In other words, the selected eigenvalues are the eigenvalues of T11
244 *> in:
245 *>
246 *> Z**T * T * Z = ( T11 T12 ) n1
247 *> ( 0 T22 ) n2
248 *> n1 n2
249 *>
250 *> where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
251 *> of Z span the specified invariant subspace of T.
252 *>
253 *> If T has been obtained from the real Schur factorization of a matrix
254 *> A = Q*T*Q**T, then the reordered real Schur factorization of A is given
255 *> by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span
256 *> the corresponding invariant subspace of A.
257 *>
258 *> The reciprocal condition number of the average of the eigenvalues of
259 *> T11 may be returned in S. S lies between 0 (very badly conditioned)
260 *> and 1 (very well conditioned). It is computed as follows. First we
261 *> compute R so that
262 *>
263 *> P = ( I R ) n1
264 *> ( 0 0 ) n2
265 *> n1 n2
266 *>
267 *> is the projector on the invariant subspace associated with T11.
268 *> R is the solution of the Sylvester equation:
269 *>
270 *> T11*R - R*T22 = T12.
271 *>
272 *> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
273 *> the two-norm of M. Then S is computed as the lower bound
274 *>
275 *> (1 + F-norm(R)**2)**(-1/2)
276 *>
277 *> on the reciprocal of 2-norm(P), the true reciprocal condition number.
278 *> S cannot underestimate 1 / 2-norm(P) by more than a factor of
279 *> sqrt(N).
280 *>
281 *> An approximate error bound for the computed average of the
282 *> eigenvalues of T11 is
283 *>
284 *> EPS * norm(T) / S
285 *>
286 *> where EPS is the machine precision.
287 *>
288 *> The reciprocal condition number of the right invariant subspace
289 *> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
290 *> SEP is defined as the separation of T11 and T22:
291 *>
292 *> sep( T11, T22 ) = sigma-min( C )
293 *>
294 *> where sigma-min(C) is the smallest singular value of the
295 *> n1*n2-by-n1*n2 matrix
296 *>
297 *> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
298 *>
299 *> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
300 *> product. We estimate sigma-min(C) by the reciprocal of an estimate of
301 *> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
302 *> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
303 *>
304 *> When SEP is small, small changes in T can cause large changes in
305 *> the invariant subspace. An approximate bound on the maximum angular
306 *> error in the computed right invariant subspace is
307 *>
308 *> EPS * norm(T) / SEP
309 *> \endverbatim
310 *>
311 * =====================================================================
312  SUBROUTINE strsen( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
313  $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
314 *
315 * -- LAPACK computational routine --
316 * -- LAPACK is a software package provided by Univ. of Tennessee, --
317 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
318 *
319 * .. Scalar Arguments ..
320  CHARACTER COMPQ, JOB
321  INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
322  REAL S, SEP
323 * ..
324 * .. Array Arguments ..
325  LOGICAL SELECT( * )
326  INTEGER IWORK( * )
327  REAL Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
328  $ wr( * )
329 * ..
330 *
331 * =====================================================================
332 *
333 * .. Parameters ..
334  REAL ZERO, ONE
335  parameter( zero = 0.0e+0, one = 1.0e+0 )
336 * ..
337 * .. Local Scalars ..
338  LOGICAL LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS,
339  $ wantsp
340  INTEGER IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2,
341  $ nn
342  REAL EST, RNORM, SCALE
343 * ..
344 * .. Local Arrays ..
345  INTEGER ISAVE( 3 )
346 * ..
347 * .. External Functions ..
348  LOGICAL LSAME
349  REAL SLANGE
350  EXTERNAL lsame, slange
351 * ..
352 * .. External Subroutines ..
353  EXTERNAL slacn2, slacpy, strexc, strsyl, xerbla
354 * ..
355 * .. Intrinsic Functions ..
356  INTRINSIC abs, max, sqrt
357 * ..
358 * .. Executable Statements ..
359 *
360 * Decode and test the input parameters
361 *
362  wantbh = lsame( job, 'B' )
363  wants = lsame( job, 'E' ) .OR. wantbh
364  wantsp = lsame( job, 'V' ) .OR. wantbh
365  wantq = lsame( compq, 'V' )
366 *
367  info = 0
368  lquery = ( lwork.EQ.-1 )
369  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
370  $ THEN
371  info = -1
372  ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
373  info = -2
374  ELSE IF( n.LT.0 ) THEN
375  info = -4
376  ELSE IF( ldt.LT.max( 1, n ) ) THEN
377  info = -6
378  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
379  info = -8
380  ELSE
381 *
382 * Set M to the dimension of the specified invariant subspace,
383 * and test LWORK and LIWORK.
384 *
385  m = 0
386  pair = .false.
387  DO 10 k = 1, n
388  IF( pair ) THEN
389  pair = .false.
390  ELSE
391  IF( k.LT.n ) THEN
392  IF( t( k+1, k ).EQ.zero ) THEN
393  IF( SELECT( k ) )
394  $ m = m + 1
395  ELSE
396  pair = .true.
397  IF( SELECT( k ) .OR. SELECT( k+1 ) )
398  $ m = m + 2
399  END IF
400  ELSE
401  IF( SELECT( n ) )
402  $ m = m + 1
403  END IF
404  END IF
405  10 CONTINUE
406 *
407  n1 = m
408  n2 = n - m
409  nn = n1*n2
410 *
411  IF( wantsp ) THEN
412  lwmin = max( 1, 2*nn )
413  liwmin = max( 1, nn )
414  ELSE IF( lsame( job, 'N' ) ) THEN
415  lwmin = max( 1, n )
416  liwmin = 1
417  ELSE IF( lsame( job, 'E' ) ) THEN
418  lwmin = max( 1, nn )
419  liwmin = 1
420  END IF
421 *
422  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
423  info = -15
424  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
425  info = -17
426  END IF
427  END IF
428 *
429  IF( info.EQ.0 ) THEN
430  work( 1 ) = lwmin
431  iwork( 1 ) = liwmin
432  END IF
433 *
434  IF( info.NE.0 ) THEN
435  CALL xerbla( 'STRSEN', -info )
436  RETURN
437  ELSE IF( lquery ) THEN
438  RETURN
439  END IF
440 *
441 * Quick return if possible.
442 *
443  IF( m.EQ.n .OR. m.EQ.0 ) THEN
444  IF( wants )
445  $ s = one
446  IF( wantsp )
447  $ sep = slange( '1', n, n, t, ldt, work )
448  GO TO 40
449  END IF
450 *
451 * Collect the selected blocks at the top-left corner of T.
452 *
453  ks = 0
454  pair = .false.
455  DO 20 k = 1, n
456  IF( pair ) THEN
457  pair = .false.
458  ELSE
459  swap = SELECT( k )
460  IF( k.LT.n ) THEN
461  IF( t( k+1, k ).NE.zero ) THEN
462  pair = .true.
463  swap = swap .OR. SELECT( k+1 )
464  END IF
465  END IF
466  IF( swap ) THEN
467  ks = ks + 1
468 *
469 * Swap the K-th block to position KS.
470 *
471  ierr = 0
472  kk = k
473  IF( k.NE.ks )
474  $ CALL strexc( compq, n, t, ldt, q, ldq, kk, ks, work,
475  $ ierr )
476  IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
477 *
478 * Blocks too close to swap: exit.
479 *
480  info = 1
481  IF( wants )
482  $ s = zero
483  IF( wantsp )
484  $ sep = zero
485  GO TO 40
486  END IF
487  IF( pair )
488  $ ks = ks + 1
489  END IF
490  END IF
491  20 CONTINUE
492 *
493  IF( wants ) THEN
494 *
495 * Solve Sylvester equation for R:
496 *
497 * T11*R - R*T22 = scale*T12
498 *
499  CALL slacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
500  CALL strsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
501  $ ldt, work, n1, scale, ierr )
502 *
503 * Estimate the reciprocal of the condition number of the cluster
504 * of eigenvalues.
505 *
506  rnorm = slange( 'F', n1, n2, work, n1, work )
507  IF( rnorm.EQ.zero ) THEN
508  s = one
509  ELSE
510  s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
511  $ sqrt( rnorm ) )
512  END IF
513  END IF
514 *
515  IF( wantsp ) THEN
516 *
517 * Estimate sep(T11,T22).
518 *
519  est = zero
520  kase = 0
521  30 CONTINUE
522  CALL slacn2( nn, work( nn+1 ), work, iwork, est, kase, isave )
523  IF( kase.NE.0 ) THEN
524  IF( kase.EQ.1 ) THEN
525 *
526 * Solve T11*R - R*T22 = scale*X.
527 *
528  CALL strsyl( 'N', 'N', -1, n1, n2, t, ldt,
529  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
530  $ ierr )
531  ELSE
532 *
533 * Solve T11**T*R - R*T22**T = scale*X.
534 *
535  CALL strsyl( 'T', 'T', -1, n1, n2, t, ldt,
536  $ t( n1+1, n1+1 ), ldt, work, n1, scale,
537  $ ierr )
538  END IF
539  GO TO 30
540  END IF
541 *
542  sep = scale / est
543  END IF
544 *
545  40 CONTINUE
546 *
547 * Store the output eigenvalues in WR and WI.
548 *
549  DO 50 k = 1, n
550  wr( k ) = t( k, k )
551  wi( k ) = zero
552  50 CONTINUE
553  DO 60 k = 1, n - 1
554  IF( t( k+1, k ).NE.zero ) THEN
555  wi( k ) = sqrt( abs( t( k, k+1 ) ) )*
556  $ sqrt( abs( t( k+1, k ) ) )
557  wi( k+1 ) = -wi( k )
558  END IF
559  60 CONTINUE
560 *
561  work( 1 ) = lwmin
562  iwork( 1 ) = liwmin
563 *
564  RETURN
565 *
566 * End of STRSEN
567 *
568  END
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:136
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
Definition: strexc.f:148
subroutine strsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
STRSEN
Definition: strsen.f:314
subroutine strsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
STRSYL
Definition: strsyl.f:164