LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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 )
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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f90:137
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 dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
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 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 dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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 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 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 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 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 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 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
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
Here is the call graph for this function:
Here is the caller graph for this function: