LAPACK 3.12.1
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 217 of file dtgex2.f.

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