LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2015
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 223 of file dtgex2.f.

223 *
224 * -- LAPACK auxiliary routine (version 3.6.0) --
225 * -- LAPACK is a software package provided by Univ. of Tennessee, --
226 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227 * November 2015
228 *
229 * .. Scalar Arguments ..
230  LOGICAL wantq, wantz
231  INTEGER info, j1, lda, ldb, ldq, ldz, lwork, n, n1, n2
232 * ..
233 * .. Array Arguments ..
234  DOUBLE PRECISION a( lda, * ), b( ldb, * ), q( ldq, * ),
235  $ work( * ), z( ldz, * )
236 * ..
237 *
238 * =====================================================================
239 * Replaced various illegal calls to DCOPY by calls to DLASET, or by DO
240 * loops. Sven Hammarling, 1/5/02.
241 *
242 * .. Parameters ..
243  DOUBLE PRECISION zero, one
244  parameter ( zero = 0.0d+0, one = 1.0d+0 )
245  DOUBLE PRECISION twenty
246  parameter ( twenty = 2.0d+01 )
247  INTEGER ldst
248  parameter ( ldst = 4 )
249  LOGICAL wands
250  parameter ( wands = .true. )
251 * ..
252 * .. Local Scalars ..
253  LOGICAL dtrong, weak
254  INTEGER i, idum, linfo, m
255  DOUBLE PRECISION bqra21, brqa21, ddum, dnorm, dscale, dsum, eps,
256  $ f, g, sa, sb, scale, smlnum, ss, thresh, ws
257 * ..
258 * .. Local Arrays ..
259  INTEGER iwork( ldst )
260  DOUBLE PRECISION ai( 2 ), ar( 2 ), be( 2 ), ir( ldst, ldst ),
261  $ ircop( ldst, ldst ), li( ldst, ldst ),
262  $ licop( ldst, ldst ), s( ldst, ldst ),
263  $ scpy( ldst, ldst ), t( ldst, ldst ),
264  $ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
265 * ..
266 * .. External Functions ..
267  DOUBLE PRECISION dlamch
268  EXTERNAL dlamch
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL dgemm, dgeqr2, dgerq2, dlacpy, dlagv2, dlartg,
272  $ dlaset, dlassq, dorg2r, dorgr2, dorm2r, dormr2,
273  $ drot, dscal, dtgsy2
274 * ..
275 * .. Intrinsic Functions ..
276  INTRINSIC abs, max, sqrt
277 * ..
278 * .. Executable Statements ..
279 *
280  info = 0
281 *
282 * Quick return if possible
283 *
284  IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
285  $ RETURN
286  IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
287  $ RETURN
288  m = n1 + n2
289  IF( lwork.LT.max( 1, n*m, m*m*2 ) ) THEN
290  info = -16
291  work( 1 ) = max( 1, n*m, m*m*2 )
292  RETURN
293  END IF
294 *
295  weak = .false.
296  dtrong = .false.
297 *
298 * Make a local copy of selected block
299 *
300  CALL dlaset( 'Full', ldst, ldst, zero, zero, li, ldst )
301  CALL dlaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
302  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
303  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
304 *
305 * Compute threshold for testing acceptance of swapping.
306 *
307  eps = dlamch( 'P' )
308  smlnum = dlamch( 'S' ) / eps
309  dscale = zero
310  dsum = one
311  CALL dlacpy( 'Full', m, m, s, ldst, work, m )
312  CALL dlassq( m*m, work, 1, dscale, dsum )
313  CALL dlacpy( 'Full', m, m, t, ldst, work, m )
314  CALL dlassq( m*m, work, 1, dscale, dsum )
315  dnorm = dscale*sqrt( dsum )
316 *
317 * THRES has been changed from
318 * THRESH = MAX( TEN*EPS*SA, SMLNUM )
319 * to
320 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
321 * on 04/01/10.
322 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
323 * Jim Demmel and Guillaume Revy. See forum post 1783.
324 *
325  thresh = max( twenty*eps*dnorm, smlnum )
326 *
327  IF( m.EQ.2 ) THEN
328 *
329 * CASE 1: Swap 1-by-1 and 1-by-1 blocks.
330 *
331 * Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
332 * using Givens rotations and perform the swap tentatively.
333 *
334  f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
335  g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
336  sb = abs( t( 2, 2 ) )
337  sa = abs( s( 2, 2 ) )
338  CALL dlartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
339  ir( 2, 1 ) = -ir( 1, 2 )
340  ir( 2, 2 ) = ir( 1, 1 )
341  CALL drot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
342  $ ir( 2, 1 ) )
343  CALL drot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
344  $ ir( 2, 1 ) )
345  IF( sa.GE.sb ) THEN
346  CALL dlartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
347  $ ddum )
348  ELSE
349  CALL dlartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350  $ ddum )
351  END IF
352  CALL drot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
353  $ li( 2, 1 ) )
354  CALL drot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
355  $ li( 2, 1 ) )
356  li( 2, 2 ) = li( 1, 1 )
357  li( 1, 2 ) = -li( 2, 1 )
358 *
359 * Weak stability test:
360 * |S21| + |T21| <= O(EPS * F-norm((S, T)))
361 *
362  ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
363  weak = ws.LE.thresh
364  IF( .NOT.weak )
365  $ GO TO 70
366 *
367  IF( wands ) THEN
368 *
369 * Strong stability test:
370 * F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,B)))
371 *
372  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
373  $ m )
374  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
375  $ work, m )
376  CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
377  $ work( m*m+1 ), m )
378  dscale = zero
379  dsum = one
380  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
381 *
382  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
383  $ m )
384  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
385  $ work, m )
386  CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
387  $ work( m*m+1 ), m )
388  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389  ss = dscale*sqrt( dsum )
390  dtrong = ss.LE.thresh
391  IF( .NOT.dtrong )
392  $ GO TO 70
393  END IF
394 *
395 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
396 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
397 *
398  CALL drot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
399  $ ir( 2, 1 ) )
400  CALL drot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
401  $ ir( 2, 1 ) )
402  CALL drot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403  $ li( 1, 1 ), li( 2, 1 ) )
404  CALL drot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
405  $ li( 1, 1 ), li( 2, 1 ) )
406 *
407 * Set N1-by-N2 (2,1) - blocks to ZERO.
408 *
409  a( j1+1, j1 ) = zero
410  b( j1+1, j1 ) = zero
411 *
412 * Accumulate transformations into Q and Z if requested.
413 *
414  IF( wantz )
415  $ CALL drot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
416  $ ir( 2, 1 ) )
417  IF( wantq )
418  $ CALL drot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
419  $ li( 2, 1 ) )
420 *
421 * Exit with INFO = 0 if swap was successfully performed.
422 *
423  RETURN
424 *
425  ELSE
426 *
427 * CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
428 * and 2-by-2 blocks.
429 *
430 * Solve the generalized Sylvester equation
431 * S11 * R - L * S22 = SCALE * S12
432 * T11 * R - L * T22 = SCALE * T12
433 * for R and L. Solutions in LI and IR.
434 *
435  CALL dlacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436  CALL dlacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
437  $ ir( n2+1, n1+1 ), ldst )
438  CALL dtgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
439  $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
440  $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
441  $ linfo )
442 *
443 * Compute orthogonal matrix QL:
444 *
445 * QL**T * LI = [ TL ]
446 * [ 0 ]
447 * where
448 * LI = [ -L ]
449 * [ SCALE * identity(N2) ]
450 *
451  DO 10 i = 1, n2
452  CALL dscal( n1, -one, li( 1, i ), 1 )
453  li( n1+i, i ) = scale
454  10 CONTINUE
455  CALL dgeqr2( m, n2, li, ldst, taul, work, linfo )
456  IF( linfo.NE.0 )
457  $ GO TO 70
458  CALL dorg2r( m, m, n2, li, ldst, taul, work, linfo )
459  IF( linfo.NE.0 )
460  $ GO TO 70
461 *
462 * Compute orthogonal matrix RQ:
463 *
464 * IR * RQ**T = [ 0 TR],
465 *
466 * where IR = [ SCALE * identity(N1), R ]
467 *
468  DO 20 i = 1, n1
469  ir( n2+i, i ) = scale
470  20 CONTINUE
471  CALL dgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
472  IF( linfo.NE.0 )
473  $ GO TO 70
474  CALL dorgr2( m, m, n1, ir, ldst, taur, work, linfo )
475  IF( linfo.NE.0 )
476  $ GO TO 70
477 *
478 * Perform the swapping tentatively:
479 *
480  CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
481  $ work, m )
482  CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
483  $ ldst )
484  CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
485  $ work, m )
486  CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
487  $ ldst )
488  CALL dlacpy( 'F', m, m, s, ldst, scpy, ldst )
489  CALL dlacpy( 'F', m, m, t, ldst, tcpy, ldst )
490  CALL dlacpy( 'F', m, m, ir, ldst, ircop, ldst )
491  CALL dlacpy( 'F', m, m, li, ldst, licop, ldst )
492 *
493 * Triangularize the B-part by an RQ factorization.
494 * Apply transformation (from left) to A-part, giving S.
495 *
496  CALL dgerq2( m, m, t, ldst, taur, work, linfo )
497  IF( linfo.NE.0 )
498  $ GO TO 70
499  CALL dormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
500  $ linfo )
501  IF( linfo.NE.0 )
502  $ GO TO 70
503  CALL dormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst, work,
504  $ linfo )
505  IF( linfo.NE.0 )
506  $ GO TO 70
507 *
508 * Compute F-norm(S21) in BRQA21. (T21 is 0.)
509 *
510  dscale = zero
511  dsum = one
512  DO 30 i = 1, n2
513  CALL dlassq( n1, s( n2+1, i ), 1, dscale, dsum )
514  30 CONTINUE
515  brqa21 = dscale*sqrt( dsum )
516 *
517 * Triangularize the B-part by a QR factorization.
518 * Apply transformation (from right) to A-part, giving S.
519 *
520  CALL dgeqr2( m, m, tcpy, ldst, taul, work, linfo )
521  IF( linfo.NE.0 )
522  $ GO TO 70
523  CALL dorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
524  $ work, info )
525  CALL dorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop, ldst,
526  $ work, info )
527  IF( linfo.NE.0 )
528  $ GO TO 70
529 *
530 * Compute F-norm(S21) in BQRA21. (T21 is 0.)
531 *
532  dscale = zero
533  dsum = one
534  DO 40 i = 1, n2
535  CALL dlassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
536  40 CONTINUE
537  bqra21 = dscale*sqrt( dsum )
538 *
539 * Decide which method to use.
540 * Weak stability test:
541 * F-norm(S21) <= O(EPS * F-norm((S, T)))
542 *
543  IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresh ) THEN
544  CALL dlacpy( 'F', m, m, scpy, ldst, s, ldst )
545  CALL dlacpy( 'F', m, m, tcpy, ldst, t, ldst )
546  CALL dlacpy( 'F', m, m, ircop, ldst, ir, ldst )
547  CALL dlacpy( 'F', m, m, licop, ldst, li, ldst )
548  ELSE IF( brqa21.GE.thresh ) THEN
549  GO TO 70
550  END IF
551 *
552 * Set lower triangle of B-part to zero
553 *
554  CALL dlaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
555 *
556  IF( wands ) THEN
557 *
558 * Strong stability test:
559 * F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))
560 *
561  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
562  $ m )
563  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
564  $ work, m )
565  CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
566  $ work( m*m+1 ), m )
567  dscale = zero
568  dsum = one
569  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
570 *
571  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
572  $ m )
573  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
574  $ work, m )
575  CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
576  $ work( m*m+1 ), m )
577  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578  ss = dscale*sqrt( dsum )
579  dtrong = ( ss.LE.thresh )
580  IF( .NOT.dtrong )
581  $ GO TO 70
582 *
583  END IF
584 *
585 * If the swap is accepted ("weakly" and "strongly"), apply the
586 * transformations and set N1-by-N2 (2,1)-block to zero.
587 *
588  CALL dlaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
589 *
590 * copy back M-by-M diagonal block starting at index J1 of (A, B)
591 *
592  CALL dlacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
593  CALL dlacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
594  CALL dlaset( 'Full', ldst, ldst, zero, zero, t, ldst )
595 *
596 * Standardize existing 2-by-2 blocks.
597 *
598  CALL dlaset( 'Full', m, m, zero, zero, work, m )
599  work( 1 ) = one
600  t( 1, 1 ) = one
601  idum = lwork - m*m - 2
602  IF( n2.GT.1 ) THEN
603  CALL dlagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
604  $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
605  work( m+1 ) = -work( 2 )
606  work( m+2 ) = work( 1 )
607  t( n2, n2 ) = t( 1, 1 )
608  t( 1, 2 ) = -t( 2, 1 )
609  END IF
610  work( m*m ) = one
611  t( m, m ) = one
612 *
613  IF( n1.GT.1 ) THEN
614  CALL dlagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
615  $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
616  $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
617  $ t( m, m-1 ) )
618  work( m*m ) = work( n2*m+n2+1 )
619  work( m*m-1 ) = -work( n2*m+n2+2 )
620  t( m, m ) = t( n2+1, n2+1 )
621  t( m-1, m ) = -t( m, m-1 )
622  END IF
623  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
624  $ lda, zero, work( m*m+1 ), n2 )
625  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
626  $ lda )
627  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
628  $ ldb, zero, work( m*m+1 ), n2 )
629  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
630  $ ldb )
631  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
632  $ work( m*m+1 ), m )
633  CALL dlacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
634  CALL dgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
635  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
636  CALL dlacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
637  CALL dgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
638  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
639  CALL dlacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
640  CALL dgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
641  $ work, m )
642  CALL dlacpy( 'Full', m, m, work, m, ir, ldst )
643 *
644 * Accumulate transformations into Q and Z if requested.
645 *
646  IF( wantq ) THEN
647  CALL dgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
648  $ ldst, zero, work, n )
649  CALL dlacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
650 *
651  END IF
652 *
653  IF( wantz ) THEN
654  CALL dgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
655  $ ldst, zero, work, n )
656  CALL dlacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
657 *
658  END IF
659 *
660 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
661 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
662 *
663  i = j1 + m
664  IF( i.LE.n ) THEN
665  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
666  $ a( j1, i ), lda, zero, work, m )
667  CALL dlacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
668  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
669  $ b( j1, i ), ldb, zero, work, m )
670  CALL dlacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
671  END IF
672  i = j1 - 1
673  IF( i.GT.0 ) THEN
674  CALL dgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
675  $ ldst, zero, work, i )
676  CALL dlacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
677  CALL dgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
678  $ ldst, zero, work, i )
679  CALL dlacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
680  END IF
681 *
682 * Exit with INFO = 0 if swap was successfully performed.
683 *
684  RETURN
685 *
686  END IF
687 *
688 * Exit with INFO = 1 if swap was rejected.
689 *
690  70 CONTINUE
691 *
692  info = 1
693  RETURN
694 *
695 * End of DTGEX2
696 *
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:112
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:125
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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: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:105
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
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:161
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
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:161
dlassq
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
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:116
dlartg
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
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:159
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:116
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:276

Here is the call graph for this function:

Here is the caller graph for this function: