LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dtgex2.f
Go to the documentation of this file.
1*> \brief \b DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgex2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgex2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgex2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
22* LDZ, J1, N1, N2, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* LOGICAL WANTQ, WANTZ
26* INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30* \$ WORK( * ), Z( LDZ, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
40*> of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
41*> (A, B) by an orthogonal equivalence transformation.
42*>
43*> (A, B) must be in generalized real Schur canonical form (as returned
44*> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
45*> diagonal blocks. B is upper triangular.
46*>
47*> Optionally, the matrices Q and Z of generalized Schur vectors are
48*> updated.
49*>
50*> Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
51*> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
52*>
53*> \endverbatim
54*
55* Arguments:
56* ==========
57*
58*> \param[in] WANTQ
59*> \verbatim
60*> WANTQ is LOGICAL
61*> .TRUE. : update the left transformation matrix Q;
62*> .FALSE.: do not update Q.
63*> \endverbatim
64*>
65*> \param[in] WANTZ
66*> \verbatim
67*> WANTZ is LOGICAL
68*> .TRUE. : update the right transformation matrix Z;
69*> .FALSE.: do not update Z.
70*> \endverbatim
71*>
72*> \param[in] N
73*> \verbatim
74*> N is INTEGER
75*> The order of the matrices A and B. N >= 0.
76*> \endverbatim
77*>
78*> \param[in,out] A
79*> \verbatim
80*> A is DOUBLE PRECISION array, dimensions (LDA,N)
81*> On entry, the matrix A in the pair (A, B).
82*> On exit, the updated matrix A.
83*> \endverbatim
84*>
85*> \param[in] LDA
86*> \verbatim
87*> LDA is INTEGER
88*> The leading dimension of the array A. LDA >= max(1,N).
89*> \endverbatim
90*>
91*> \param[in,out] B
92*> \verbatim
93*> B is DOUBLE PRECISION array, dimensions (LDB,N)
94*> On entry, the matrix B in the pair (A, B).
95*> On exit, the updated matrix B.
96*> \endverbatim
97*>
98*> \param[in] LDB
99*> \verbatim
100*> LDB is INTEGER
101*> The leading dimension of the array B. LDB >= max(1,N).
102*> \endverbatim
103*>
104*> \param[in,out] Q
105*> \verbatim
106*> Q is DOUBLE PRECISION array, dimension (LDQ,N)
107*> On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
108*> On exit, the updated matrix Q.
109*> Not referenced if WANTQ = .FALSE..
110*> \endverbatim
111*>
112*> \param[in] LDQ
113*> \verbatim
114*> LDQ is INTEGER
115*> The leading dimension of the array Q. LDQ >= 1.
116*> If WANTQ = .TRUE., LDQ >= N.
117*> \endverbatim
118*>
119*> \param[in,out] Z
120*> \verbatim
121*> Z is DOUBLE PRECISION array, dimension (LDZ,N)
122*> On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
123*> On exit, the updated matrix Z.
124*> Not referenced if WANTZ = .FALSE..
125*> \endverbatim
126*>
127*> \param[in] LDZ
128*> \verbatim
129*> LDZ is INTEGER
130*> The leading dimension of the array Z. LDZ >= 1.
131*> If WANTZ = .TRUE., LDZ >= N.
132*> \endverbatim
133*>
134*> \param[in] J1
135*> \verbatim
136*> J1 is INTEGER
137*> The index to the first block (A11, B11). 1 <= J1 <= N.
138*> \endverbatim
139*>
140*> \param[in] N1
141*> \verbatim
142*> N1 is INTEGER
143*> The order of the first block (A11, B11). N1 = 0, 1 or 2.
144*> \endverbatim
145*>
146*> \param[in] N2
147*> \verbatim
148*> N2 is INTEGER
149*> The order of the second block (A22, B22). N2 = 0, 1 or 2.
150*> \endverbatim
151*>
152*> \param[out] WORK
153*> \verbatim
154*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)).
155*> \endverbatim
156*>
157*> \param[in] LWORK
158*> \verbatim
159*> LWORK is INTEGER
160*> The dimension of the array WORK.
161*> LWORK >= MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 )
162*> \endverbatim
163*>
164*> \param[out] INFO
165*> \verbatim
166*> INFO is INTEGER
167*> =0: Successful exit
168*> >0: If INFO = 1, the transformed matrix (A, B) would be
169*> too far from generalized Schur form; the blocks are
170*> not swapped and (A, B) and (Q, Z) are unchanged.
171*> The problem of swapping is too ill-conditioned.
172*> <0: If INFO = -16: LWORK is too small. Appropriate value
173*> for LWORK is returned in WORK(1).
174*> \endverbatim
175*
176* Authors:
177* ========
178*
179*> \author Univ. of Tennessee
180*> \author Univ. of California Berkeley
181*> \author Univ. of Colorado Denver
182*> \author NAG Ltd.
183*
184*> \ingroup tgex2
185*
186*> \par Further Details:
187* =====================
188*>
189*> In the current code both weak and strong stability tests are
190*> performed. The user can omit the strong stability test by changing
191*> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
192*> details.
193*
194*> \par Contributors:
195* ==================
196*>
197*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
198*> Umea University, S-901 87 Umea, Sweden.
199*
200*> \par References:
201* ================
202*>
203*> \verbatim
204*>
205*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
206*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
207*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
208*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
209*>
210*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
211*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
212*> Estimation: Theory, Algorithms and Software,
213*> Report UMINF - 94.04, Department of Computing Science, Umea
214*> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
215*> Note 87. To appear in Numerical Algorithms, 1996.
216*> \endverbatim
217*>
218* =====================================================================
219 SUBROUTINE dtgex2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
220 \$ LDZ, J1, N1, N2, WORK, LWORK, INFO )
221*
222* -- LAPACK auxiliary routine --
223* -- LAPACK is a software package provided by Univ. of Tennessee, --
224* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225*
226* .. Scalar Arguments ..
227 LOGICAL WANTQ, WANTZ
228 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2
229* ..
230* .. Array Arguments ..
231 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
232 \$ work( * ), z( ldz, * )
233* ..
234*
235* =====================================================================
236* Replaced various illegal calls to DCOPY by calls to DLASET, or by DO
237* loops. Sven Hammarling, 1/5/02.
238*
239* .. Parameters ..
240 DOUBLE PRECISION ZERO, ONE
241 parameter( zero = 0.0d+0, one = 1.0d+0 )
242 DOUBLE PRECISION TWENTY
243 parameter( twenty = 2.0d+01 )
244 INTEGER LDST
245 parameter( ldst = 4 )
246 LOGICAL WANDS
247 parameter( wands = .true. )
248* ..
249* .. Local Scalars ..
250 LOGICAL STRONG, WEAK
251 INTEGER I, IDUM, LINFO, M
252 DOUBLE PRECISION BQRA21, BRQA21, DDUM, DNORMA, DNORMB, DSCALE,
253 \$ dsum, eps, f, g, sa, sb, scale, smlnum,
254 \$ thresha, threshb
255* ..
256* .. Local Arrays ..
257 INTEGER IWORK( LDST + 2 )
258 DOUBLE PRECISION AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
259 \$ ircop( ldst, ldst ), li( ldst, ldst ),
260 \$ licop( ldst, ldst ), s( ldst, ldst ),
261 \$ scpy( ldst, ldst ), t( ldst, ldst ),
262 \$ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
263* ..
264* .. External Functions ..
265 DOUBLE PRECISION DLAMCH
266 EXTERNAL dlamch
267* ..
268* .. External Subroutines ..
269 EXTERNAL dgemm, dgeqr2, dgerq2, dlacpy, dlagv2, dlartg,
271 \$ drot, dscal, dtgsy2
272* ..
273* .. Intrinsic Functions ..
274 INTRINSIC abs, max, sqrt
275* ..
276* .. Executable Statements ..
277*
278 info = 0
279*
280* Quick return if possible
281*
282 IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
283 \$ RETURN
284 IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
285 \$ RETURN
286 m = n1 + n2
287 IF( lwork.LT.max( 1, n*m, m*m*2 ) ) THEN
288 info = -16
289 work( 1 ) = max( 1, n*m, m*m*2 )
290 RETURN
291 END IF
292*
293 weak = .false.
294 strong = .false.
295*
296* Make a local copy of selected block
297*
298 CALL dlaset( 'Full', ldst, ldst, zero, zero, li, ldst )
299 CALL dlaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
300 CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
301 CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
302*
303* Compute threshold for testing acceptance of swapping.
304*
305 eps = dlamch( 'P' )
306 smlnum = dlamch( 'S' ) / eps
307 dscale = zero
308 dsum = one
309 CALL dlacpy( 'Full', m, m, s, ldst, work, m )
310 CALL dlassq( m*m, work, 1, dscale, dsum )
311 dnorma = dscale*sqrt( dsum )
312 dscale = zero
313 dsum = one
314 CALL dlacpy( 'Full', m, m, t, ldst, work, m )
315 CALL dlassq( m*m, work, 1, dscale, dsum )
316 dnormb = dscale*sqrt( dsum )
317*
318* THRES has been changed from
319* THRESH = MAX( TEN*EPS*SA, SMLNUM )
320* to
321* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
322* on 04/01/10.
323* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
324* Jim Demmel and Guillaume Revy. See forum post 1783.
325*
326 thresha = max( twenty*eps*dnorma, smlnum )
327 threshb = max( twenty*eps*dnormb, smlnum )
328*
329 IF( m.EQ.2 ) THEN
330*
331* CASE 1: Swap 1-by-1 and 1-by-1 blocks.
332*
333* Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
334* using Givens rotations and perform the swap tentatively.
335*
336 f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
337 g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
338 sa = abs( s( 2, 2 ) ) * abs( t( 1, 1 ) )
339 sb = abs( s( 1, 1 ) ) * abs( t( 2, 2 ) )
340 CALL dlartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
341 ir( 2, 1 ) = -ir( 1, 2 )
342 ir( 2, 2 ) = ir( 1, 1 )
343 CALL drot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
344 \$ ir( 2, 1 ) )
345 CALL drot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
346 \$ ir( 2, 1 ) )
347 IF( sa.GE.sb ) THEN
348 CALL dlartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
349 \$ ddum )
350 ELSE
351 CALL dlartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
352 \$ ddum )
353 END IF
354 CALL drot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
355 \$ li( 2, 1 ) )
356 CALL drot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
357 \$ li( 2, 1 ) )
358 li( 2, 2 ) = li( 1, 1 )
359 li( 1, 2 ) = -li( 2, 1 )
360*
361* Weak stability test: |S21| <= O(EPS F-norm((A)))
362* and |T21| <= O(EPS F-norm((B)))
363*
364 weak = abs( s( 2, 1 ) ) .LE. thresha .AND.
365 \$ abs( t( 2, 1 ) ) .LE. threshb
366 IF( .NOT.weak )
367 \$ GO TO 70
368*
369 IF( wands ) THEN
370*
371* Strong stability test:
372* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
373* and
374* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
375*
376 CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
377 \$ m )
378 CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
379 \$ work, m )
380 CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
381 \$ work( m*m+1 ), m )
382 dscale = zero
383 dsum = one
384 CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
385 sa = dscale*sqrt( dsum )
386*
387 CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
388 \$ m )
389 CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
390 \$ work, m )
391 CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
392 \$ work( m*m+1 ), m )
393 dscale = zero
394 dsum = one
395 CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
396 sb = dscale*sqrt( dsum )
397 strong = sa.LE.thresha .AND. sb.LE.threshb
398 IF( .NOT.strong )
399 \$ GO TO 70
400 END IF
401*
402* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
403* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
404*
405 CALL drot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
406 \$ ir( 2, 1 ) )
407 CALL drot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
408 \$ ir( 2, 1 ) )
409 CALL drot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
410 \$ li( 1, 1 ), li( 2, 1 ) )
411 CALL drot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
412 \$ li( 1, 1 ), li( 2, 1 ) )
413*
414* Set N1-by-N2 (2,1) - blocks to ZERO.
415*
416 a( j1+1, j1 ) = zero
417 b( j1+1, j1 ) = zero
418*
419* Accumulate transformations into Q and Z if requested.
420*
421 IF( wantz )
422 \$ CALL drot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
423 \$ ir( 2, 1 ) )
424 IF( wantq )
425 \$ CALL drot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
426 \$ li( 2, 1 ) )
427*
428* Exit with INFO = 0 if swap was successfully performed.
429*
430 RETURN
431*
432 ELSE
433*
434* CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
435* and 2-by-2 blocks.
436*
437* Solve the generalized Sylvester equation
438* S11 * R - L * S22 = SCALE * S12
439* T11 * R - L * T22 = SCALE * T12
440* for R and L. Solutions in LI and IR.
441*
442 CALL dlacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
443 CALL dlacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
444 \$ ir( n2+1, n1+1 ), ldst )
445 CALL dtgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
446 \$ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
447 \$ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
448 \$ linfo )
449 IF( linfo.NE.0 )
450 \$ GO TO 70
451*
452* Compute orthogonal matrix QL:
453*
454* QL**T * LI = [ TL ]
455* [ 0 ]
456* where
457* LI = [ -L ]
458* [ SCALE * identity(N2) ]
459*
460 DO 10 i = 1, n2
461 CALL dscal( n1, -one, li( 1, i ), 1 )
462 li( n1+i, i ) = scale
463 10 CONTINUE
464 CALL dgeqr2( m, n2, li, ldst, taul, work, linfo )
465 IF( linfo.NE.0 )
466 \$ GO TO 70
467 CALL dorg2r( m, m, n2, li, ldst, taul, work, linfo )
468 IF( linfo.NE.0 )
469 \$ GO TO 70
470*
471* Compute orthogonal matrix RQ:
472*
473* IR * RQ**T = [ 0 TR],
474*
475* where IR = [ SCALE * identity(N1), R ]
476*
477 DO 20 i = 1, n1
478 ir( n2+i, i ) = scale
479 20 CONTINUE
480 CALL dgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
481 IF( linfo.NE.0 )
482 \$ GO TO 70
483 CALL dorgr2( m, m, n1, ir, ldst, taur, work, linfo )
484 IF( linfo.NE.0 )
485 \$ GO TO 70
486*
487* Perform the swapping tentatively:
488*
489 CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
490 \$ work, m )
491 CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
492 \$ ldst )
493 CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
494 \$ work, m )
495 CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
496 \$ ldst )
497 CALL dlacpy( 'F', m, m, s, ldst, scpy, ldst )
498 CALL dlacpy( 'F', m, m, t, ldst, tcpy, ldst )
499 CALL dlacpy( 'F', m, m, ir, ldst, ircop, ldst )
500 CALL dlacpy( 'F', m, m, li, ldst, licop, ldst )
501*
502* Triangularize the B-part by an RQ factorization.
503* Apply transformation (from left) to A-part, giving S.
504*
505 CALL dgerq2( m, m, t, ldst, taur, work, linfo )
506 IF( linfo.NE.0 )
507 \$ GO TO 70
508 CALL dormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
509 \$ linfo )
510 IF( linfo.NE.0 )
511 \$ GO TO 70
512 CALL dormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst, work,
513 \$ linfo )
514 IF( linfo.NE.0 )
515 \$ GO TO 70
516*
517* Compute F-norm(S21) in BRQA21. (T21 is 0.)
518*
519 dscale = zero
520 dsum = one
521 DO 30 i = 1, n2
522 CALL dlassq( n1, s( n2+1, i ), 1, dscale, dsum )
523 30 CONTINUE
524 brqa21 = dscale*sqrt( dsum )
525*
526* Triangularize the B-part by a QR factorization.
527* Apply transformation (from right) to A-part, giving S.
528*
529 CALL dgeqr2( m, m, tcpy, ldst, taul, work, linfo )
530 IF( linfo.NE.0 )
531 \$ GO TO 70
532 CALL dorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
533 \$ work, info )
534 CALL dorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop, ldst,
535 \$ work, info )
536 IF( linfo.NE.0 )
537 \$ GO TO 70
538*
539* Compute F-norm(S21) in BQRA21. (T21 is 0.)
540*
541 dscale = zero
542 dsum = one
543 DO 40 i = 1, n2
544 CALL dlassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
545 40 CONTINUE
546 bqra21 = dscale*sqrt( dsum )
547*
548* Decide which method to use.
549* Weak stability test:
550* F-norm(S21) <= O(EPS * F-norm((S)))
551*
552 IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresha ) THEN
553 CALL dlacpy( 'F', m, m, scpy, ldst, s, ldst )
554 CALL dlacpy( 'F', m, m, tcpy, ldst, t, ldst )
555 CALL dlacpy( 'F', m, m, ircop, ldst, ir, ldst )
556 CALL dlacpy( 'F', m, m, licop, ldst, li, ldst )
557 ELSE IF( brqa21.GE.thresha ) THEN
558 GO TO 70
559 END IF
560*
561* Set lower triangle of B-part to zero
562*
563 CALL dlaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
564*
565 IF( wands ) THEN
566*
567* Strong stability test:
568* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
569* and
570* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
571*
572 CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
573 \$ m )
574 CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
575 \$ work, m )
576 CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
577 \$ work( m*m+1 ), m )
578 dscale = zero
579 dsum = one
580 CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
581 sa = dscale*sqrt( dsum )
582*
583 CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
584 \$ m )
585 CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
586 \$ work, m )
587 CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
588 \$ work( m*m+1 ), m )
589 dscale = zero
590 dsum = one
591 CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
592 sb = dscale*sqrt( dsum )
593 strong = sa.LE.thresha .AND. sb.LE.threshb
594 IF( .NOT.strong )
595 \$ GO TO 70
596*
597 END IF
598*
599* If the swap is accepted ("weakly" and "strongly"), apply the
600* transformations and set N1-by-N2 (2,1)-block to zero.
601*
602 CALL dlaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
603*
604* copy back M-by-M diagonal block starting at index J1 of (A, B)
605*
606 CALL dlacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
607 CALL dlacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
608 CALL dlaset( 'Full', ldst, ldst, zero, zero, t, ldst )
609*
610* Standardize existing 2-by-2 blocks.
611*
612 CALL dlaset( 'Full', m, m, zero, zero, work, m )
613 work( 1 ) = one
614 t( 1, 1 ) = one
615 idum = lwork - m*m - 2
616 IF( n2.GT.1 ) THEN
617 CALL dlagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
618 \$ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
619 work( m+1 ) = -work( 2 )
620 work( m+2 ) = work( 1 )
621 t( n2, n2 ) = t( 1, 1 )
622 t( 1, 2 ) = -t( 2, 1 )
623 END IF
624 work( m*m ) = one
625 t( m, m ) = one
626*
627 IF( n1.GT.1 ) THEN
628 CALL dlagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
629 \$ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
630 \$ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
631 \$ t( m, m-1 ) )
632 work( m*m ) = work( n2*m+n2+1 )
633 work( m*m-1 ) = -work( n2*m+n2+2 )
634 t( m, m ) = t( n2+1, n2+1 )
635 t( m-1, m ) = -t( m, m-1 )
636 END IF
637 CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
638 \$ lda, zero, work( m*m+1 ), n2 )
639 CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
640 \$ lda )
641 CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
642 \$ ldb, zero, work( m*m+1 ), n2 )
643 CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
644 \$ ldb )
645 CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
646 \$ work( m*m+1 ), m )
647 CALL dlacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
648 CALL dgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
649 \$ t( n2+1, n2+1 ), ldst, zero, work, n2 )
650 CALL dlacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
651 CALL dgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
652 \$ t( n2+1, n2+1 ), ldst, zero, work, n2 )
653 CALL dlacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
654 CALL dgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
655 \$ work, m )
656 CALL dlacpy( 'Full', m, m, work, m, ir, ldst )
657*
658* Accumulate transformations into Q and Z if requested.
659*
660 IF( wantq ) THEN
661 CALL dgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
662 \$ ldst, zero, work, n )
663 CALL dlacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
664*
665 END IF
666*
667 IF( wantz ) THEN
668 CALL dgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
669 \$ ldst, zero, work, n )
670 CALL dlacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
671*
672 END IF
673*
674* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
675* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
676*
677 i = j1 + m
678 IF( i.LE.n ) THEN
679 CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
680 \$ a( j1, i ), lda, zero, work, m )
681 CALL dlacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
682 CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
683 \$ b( j1, i ), ldb, zero, work, m )
684 CALL dlacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
685 END IF
686 i = j1 - 1
687 IF( i.GT.0 ) THEN
688 CALL dgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
689 \$ ldst, zero, work, i )
690 CALL dlacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
691 CALL dgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
692 \$ ldst, zero, work, i )
693 CALL dlacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
694 END IF
695*
696* Exit with INFO = 0 if swap was successfully performed.
697*
698 RETURN
699*
700 END IF
701*
702* Exit with INFO = 1 if swap was rejected.
703*
704 70 CONTINUE
705*
706 info = 1
707 RETURN
708*
709* End of DTGEX2
710*
711 END
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dgeqr2(m, n, a, lda, tau, work, info)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition dgeqr2.f:130
subroutine dgerq2(m, n, a, lda, tau, work, info)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.
Definition dgerq2.f:123
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 dlagv2(a, lda, b, ldb, alphar, alphai, beta, csl, snl, csr, snr)
DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,...
Definition dlagv2.f:157
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:124
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dtgex2(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, n1, n2, work, lwork, info)
DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition dtgex2.f:221
subroutine dtgsy2(trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, iwork, pq, info)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition dtgsy2.f:274
subroutine dorg2r(m, n, k, a, lda, tau, work, info)
DORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
Definition dorg2r.f:114
subroutine dorgr2(m, n, k, a, lda, tau, work, info)
DORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...
Definition dorgr2.f:114
subroutine dorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition dorm2r.f:159
subroutine dormr2(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
DORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
Definition dormr2.f:159