LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctgsen.f
Go to the documentation of this file.
1*> \brief \b CTGSEN
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTGSEN + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsen.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsen.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsen.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
22* ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
23* WORK, LWORK, IWORK, LIWORK, INFO )
24*
25* .. Scalar Arguments ..
26* LOGICAL WANTQ, WANTZ
27* INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
28* $ M, N
29* REAL PL, PR
30* ..
31* .. Array Arguments ..
32* LOGICAL SELECT( * )
33* INTEGER IWORK( * )
34* REAL DIF( * )
35* COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
36* $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
37* ..
38*
39*
40*> \par Purpose:
41* =============
42*>
43*> \verbatim
44*>
45*> CTGSEN reorders the generalized Schur decomposition of a complex
46*> matrix pair (A, B) (in terms of an unitary equivalence trans-
47*> formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
48*> appears in the leading diagonal blocks of the pair (A,B). The leading
49*> columns of Q and Z form unitary bases of the corresponding left and
50*> right eigenspaces (deflating subspaces). (A, B) must be in
51*> generalized Schur canonical form, that is, A and B are both upper
52*> triangular.
53*>
54*> CTGSEN also computes the generalized eigenvalues
55*>
56*> w(j)= ALPHA(j) / BETA(j)
57*>
58*> of the reordered matrix pair (A, B).
59*>
60*> Optionally, the routine computes estimates of reciprocal condition
61*> numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
62*> (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
63*> between the matrix pairs (A11, B11) and (A22,B22) that correspond to
64*> the selected cluster and the eigenvalues outside the cluster, resp.,
65*> and norms of "projections" onto left and right eigenspaces w.r.t.
66*> the selected cluster in the (1,1)-block.
67*>
68*> \endverbatim
69*
70* Arguments:
71* ==========
72*
73*> \param[in] IJOB
74*> \verbatim
75*> IJOB is INTEGER
76*> Specifies whether condition numbers are required for the
77*> cluster of eigenvalues (PL and PR) or the deflating subspaces
78*> (Difu and Difl):
79*> =0: Only reorder w.r.t. SELECT. No extras.
80*> =1: Reciprocal of norms of "projections" onto left and right
81*> eigenspaces w.r.t. the selected cluster (PL and PR).
82*> =2: Upper bounds on Difu and Difl. F-norm-based estimate
83*> (DIF(1:2)).
84*> =3: Estimate of Difu and Difl. 1-norm-based estimate
85*> (DIF(1:2)).
86*> About 5 times as expensive as IJOB = 2.
87*> =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
88*> version to get it all.
89*> =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
90*> \endverbatim
91*>
92*> \param[in] WANTQ
93*> \verbatim
94*> WANTQ is LOGICAL
95*> .TRUE. : update the left transformation matrix Q;
96*> .FALSE.: do not update Q.
97*> \endverbatim
98*>
99*> \param[in] WANTZ
100*> \verbatim
101*> WANTZ is LOGICAL
102*> .TRUE. : update the right transformation matrix Z;
103*> .FALSE.: do not update Z.
104*> \endverbatim
105*>
106*> \param[in] SELECT
107*> \verbatim
108*> SELECT is LOGICAL array, dimension (N)
109*> SELECT specifies the eigenvalues in the selected cluster. To
110*> select an eigenvalue w(j), SELECT(j) must be set to
111*> .TRUE..
112*> \endverbatim
113*>
114*> \param[in] N
115*> \verbatim
116*> N is INTEGER
117*> The order of the matrices A and B. N >= 0.
118*> \endverbatim
119*>
120*> \param[in,out] A
121*> \verbatim
122*> A is COMPLEX array, dimension(LDA,N)
123*> On entry, the upper triangular matrix A, in generalized
124*> Schur canonical form.
125*> On exit, A is overwritten by the reordered matrix A.
126*> \endverbatim
127*>
128*> \param[in] LDA
129*> \verbatim
130*> LDA is INTEGER
131*> The leading dimension of the array A. LDA >= max(1,N).
132*> \endverbatim
133*>
134*> \param[in,out] B
135*> \verbatim
136*> B is COMPLEX array, dimension(LDB,N)
137*> On entry, the upper triangular matrix B, in generalized
138*> Schur canonical form.
139*> On exit, B is overwritten by the reordered matrix B.
140*> \endverbatim
141*>
142*> \param[in] LDB
143*> \verbatim
144*> LDB is INTEGER
145*> The leading dimension of the array B. LDB >= max(1,N).
146*> \endverbatim
147*>
148*> \param[out] ALPHA
149*> \verbatim
150*> ALPHA is COMPLEX array, dimension (N)
151*> \endverbatim
152*>
153*> \param[out] BETA
154*> \verbatim
155*> BETA is COMPLEX array, dimension (N)
156*>
157*> The diagonal elements of A and B, respectively,
158*> when the pair (A,B) has been reduced to generalized Schur
159*> form. ALPHA(i)/BETA(i) i=1,...,N are the generalized
160*> eigenvalues.
161*> \endverbatim
162*>
163*> \param[in,out] Q
164*> \verbatim
165*> Q is COMPLEX array, dimension (LDQ,N)
166*> On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
167*> On exit, Q has been postmultiplied by the left unitary
168*> transformation matrix which reorder (A, B); The leading M
169*> columns of Q form orthonormal bases for the specified pair of
170*> left eigenspaces (deflating subspaces).
171*> If WANTQ = .FALSE., Q is not referenced.
172*> \endverbatim
173*>
174*> \param[in] LDQ
175*> \verbatim
176*> LDQ is INTEGER
177*> The leading dimension of the array Q. LDQ >= 1.
178*> If WANTQ = .TRUE., LDQ >= N.
179*> \endverbatim
180*>
181*> \param[in,out] Z
182*> \verbatim
183*> Z is COMPLEX array, dimension (LDZ,N)
184*> On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
185*> On exit, Z has been postmultiplied by the left unitary
186*> transformation matrix which reorder (A, B); The leading M
187*> columns of Z form orthonormal bases for the specified pair of
188*> left eigenspaces (deflating subspaces).
189*> If WANTZ = .FALSE., Z is not referenced.
190*> \endverbatim
191*>
192*> \param[in] LDZ
193*> \verbatim
194*> LDZ is INTEGER
195*> The leading dimension of the array Z. LDZ >= 1.
196*> If WANTZ = .TRUE., LDZ >= N.
197*> \endverbatim
198*>
199*> \param[out] M
200*> \verbatim
201*> M is INTEGER
202*> The dimension of the specified pair of left and right
203*> eigenspaces, (deflating subspaces) 0 <= M <= N.
204*> \endverbatim
205*>
206*> \param[out] PL
207*> \verbatim
208*> PL is REAL
209*> \endverbatim
210*>
211*> \param[out] PR
212*> \verbatim
213*> PR is REAL
214*>
215*> If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
216*> reciprocal of the norm of "projections" onto left and right
217*> eigenspace with respect to the selected cluster.
218*> 0 < PL, PR <= 1.
219*> If M = 0 or M = N, PL = PR = 1.
220*> If IJOB = 0, 2 or 3 PL, PR are not referenced.
221*> \endverbatim
222*>
223*> \param[out] DIF
224*> \verbatim
225*> DIF is REAL array, dimension (2).
226*> If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
227*> If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
228*> Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
229*> estimates of Difu and Difl, computed using reversed
230*> communication with CLACN2.
231*> If M = 0 or N, DIF(1:2) = F-norm([A, B]).
232*> If IJOB = 0 or 1, DIF is not referenced.
233*> \endverbatim
234*>
235*> \param[out] WORK
236*> \verbatim
237*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
238*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
239*> \endverbatim
240*>
241*> \param[in] LWORK
242*> \verbatim
243*> LWORK is INTEGER
244*> The dimension of the array WORK. LWORK >= 1
245*> If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M)
246*> If IJOB = 3 or 5, LWORK >= 4*M*(N-M)
247*>
248*> If LWORK = -1, then a workspace query is assumed; the routine
249*> only calculates the optimal size of the WORK array, returns
250*> this value as the first entry of the WORK array, and no error
251*> message related to LWORK is issued by XERBLA.
252*> \endverbatim
253*>
254*> \param[out] IWORK
255*> \verbatim
256*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
257*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
258*> \endverbatim
259*>
260*> \param[in] LIWORK
261*> \verbatim
262*> LIWORK is INTEGER
263*> The dimension of the array IWORK. LIWORK >= 1.
264*> If IJOB = 1, 2 or 4, LIWORK >= N+2;
265*> If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M));
266*>
267*> If LIWORK = -1, then a workspace query is assumed; the
268*> routine only calculates the optimal size of the IWORK array,
269*> returns this value as the first entry of the IWORK array, and
270*> no error message related to LIWORK is issued by XERBLA.
271*> \endverbatim
272*>
273*> \param[out] INFO
274*> \verbatim
275*> INFO is INTEGER
276*> =0: Successful exit.
277*> <0: If INFO = -i, the i-th argument had an illegal value.
278*> =1: Reordering of (A, B) failed because the transformed
279*> matrix pair (A, B) would be too far from generalized
280*> Schur form; the problem is very ill-conditioned.
281*> (A, B) may have been partially reordered.
282*> If requested, 0 is returned in DIF(*), PL and PR.
283*> \endverbatim
284*
285* Authors:
286* ========
287*
288*> \author Univ. of Tennessee
289*> \author Univ. of California Berkeley
290*> \author Univ. of Colorado Denver
291*> \author NAG Ltd.
292*
293*> \ingroup tgsen
294*
295*> \par Further Details:
296* =====================
297*>
298*> \verbatim
299*>
300*> CTGSEN first collects the selected eigenvalues by computing unitary
301*> U and W that move them to the top left corner of (A, B). In other
302*> words, the selected eigenvalues are the eigenvalues of (A11, B11) in
303*>
304*> U**H*(A, B)*W = (A11 A12) (B11 B12) n1
305*> ( 0 A22),( 0 B22) n2
306*> n1 n2 n1 n2
307*>
308*> where N = n1+n2 and U**H means the conjugate transpose of U. The first
309*> n1 columns of U and W span the specified pair of left and right
310*> eigenspaces (deflating subspaces) of (A, B).
311*>
312*> If (A, B) has been obtained from the generalized real Schur
313*> decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
314*> reordered generalized Schur form of (C, D) is given by
315*>
316*> (C, D) = (Q*U)*(U**H *(A, B)*W)*(Z*W)**H,
317*>
318*> and the first n1 columns of Q*U and Z*W span the corresponding
319*> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
320*>
321*> Note that if the selected eigenvalue is sufficiently ill-conditioned,
322*> then its value may differ significantly from its value before
323*> reordering.
324*>
325*> The reciprocal condition numbers of the left and right eigenspaces
326*> spanned by the first n1 columns of U and W (or Q*U and Z*W) may
327*> be returned in DIF(1:2), corresponding to Difu and Difl, resp.
328*>
329*> The Difu and Difl are defined as:
330*>
331*> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
332*> and
333*> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
334*>
335*> where sigma-min(Zu) is the smallest singular value of the
336*> (2*n1*n2)-by-(2*n1*n2) matrix
337*>
338*> Zu = [ kron(In2, A11) -kron(A22**H, In1) ]
339*> [ kron(In2, B11) -kron(B22**H, In1) ].
340*>
341*> Here, Inx is the identity matrix of size nx and A22**H is the
342*> conjugate transpose of A22. kron(X, Y) is the Kronecker product between
343*> the matrices X and Y.
344*>
345*> When DIF(2) is small, small changes in (A, B) can cause large changes
346*> in the deflating subspace. An approximate (asymptotic) bound on the
347*> maximum angular error in the computed deflating subspaces is
348*>
349*> EPS * norm((A, B)) / DIF(2),
350*>
351*> where EPS is the machine precision.
352*>
353*> The reciprocal norm of the projectors on the left and right
354*> eigenspaces associated with (A11, B11) may be returned in PL and PR.
355*> They are computed as follows. First we compute L and R so that
356*> P*(A, B)*Q is block diagonal, where
357*>
358*> P = ( I -L ) n1 Q = ( I R ) n1
359*> ( 0 I ) n2 and ( 0 I ) n2
360*> n1 n2 n1 n2
361*>
362*> and (L, R) is the solution to the generalized Sylvester equation
363*>
364*> A11*R - L*A22 = -A12
365*> B11*R - L*B22 = -B12
366*>
367*> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
368*> An approximate (asymptotic) bound on the average absolute error of
369*> the selected eigenvalues is
370*>
371*> EPS * norm((A, B)) / PL.
372*>
373*> There are also global error bounds which valid for perturbations up
374*> to a certain restriction: A lower bound (x) on the smallest
375*> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
376*> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
377*> (i.e. (A + E, B + F), is
378*>
379*> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
380*>
381*> An approximate bound on x can be computed from DIF(1:2), PL and PR.
382*>
383*> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
384*> (L', R') and unperturbed (L, R) left and right deflating subspaces
385*> associated with the selected cluster in the (1,1)-blocks can be
386*> bounded as
387*>
388*> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
389*> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
390*>
391*> See LAPACK User's Guide section 4.11 or the following references
392*> for more information.
393*>
394*> Note that if the default method for computing the Frobenius-norm-
395*> based estimate DIF is not wanted (see CLATDF), then the parameter
396*> IDIFJB (see below) should be changed from 3 to 4 (routine CLATDF
397*> (IJOB = 2 will be used)). See CTGSYL for more details.
398*> \endverbatim
399*
400*> \par Contributors:
401* ==================
402*>
403*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
404*> Umea University, S-901 87 Umea, Sweden.
405*
406*> \par References:
407* ================
408*>
409*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
410*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
411*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
412*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
413*> \n
414*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
415*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
416*> Estimation: Theory, Algorithms and Software, Report
417*> UMINF - 94.04, Department of Computing Science, Umea University,
418*> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
419*> To appear in Numerical Algorithms, 1996.
420*> \n
421*> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
422*> for Solving the Generalized Sylvester Equation and Estimating the
423*> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
424*> Department of Computing Science, Umea University, S-901 87 Umea,
425*> Sweden, December 1993, Revised April 1994, Also as LAPACK working
426*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
427*> 1996.
428*>
429* =====================================================================
430 SUBROUTINE ctgsen( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
431 $ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
432 $ WORK, LWORK, IWORK, LIWORK, INFO )
433*
434* -- LAPACK computational routine --
435* -- LAPACK is a software package provided by Univ. of Tennessee, --
436* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
437*
438* .. Scalar Arguments ..
439 LOGICAL WANTQ, WANTZ
440 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
441 $ m, n
442 REAL PL, PR
443* ..
444* .. Array Arguments ..
445 LOGICAL SELECT( * )
446 INTEGER IWORK( * )
447 REAL DIF( * )
448 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
449 $ beta( * ), q( ldq, * ), work( * ), z( ldz, * )
450* ..
451*
452* =====================================================================
453*
454* .. Parameters ..
455 INTEGER IDIFJB
456 PARAMETER ( IDIFJB = 3 )
457 REAL ZERO, ONE
458 parameter( zero = 0.0e+0, one = 1.0e+0 )
459* ..
460* .. Local Scalars ..
461 LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
462 INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
463 $ n1, n2
464 REAL DSCALE, DSUM, RDSCAL, SAFMIN
465 COMPLEX TEMP1, TEMP2
466* ..
467* .. Local Arrays ..
468 INTEGER ISAVE( 3 )
469* ..
470* .. External Functions ..
471 REAL SROUNDUP_LWORK
472 EXTERNAL SROUNDUP_LWORK
473* ..
474* .. External Subroutines ..
475 REAL SLAMCH
476 EXTERNAL CLACN2, CLACPY, CLASSQ, CSCAL, CTGEXC, CTGSYL,
477 $ slamch, xerbla
478* ..
479* .. Intrinsic Functions ..
480 INTRINSIC abs, cmplx, conjg, max, sqrt
481* ..
482* .. Executable Statements ..
483*
484* Decode and test the input parameters
485*
486 info = 0
487 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
488*
489 IF( ijob.LT.0 .OR. ijob.GT.5 ) THEN
490 info = -1
491 ELSE IF( n.LT.0 ) THEN
492 info = -5
493 ELSE IF( lda.LT.max( 1, n ) ) THEN
494 info = -7
495 ELSE IF( ldb.LT.max( 1, n ) ) THEN
496 info = -9
497 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
498 info = -13
499 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
500 info = -15
501 END IF
502*
503 IF( info.NE.0 ) THEN
504 CALL xerbla( 'CTGSEN', -info )
505 RETURN
506 END IF
507*
508 ierr = 0
509*
510 wantp = ijob.EQ.1 .OR. ijob.GE.4
511 wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
512 wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
513 wantd = wantd1 .OR. wantd2
514*
515* Set M to the dimension of the specified pair of deflating
516* subspaces.
517*
518 m = 0
519 IF( .NOT.lquery .OR. ijob.NE.0 ) THEN
520 DO 10 k = 1, n
521 alpha( k ) = a( k, k )
522 beta( k ) = b( k, k )
523 IF( k.LT.n ) THEN
524 IF( SELECT( k ) )
525 $ m = m + 1
526 ELSE
527 IF( SELECT( n ) )
528 $ m = m + 1
529 END IF
530 10 CONTINUE
531 END IF
532*
533 IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
534 lwmin = max( 1, 2*m*(n-m) )
535 liwmin = max( 1, n+2 )
536 ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 ) THEN
537 lwmin = max( 1, 4*m*(n-m) )
538 liwmin = max( 1, 2*m*(n-m), n+2 )
539 ELSE
540 lwmin = 1
541 liwmin = 1
542 END IF
543*
544 work( 1 ) = sroundup_lwork(lwmin)
545 iwork( 1 ) = liwmin
546*
547 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
548 info = -21
549 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
550 info = -23
551 END IF
552*
553 IF( info.NE.0 ) THEN
554 CALL xerbla( 'CTGSEN', -info )
555 RETURN
556 ELSE IF( lquery ) THEN
557 RETURN
558 END IF
559*
560* Quick return if possible.
561*
562 IF( m.EQ.n .OR. m.EQ.0 ) THEN
563 IF( wantp ) THEN
564 pl = one
565 pr = one
566 END IF
567 IF( wantd ) THEN
568 dscale = zero
569 dsum = one
570 DO 20 i = 1, n
571 CALL classq( n, a( 1, i ), 1, dscale, dsum )
572 CALL classq( n, b( 1, i ), 1, dscale, dsum )
573 20 CONTINUE
574 dif( 1 ) = dscale*sqrt( dsum )
575 dif( 2 ) = dif( 1 )
576 END IF
577 GO TO 70
578 END IF
579*
580* Get machine constant
581*
582 safmin = slamch( 'S' )
583*
584* Collect the selected blocks at the top-left corner of (A, B).
585*
586 ks = 0
587 DO 30 k = 1, n
588 swap = SELECT( k )
589 IF( swap ) THEN
590 ks = ks + 1
591*
592* Swap the K-th block to position KS. Compute unitary Q
593* and Z that will swap adjacent diagonal blocks in (A, B).
594*
595 IF( k.NE.ks )
596 $ CALL ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
597 $ ldz, k, ks, ierr )
598*
599 IF( ierr.GT.0 ) THEN
600*
601* Swap is rejected: exit.
602*
603 info = 1
604 IF( wantp ) THEN
605 pl = zero
606 pr = zero
607 END IF
608 IF( wantd ) THEN
609 dif( 1 ) = zero
610 dif( 2 ) = zero
611 END IF
612 GO TO 70
613 END IF
614 END IF
615 30 CONTINUE
616 IF( wantp ) THEN
617*
618* Solve generalized Sylvester equation for R and L:
619* A11 * R - L * A22 = A12
620* B11 * R - L * B22 = B12
621*
622 n1 = m
623 n2 = n - m
624 i = n1 + 1
625 CALL clacpy( 'Full', n1, n2, a( 1, i ), lda, work, n1 )
626 CALL clacpy( 'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
627 $ n1 )
628 ijb = 0
629 CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
630 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
631 $ dscale, dif( 1 ), work( n1*n2*2+1 ),
632 $ lwork-2*n1*n2, iwork, ierr )
633*
634* Estimate the reciprocal of norms of "projections" onto
635* left and right eigenspaces
636*
637 rdscal = zero
638 dsum = one
639 CALL classq( n1*n2, work, 1, rdscal, dsum )
640 pl = rdscal*sqrt( dsum )
641 IF( pl.EQ.zero ) THEN
642 pl = one
643 ELSE
644 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
645 END IF
646 rdscal = zero
647 dsum = one
648 CALL classq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
649 pr = rdscal*sqrt( dsum )
650 IF( pr.EQ.zero ) THEN
651 pr = one
652 ELSE
653 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
654 END IF
655 END IF
656 IF( wantd ) THEN
657*
658* Compute estimates Difu and Difl.
659*
660 IF( wantd1 ) THEN
661 n1 = m
662 n2 = n - m
663 i = n1 + 1
664 ijb = idifjb
665*
666* Frobenius norm-based Difu estimate.
667*
668 CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
669 $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
670 $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
671 $ lwork-2*n1*n2, iwork, ierr )
672*
673* Frobenius norm-based Difl estimate.
674*
675 CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
676 $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
677 $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
678 $ lwork-2*n1*n2, iwork, ierr )
679 ELSE
680*
681* Compute 1-norm-based estimates of Difu and Difl using
682* reversed communication with CLACN2. In each step a
683* generalized Sylvester equation or a transposed variant
684* is solved.
685*
686 kase = 0
687 n1 = m
688 n2 = n - m
689 i = n1 + 1
690 ijb = 0
691 mn2 = 2*n1*n2
692*
693* 1-norm-based estimate of Difu.
694*
695 40 CONTINUE
696 CALL clacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
697 $ isave )
698 IF( kase.NE.0 ) THEN
699 IF( kase.EQ.1 ) THEN
700*
701* Solve generalized Sylvester equation
702*
703 CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,
704 $ work, n1, b, ldb, b( i, i ), ldb,
705 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
706 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
707 $ ierr )
708 ELSE
709*
710* Solve the transposed variant.
711*
712 CALL ctgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,
713 $ work, n1, b, ldb, b( i, i ), ldb,
714 $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
715 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
716 $ ierr )
717 END IF
718 GO TO 40
719 END IF
720 dif( 1 ) = dscale / dif( 1 )
721*
722* 1-norm-based estimate of Difl.
723*
724 50 CONTINUE
725 CALL clacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
726 $ isave )
727 IF( kase.NE.0 ) THEN
728 IF( kase.EQ.1 ) THEN
729*
730* Solve generalized Sylvester equation
731*
732 CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,
733 $ work, n2, b( i, i ), ldb, b, ldb,
734 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
735 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
736 $ ierr )
737 ELSE
738*
739* Solve the transposed variant.
740*
741 CALL ctgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,
742 $ work, n2, b, ldb, b( i, i ), ldb,
743 $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
744 $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
745 $ ierr )
746 END IF
747 GO TO 50
748 END IF
749 dif( 2 ) = dscale / dif( 2 )
750 END IF
751 END IF
752*
753* If B(K,K) is complex, make it real and positive (normalization
754* of the generalized Schur form) and Store the generalized
755* eigenvalues of reordered pair (A, B)
756*
757 DO 60 k = 1, n
758 dscale = abs( b( k, k ) )
759 IF( dscale.GT.safmin ) THEN
760 temp1 = conjg( b( k, k ) / dscale )
761 temp2 = b( k, k ) / dscale
762 b( k, k ) = dscale
763 CALL cscal( n-k, temp1, b( k, k+1 ), ldb )
764 CALL cscal( n-k+1, temp1, a( k, k ), lda )
765 IF( wantq )
766 $ CALL cscal( n, temp2, q( 1, k ), 1 )
767 ELSE
768 b( k, k ) = cmplx( zero, zero )
769 END IF
770*
771 alpha( k ) = a( k, k )
772 beta( k ) = b( k, k )
773*
774 60 CONTINUE
775*
776 70 CONTINUE
777*
778 work( 1 ) = sroundup_lwork(lwmin)
779 iwork( 1 ) = liwmin
780*
781 RETURN
782*
783* End of CTGSEN
784*
785 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ctgsen(ijob, wantq, wantz, select, n, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, m, pl, pr, dif, work, lwork, iwork, liwork, info)
CTGSEN
Definition ctgsen.f:433