LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztrsen.f
Go to the documentation of this file.
1*> \brief \b ZTRSEN
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZTRSEN + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsen.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsen.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsen.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
22* SEP, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER COMPQ, JOB
26* INTEGER INFO, LDQ, LDT, LWORK, M, N
27* DOUBLE PRECISION S, SEP
28* ..
29* .. Array Arguments ..
30* LOGICAL SELECT( * )
31* COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> ZTRSEN reorders the Schur factorization of a complex matrix
41*> A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
42*> the leading positions on the diagonal of the upper triangular matrix
43*> T, and the leading columns of Q form an orthonormal basis of the
44*> corresponding right invariant subspace.
45*>
46*> Optionally the routine computes the reciprocal condition numbers of
47*> the cluster of eigenvalues and/or the invariant subspace.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] JOB
54*> \verbatim
55*> JOB is CHARACTER*1
56*> Specifies whether condition numbers are required for the
57*> cluster of eigenvalues (S) or the invariant subspace (SEP):
58*> = 'N': none;
59*> = 'E': for eigenvalues only (S);
60*> = 'V': for invariant subspace only (SEP);
61*> = 'B': for both eigenvalues and invariant subspace (S and
62*> SEP).
63*> \endverbatim
64*>
65*> \param[in] COMPQ
66*> \verbatim
67*> COMPQ is CHARACTER*1
68*> = 'V': update the matrix Q of Schur vectors;
69*> = 'N': do not update Q.
70*> \endverbatim
71*>
72*> \param[in] SELECT
73*> \verbatim
74*> SELECT is LOGICAL array, dimension (N)
75*> SELECT specifies the eigenvalues in the selected cluster. To
76*> select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
77*> \endverbatim
78*>
79*> \param[in] N
80*> \verbatim
81*> N is INTEGER
82*> The order of the matrix T. N >= 0.
83*> \endverbatim
84*>
85*> \param[in,out] T
86*> \verbatim
87*> T is COMPLEX*16 array, dimension (LDT,N)
88*> On entry, the upper triangular matrix T.
89*> On exit, T is overwritten by the reordered matrix T, with the
90*> selected eigenvalues as the leading diagonal elements.
91*> \endverbatim
92*>
93*> \param[in] LDT
94*> \verbatim
95*> LDT is INTEGER
96*> The leading dimension of the array T. LDT >= max(1,N).
97*> \endverbatim
98*>
99*> \param[in,out] Q
100*> \verbatim
101*> Q is COMPLEX*16 array, dimension (LDQ,N)
102*> On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
103*> On exit, if COMPQ = 'V', Q has been postmultiplied by the
104*> unitary transformation matrix which reorders T; the leading M
105*> columns of Q form an orthonormal basis for the specified
106*> invariant subspace.
107*> If COMPQ = 'N', Q is not referenced.
108*> \endverbatim
109*>
110*> \param[in] LDQ
111*> \verbatim
112*> LDQ is INTEGER
113*> The leading dimension of the array Q.
114*> LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
115*> \endverbatim
116*>
117*> \param[out] W
118*> \verbatim
119*> W is COMPLEX*16 array, dimension (N)
120*> The reordered eigenvalues of T, in the same order as they
121*> appear on the diagonal of T.
122*> \endverbatim
123*>
124*> \param[out] M
125*> \verbatim
126*> M is INTEGER
127*> The dimension of the specified invariant subspace.
128*> 0 <= M <= N.
129*> \endverbatim
130*>
131*> \param[out] S
132*> \verbatim
133*> S is DOUBLE PRECISION
134*> If JOB = 'E' or 'B', S is a lower bound on the reciprocal
135*> condition number for the selected cluster of eigenvalues.
136*> S cannot underestimate the true reciprocal condition number
137*> by more than a factor of sqrt(N). If M = 0 or N, S = 1.
138*> If JOB = 'N' or 'V', S is not referenced.
139*> \endverbatim
140*>
141*> \param[out] SEP
142*> \verbatim
143*> SEP is DOUBLE PRECISION
144*> If JOB = 'V' or 'B', SEP is the estimated reciprocal
145*> condition number of the specified invariant subspace. If
146*> M = 0 or N, SEP = norm(T).
147*> If JOB = 'N' or 'E', SEP is not referenced.
148*> \endverbatim
149*>
150*> \param[out] WORK
151*> \verbatim
152*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
153*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
154*> \endverbatim
155*>
156*> \param[in] LWORK
157*> \verbatim
158*> LWORK is INTEGER
159*> The dimension of the array WORK.
160*> If JOB = 'N', LWORK >= 1;
161*> if JOB = 'E', LWORK = max(1,M*(N-M));
162*> if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
163*>
164*> If LWORK = -1, then a workspace query is assumed; the routine
165*> only calculates the optimal size of the WORK array, returns
166*> this value as the first entry of the WORK array, and no error
167*> message related to LWORK is issued by XERBLA.
168*> \endverbatim
169*>
170*> \param[out] INFO
171*> \verbatim
172*> INFO is INTEGER
173*> = 0: successful exit
174*> < 0: if INFO = -i, the i-th argument had an illegal value
175*> \endverbatim
176*
177* Authors:
178* ========
179*
180*> \author Univ. of Tennessee
181*> \author Univ. of California Berkeley
182*> \author Univ. of Colorado Denver
183*> \author NAG Ltd.
184*
185*> \ingroup trsen
186*
187*> \par Further Details:
188* =====================
189*>
190*> \verbatim
191*>
192*> ZTRSEN first collects the selected eigenvalues by computing a unitary
193*> transformation Z to move them to the top left corner of T. In other
194*> words, the selected eigenvalues are the eigenvalues of T11 in:
195*>
196*> Z**H * T * Z = ( T11 T12 ) n1
197*> ( 0 T22 ) n2
198*> n1 n2
199*>
200*> where N = n1+n2. The first
201*> n1 columns of Z span the specified invariant subspace of T.
202*>
203*> If T has been obtained from the Schur factorization of a matrix
204*> A = Q*T*Q**H, then the reordered Schur factorization of A is given by
205*> A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
206*> corresponding invariant subspace of A.
207*>
208*> The reciprocal condition number of the average of the eigenvalues of
209*> T11 may be returned in S. S lies between 0 (very badly conditioned)
210*> and 1 (very well conditioned). It is computed as follows. First we
211*> compute R so that
212*>
213*> P = ( I R ) n1
214*> ( 0 0 ) n2
215*> n1 n2
216*>
217*> is the projector on the invariant subspace associated with T11.
218*> R is the solution of the Sylvester equation:
219*>
220*> T11*R - R*T22 = T12.
221*>
222*> Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
223*> the two-norm of M. Then S is computed as the lower bound
224*>
225*> (1 + F-norm(R)**2)**(-1/2)
226*>
227*> on the reciprocal of 2-norm(P), the true reciprocal condition number.
228*> S cannot underestimate 1 / 2-norm(P) by more than a factor of
229*> sqrt(N).
230*>
231*> An approximate error bound for the computed average of the
232*> eigenvalues of T11 is
233*>
234*> EPS * norm(T) / S
235*>
236*> where EPS is the machine precision.
237*>
238*> The reciprocal condition number of the right invariant subspace
239*> spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
240*> SEP is defined as the separation of T11 and T22:
241*>
242*> sep( T11, T22 ) = sigma-min( C )
243*>
244*> where sigma-min(C) is the smallest singular value of the
245*> n1*n2-by-n1*n2 matrix
246*>
247*> C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
248*>
249*> I(m) is an m by m identity matrix, and kprod denotes the Kronecker
250*> product. We estimate sigma-min(C) by the reciprocal of an estimate of
251*> the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
252*> cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
253*>
254*> When SEP is small, small changes in T can cause large changes in
255*> the invariant subspace. An approximate bound on the maximum angular
256*> error in the computed right invariant subspace is
257*>
258*> EPS * norm(T) / SEP
259*> \endverbatim
260*>
261* =====================================================================
262 SUBROUTINE ztrsen( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
263 $ SEP, WORK, LWORK, INFO )
264*
265* -- LAPACK computational routine --
266* -- LAPACK is a software package provided by Univ. of Tennessee, --
267* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268*
269* .. Scalar Arguments ..
270 CHARACTER COMPQ, JOB
271 INTEGER INFO, LDQ, LDT, LWORK, M, N
272 DOUBLE PRECISION S, SEP
273* ..
274* .. Array Arguments ..
275 LOGICAL SELECT( * )
276 COMPLEX*16 Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
277* ..
278*
279* =====================================================================
280*
281* .. Parameters ..
282 DOUBLE PRECISION ZERO, ONE
283 parameter( zero = 0.0d+0, one = 1.0d+0 )
284* ..
285* .. Local Scalars ..
286 LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
287 INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
288 DOUBLE PRECISION EST, RNORM, SCALE
289* ..
290* .. Local Arrays ..
291 INTEGER ISAVE( 3 )
292 DOUBLE PRECISION RWORK( 1 )
293* ..
294* .. External Functions ..
295 LOGICAL LSAME
296 DOUBLE PRECISION ZLANGE
297 EXTERNAL lsame, zlange
298* ..
299* .. External Subroutines ..
300 EXTERNAL xerbla, zlacn2, zlacpy, ztrexc, ztrsyl
301* ..
302* .. Intrinsic Functions ..
303 INTRINSIC max, sqrt
304* ..
305* .. Executable Statements ..
306*
307* Decode and test the input parameters.
308*
309 wantbh = lsame( job, 'B' )
310 wants = lsame( job, 'E' ) .OR. wantbh
311 wantsp = lsame( job, 'V' ) .OR. wantbh
312 wantq = lsame( compq, 'V' )
313*
314* Set M to the number of selected eigenvalues.
315*
316 m = 0
317 DO 10 k = 1, n
318 IF( SELECT( k ) )
319 $ m = m + 1
320 10 CONTINUE
321*
322 n1 = m
323 n2 = n - m
324 nn = n1*n2
325*
326 info = 0
327 lquery = ( lwork.EQ.-1 )
328*
329 IF( wantsp ) THEN
330 lwmin = max( 1, 2*nn )
331 ELSE IF( lsame( job, 'N' ) ) THEN
332 lwmin = 1
333 ELSE IF( lsame( job, 'E' ) ) THEN
334 lwmin = max( 1, nn )
335 END IF
336*
337 IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
338 $ THEN
339 info = -1
340 ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
341 info = -2
342 ELSE IF( n.LT.0 ) THEN
343 info = -4
344 ELSE IF( ldt.LT.max( 1, n ) ) THEN
345 info = -6
346 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
347 info = -8
348 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
349 info = -14
350 END IF
351*
352 IF( info.EQ.0 ) THEN
353 work( 1 ) = lwmin
354 END IF
355*
356 IF( info.NE.0 ) THEN
357 CALL xerbla( 'ZTRSEN', -info )
358 RETURN
359 ELSE IF( lquery ) THEN
360 RETURN
361 END IF
362*
363* Quick return if possible
364*
365 IF( m.EQ.n .OR. m.EQ.0 ) THEN
366 IF( wants )
367 $ s = one
368 IF( wantsp )
369 $ sep = zlange( '1', n, n, t, ldt, rwork )
370 GO TO 40
371 END IF
372*
373* Collect the selected eigenvalues at the top left corner of T.
374*
375 ks = 0
376 DO 20 k = 1, n
377 IF( SELECT( k ) ) THEN
378 ks = ks + 1
379*
380* Swap the K-th eigenvalue to position KS.
381*
382 IF( k.NE.ks )
383 $ CALL ztrexc( compq, n, t, ldt, q, ldq, k, ks, ierr )
384 END IF
385 20 CONTINUE
386*
387 IF( wants ) THEN
388*
389* Solve the Sylvester equation for R:
390*
391* T11*R - R*T22 = scale*T12
392*
393 CALL zlacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
394 CALL ztrsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
395 $ ldt, work, n1, scale, ierr )
396*
397* Estimate the reciprocal of the condition number of the cluster
398* of eigenvalues.
399*
400 rnorm = zlange( 'F', n1, n2, work, n1, rwork )
401 IF( rnorm.EQ.zero ) THEN
402 s = one
403 ELSE
404 s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
405 $ sqrt( rnorm ) )
406 END IF
407 END IF
408*
409 IF( wantsp ) THEN
410*
411* Estimate sep(T11,T22).
412*
413 est = zero
414 kase = 0
415 30 CONTINUE
416 CALL zlacn2( nn, work( nn+1 ), work, est, kase, isave )
417 IF( kase.NE.0 ) THEN
418 IF( kase.EQ.1 ) THEN
419*
420* Solve T11*R - R*T22 = scale*X.
421*
422 CALL ztrsyl( 'N', 'N', -1, n1, n2, t, ldt,
423 $ t( n1+1, n1+1 ), ldt, work, n1, scale,
424 $ ierr )
425 ELSE
426*
427* Solve T11**H*R - R*T22**H = scale*X.
428*
429 CALL ztrsyl( 'C', 'C', -1, n1, n2, t, ldt,
430 $ t( n1+1, n1+1 ), ldt, work, n1, scale,
431 $ ierr )
432 END IF
433 GO TO 30
434 END IF
435*
436 sep = scale / est
437 END IF
438*
439 40 CONTINUE
440*
441* Copy reordered eigenvalues to W.
442*
443 DO 50 k = 1, n
444 w( k ) = t( k, k )
445 50 CONTINUE
446*
447 work( 1 ) = lwmin
448*
449 RETURN
450*
451* End of ZTRSEN
452*
453 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:133
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
subroutine ztrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
ZTREXC
Definition ztrexc.f:126
subroutine ztrsen(job, compq, select, n, t, ldt, q, ldq, w, m, s, sep, work, lwork, info)
ZTRSEN
Definition ztrsen.f:264
subroutine ztrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
ZTRSYL
Definition ztrsyl.f:157