LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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