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