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