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

◆ stgex2()

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

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

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

Purpose:
!>
!> STGEX2 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 SGGES), 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 REAL array, dimension (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 REAL array, dimension (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 REAL 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 REAL 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 REAL array, dimension (MAX(1,LWORK)).
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          LWORK >=  MAX( 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 stgex2.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 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
230 $ WORK( * ), Z( LDZ, * )
231* ..
232*
233* =====================================================================
234* Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
235* loops. Sven Hammarling, 1/5/02.
236*
237* .. Parameters ..
238 REAL ZERO, ONE
239 parameter( zero = 0.0e+0, one = 1.0e+0 )
240 REAL TWENTY
241 parameter( twenty = 2.0e+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 REAL BQRA21, BRQA21, DDUM, DNORMA, DNORMB,
251 $ DSCALE,
252 $ DSUM, EPS, F, G, SA, SB, SCALE, SMLNUM,
253 $ THRESHA, THRESHB
254* ..
255* .. Local Arrays ..
256 INTEGER IWORK( LDST + 2 )
257 REAL AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
258 $ IRCOP( LDST, LDST ), LI( LDST, LDST ),
259 $ LICOP( LDST, LDST ), S( LDST, LDST ),
260 $ SCPY( LDST, LDST ), T( LDST, LDST ),
261 $ TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
262* ..
263* .. External Functions ..
264 REAL SLAMCH
265 EXTERNAL slamch
266* ..
267* .. External Subroutines ..
268 EXTERNAL sgemm, sgeqr2, sgerq2, slacpy, slagv2,
269 $ slartg,
270 $ slaset, slassq, sorg2r, sorgr2, sorm2r, sormr2,
271 $ srot, sscal, stgsy2
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( n*m, m*m*2 ) ) THEN
288 info = -16
289 work( 1 ) = real( max( 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 slaset( 'Full', ldst, ldst, zero, zero, li, ldst )
299 CALL slaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
300 CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
301 CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
302*
303* Compute threshold for testing acceptance of swapping.
304*
305 eps = slamch( 'P' )
306 smlnum = slamch( 'S' ) / eps
307 dscale = zero
308 dsum = one
309 CALL slacpy( 'Full', m, m, s, ldst, work, m )
310 CALL slassq( m*m, work, 1, dscale, dsum )
311 dnorma = dscale*sqrt( dsum )
312 dscale = zero
313 dsum = one
314 CALL slacpy( 'Full', m, m, t, ldst, work, m )
315 CALL slassq( 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 slartg( 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 srot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
344 $ ir( 2, 1 ) )
345 CALL srot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
346 $ ir( 2, 1 ) )
347 IF( sa.GE.sb ) THEN
348 CALL slartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2,
349 $ 1 ),
350 $ ddum )
351 ELSE
352 CALL slartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2,
353 $ 1 ),
354 $ ddum )
355 END IF
356 CALL srot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
357 $ li( 2, 1 ) )
358 CALL srot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
359 $ li( 2, 1 ) )
360 li( 2, 2 ) = li( 1, 1 )
361 li( 1, 2 ) = -li( 2, 1 )
362*
363* Weak stability test: |S21| <= O(EPS F-norm((A)))
364* and |T21| <= O(EPS F-norm((B)))
365*
366 weak = abs( s( 2, 1 ) ) .LE. thresha .AND.
367 $ abs( t( 2, 1 ) ) .LE. threshb
368 IF( .NOT.weak )
369 $ GO TO 70
370*
371 IF( wands ) THEN
372*
373* Strong stability test:
374* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
375* and
376* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
377*
378 CALL slacpy( 'Full', m, m, a( j1, j1 ), lda,
379 $ work( m*m+1 ),
380 $ m )
381 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst,
382 $ zero,
383 $ work, m )
384 CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst,
385 $ one,
386 $ work( m*m+1 ), m )
387 dscale = zero
388 dsum = one
389 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
390 sa = dscale*sqrt( dsum )
391*
392 CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb,
393 $ work( m*m+1 ),
394 $ m )
395 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst,
396 $ zero,
397 $ work, m )
398 CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst,
399 $ one,
400 $ work( m*m+1 ), m )
401 dscale = zero
402 dsum = one
403 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
404 sb = dscale*sqrt( dsum )
405 strong = sa.LE.thresha .AND. sb.LE.threshb
406 IF( .NOT.strong )
407 $ GO TO 70
408 END IF
409*
410* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
411* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
412*
413 CALL srot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
414 $ ir( 2, 1 ) )
415 CALL srot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
416 $ ir( 2, 1 ) )
417 CALL srot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
418 $ li( 1, 1 ), li( 2, 1 ) )
419 CALL srot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
420 $ li( 1, 1 ), li( 2, 1 ) )
421*
422* Set N1-by-N2 (2,1) - blocks to ZERO.
423*
424 a( j1+1, j1 ) = zero
425 b( j1+1, j1 ) = zero
426*
427* Accumulate transformations into Q and Z if requested.
428*
429 IF( wantz )
430 $ CALL srot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
431 $ ir( 2, 1 ) )
432 IF( wantq )
433 $ CALL srot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
434 $ li( 2, 1 ) )
435*
436* Exit with INFO = 0 if swap was successfully performed.
437*
438 RETURN
439*
440 ELSE
441*
442* CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
443* and 2-by-2 blocks.
444*
445* Solve the generalized Sylvester equation
446* S11 * R - L * S22 = SCALE * S12
447* T11 * R - L * T22 = SCALE * T12
448* for R and L. Solutions in LI and IR.
449*
450 CALL slacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
451 CALL slacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
452 $ ir( n2+1, n1+1 ), ldst )
453 CALL stgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
454 $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
455 $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
456 $ linfo )
457 IF( linfo.NE.0 )
458 $ GO TO 70
459*
460* Compute orthogonal matrix QL:
461*
462* QL**T * LI = [ TL ]
463* [ 0 ]
464* where
465* LI = [ -L ]
466* [ SCALE * identity(N2) ]
467*
468 DO 10 i = 1, n2
469 CALL sscal( n1, -one, li( 1, i ), 1 )
470 li( n1+i, i ) = scale
471 10 CONTINUE
472 CALL sgeqr2( m, n2, li, ldst, taul, work, linfo )
473 IF( linfo.NE.0 )
474 $ GO TO 70
475 CALL sorg2r( m, m, n2, li, ldst, taul, work, linfo )
476 IF( linfo.NE.0 )
477 $ GO TO 70
478*
479* Compute orthogonal matrix RQ:
480*
481* IR * RQ**T = [ 0 TR],
482*
483* where IR = [ SCALE * identity(N1), R ]
484*
485 DO 20 i = 1, n1
486 ir( n2+i, i ) = scale
487 20 CONTINUE
488 CALL sgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
489 IF( linfo.NE.0 )
490 $ GO TO 70
491 CALL sorgr2( m, m, n1, ir, ldst, taur, work, linfo )
492 IF( linfo.NE.0 )
493 $ GO TO 70
494*
495* Perform the swapping tentatively:
496*
497 CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
498 $ work, m )
499 CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero,
500 $ s,
501 $ ldst )
502 CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
503 $ work, m )
504 CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero,
505 $ t,
506 $ ldst )
507 CALL slacpy( 'F', m, m, s, ldst, scpy, ldst )
508 CALL slacpy( 'F', m, m, t, ldst, tcpy, ldst )
509 CALL slacpy( 'F', m, m, ir, ldst, ircop, ldst )
510 CALL slacpy( 'F', m, m, li, ldst, licop, ldst )
511*
512* Triangularize the B-part by an RQ factorization.
513* Apply transformation (from left) to A-part, giving S.
514*
515 CALL sgerq2( m, m, t, ldst, taur, work, linfo )
516 IF( linfo.NE.0 )
517 $ GO TO 70
518 CALL sormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst,
519 $ work,
520 $ linfo )
521 IF( linfo.NE.0 )
522 $ GO TO 70
523 CALL sormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst,
524 $ work,
525 $ linfo )
526 IF( linfo.NE.0 )
527 $ GO TO 70
528*
529* Compute F-norm(S21) in BRQA21. (T21 is 0.)
530*
531 dscale = zero
532 dsum = one
533 DO 30 i = 1, n2
534 CALL slassq( n1, s( n2+1, i ), 1, dscale, dsum )
535 30 CONTINUE
536 brqa21 = dscale*sqrt( dsum )
537*
538* Triangularize the B-part by a QR factorization.
539* Apply transformation (from right) to A-part, giving S.
540*
541 CALL sgeqr2( m, m, tcpy, ldst, taul, work, linfo )
542 IF( linfo.NE.0 )
543 $ GO TO 70
544 CALL sorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy,
545 $ ldst,
546 $ work, info )
547 CALL sorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop,
548 $ ldst,
549 $ work, info )
550 IF( linfo.NE.0 )
551 $ GO TO 70
552*
553* Compute F-norm(S21) in BQRA21. (T21 is 0.)
554*
555 dscale = zero
556 dsum = one
557 DO 40 i = 1, n2
558 CALL slassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
559 40 CONTINUE
560 bqra21 = dscale*sqrt( dsum )
561*
562* Decide which method to use.
563* Weak stability test:
564* F-norm(S21) <= O(EPS * F-norm((S)))
565*
566 IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresha ) THEN
567 CALL slacpy( 'F', m, m, scpy, ldst, s, ldst )
568 CALL slacpy( 'F', m, m, tcpy, ldst, t, ldst )
569 CALL slacpy( 'F', m, m, ircop, ldst, ir, ldst )
570 CALL slacpy( 'F', m, m, licop, ldst, li, ldst )
571 ELSE IF( brqa21.GE.thresha ) THEN
572 GO TO 70
573 END IF
574*
575* Set lower triangle of B-part to zero
576*
577 CALL slaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
578*
579 IF( wands ) THEN
580*
581* Strong stability test:
582* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
583* and
584* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
585*
586 CALL slacpy( 'Full', m, m, a( j1, j1 ), lda,
587 $ work( m*m+1 ),
588 $ m )
589 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst,
590 $ zero,
591 $ work, m )
592 CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst,
593 $ one,
594 $ work( m*m+1 ), m )
595 dscale = zero
596 dsum = one
597 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
598 sa = dscale*sqrt( dsum )
599*
600 CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb,
601 $ work( m*m+1 ),
602 $ m )
603 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst,
604 $ zero,
605 $ work, m )
606 CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst,
607 $ one,
608 $ work( m*m+1 ), m )
609 dscale = zero
610 dsum = one
611 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
612 sb = dscale*sqrt( dsum )
613 strong = sa.LE.thresha .AND. sb.LE.threshb
614 IF( .NOT.strong )
615 $ GO TO 70
616*
617 END IF
618*
619* If the swap is accepted ("weakly" and "strongly"), apply the
620* transformations and set N1-by-N2 (2,1)-block to zero.
621*
622 CALL slaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
623*
624* copy back M-by-M diagonal block starting at index J1 of (A, B)
625*
626 CALL slacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
627 CALL slacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
628 CALL slaset( 'Full', ldst, ldst, zero, zero, t, ldst )
629*
630* Standardize existing 2-by-2 blocks.
631*
632 CALL slaset( 'Full', m, m, zero, zero, work, m )
633 work( 1 ) = one
634 t( 1, 1 ) = one
635 idum = lwork - m*m - 2
636 IF( n2.GT.1 ) THEN
637 CALL slagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai,
638 $ be,
639 $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
640 work( m+1 ) = -work( 2 )
641 work( m+2 ) = work( 1 )
642 t( n2, n2 ) = t( 1, 1 )
643 t( 1, 2 ) = -t( 2, 1 )
644 END IF
645 work( m*m ) = one
646 t( m, m ) = one
647*
648 IF( n1.GT.1 ) THEN
649 CALL slagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ),
650 $ ldb,
651 $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
652 $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
653 $ t( m, m-1 ) )
654 work( m*m ) = work( n2*m+n2+1 )
655 work( m*m-1 ) = -work( n2*m+n2+2 )
656 t( m, m ) = t( n2+1, n2+1 )
657 t( m-1, m ) = -t( m, m-1 )
658 END IF
659 CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1,
660 $ j1+n2 ),
661 $ lda, zero, work( m*m+1 ), n2 )
662 CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1,
663 $ j1+n2 ),
664 $ lda )
665 CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1,
666 $ j1+n2 ),
667 $ ldb, zero, work( m*m+1 ), n2 )
668 CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1,
669 $ j1+n2 ),
670 $ ldb )
671 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
672 $ work( m*m+1 ), m )
673 CALL slacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
674 CALL sgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
675 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
676 CALL slacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
677 CALL sgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
678 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
679 CALL slacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
680 CALL sgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
681 $ work, m )
682 CALL slacpy( 'Full', m, m, work, m, ir, ldst )
683*
684* Accumulate transformations into Q and Z if requested.
685*
686 IF( wantq ) THEN
687 CALL sgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
688 $ ldst, zero, work, n )
689 CALL slacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
690*
691 END IF
692*
693 IF( wantz ) THEN
694 CALL sgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
695 $ ldst, zero, work, n )
696 CALL slacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
697*
698 END IF
699*
700* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
701* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
702*
703 i = j1 + m
704 IF( i.LE.n ) THEN
705 CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
706 $ a( j1, i ), lda, zero, work, m )
707 CALL slacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
708 CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
709 $ b( j1, i ), ldb, zero, work, m )
710 CALL slacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
711 END IF
712 i = j1 - 1
713 IF( i.GT.0 ) THEN
714 CALL sgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
715 $ ldst, zero, work, i )
716 CALL slacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
717 CALL sgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
718 $ ldst, zero, work, i )
719 CALL slacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
720 END IF
721*
722* Exit with INFO = 0 if swap was successfully performed.
723*
724 RETURN
725*
726 END IF
727*
728* Exit with INFO = 1 if swap was rejected.
729*
730 70 CONTINUE
731*
732 info = 1
733 RETURN
734*
735* End of STGEX2
736*
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
sgeqr2
subroutine sgeqr2(m, n, a, lda, tau, work, info)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
Definition sgeqr2.f:128
sgerq2
subroutine sgerq2(m, n, a, lda, tau, work, info)
SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm.
Definition sgerq2.f:121
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slagv2
subroutine slagv2(a, lda, b, ldb, alphar, alphai, beta, csl, snl, csr, snr)
SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,...
Definition slagv2.f:156
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
slassq
subroutine slassq(n, x, incx, scale, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:122
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
stgsy2
subroutine stgsy2(trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, iwork, pq, info)
STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition stgsy2.f:273
sorg2r
subroutine sorg2r(m, n, k, a, lda, tau, work, info)
SORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
Definition sorg2r.f:112
sorgr2
subroutine sorgr2(m, n, k, a, lda, tau, work, info)
SORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...
Definition sorgr2.f:112
sorm2r
subroutine sorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition sorm2r.f:157
sormr2
subroutine sormr2(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
SORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
Definition sormr2.f:157
Here is the call graph for this function:
Here is the caller graph for this function: