LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dtgex2()

subroutine dtgex2 ( logical  wantq,
logical  wantz,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision, dimension( ldq, * )  q,
integer  ldq,
double precision, dimension( ldz, * )  z,
integer  ldz,
integer  j1,
integer  n1,
integer  n2,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation.

Download DTGEX2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
 of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
 (A, B) by an orthogonal equivalence transformation.

 (A, B) must be in generalized real Schur canonical form (as returned
 by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
 diagonal blocks. B is upper triangular.

 Optionally, the matrices Q and Z of generalized Schur vectors are
 updated.

        Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
        Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
Parameters
[in]WANTQ
          WANTQ is LOGICAL
          .TRUE. : update the left transformation matrix Q;
          .FALSE.: do not update Q.
[in]WANTZ
          WANTZ is LOGICAL
          .TRUE. : update the right transformation matrix Z;
          .FALSE.: do not update Z.
[in]N
          N is INTEGER
          The order of the matrices A and B. N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimensions (LDA,N)
          On entry, the matrix A in the pair (A, B).
          On exit, the updated matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,N).
[in,out]B
          B is DOUBLE PRECISION array, dimensions (LDB,N)
          On entry, the matrix B in the pair (A, B).
          On exit, the updated matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[in,out]Q
          Q is DOUBLE PRECISION array, dimension (LDQ,N)
          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
          On exit, the updated matrix Q.
          Not referenced if WANTQ = .FALSE..
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q. LDQ >= 1.
          If WANTQ = .TRUE., LDQ >= N.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ,N)
          On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
          On exit, the updated matrix Z.
          Not referenced if WANTZ = .FALSE..
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z. LDZ >= 1.
          If WANTZ = .TRUE., LDZ >= N.
[in]J1
          J1 is INTEGER
          The index to the first block (A11, B11). 1 <= J1 <= N.
[in]N1
          N1 is INTEGER
          The order of the first block (A11, B11). N1 = 0, 1 or 2.
[in]N2
          N2 is INTEGER
          The order of the second block (A22, B22). N2 = 0, 1 or 2.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)).
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >=  MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 )
[out]INFO
          INFO is INTEGER
            =0: Successful exit
            >0: If INFO = 1, the transformed matrix (A, B) would be
                too far from generalized Schur form; the blocks are
                not swapped and (A, B) and (Q, Z) are unchanged.
                The problem of swapping is too ill-conditioned.
            <0: If INFO = -16: LWORK is too small. Appropriate value
                for LWORK is returned in WORK(1).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
In the current code both weak and strong stability tests are performed. The user can omit the strong stability test by changing the internal logical parameter WANDS to .FALSE.. See ref. [2] for details.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.

  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
      Estimation: Theory, Algorithms and Software,
      Report UMINF - 94.04, Department of Computing Science, Umea
      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
      Note 87. To appear in Numerical Algorithms, 1996.

Definition at line 219 of file dtgex2.f.

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,
270 $ dlaset, dlassq, dorg2r, dorgr2, dorm2r, dormr2,
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*
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgeqr2
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
dgerq2
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
dlacpy
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
dlagv2
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlaset
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
dlassq
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:124
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtgsy2
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
dorg2r
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
dorgr2
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
dorm2r
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
dormr2
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
Here is the call graph for this function:
Here is the caller graph for this function: