LAPACK 3.12.0
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 219 of file stgex2.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 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
232 $ WORK( * ), Z( LDZ, * )
233* ..
234*
235* =====================================================================
236* Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
237* loops. Sven Hammarling, 1/5/02.
238*
239* .. Parameters ..
240 REAL ZERO, ONE
241 parameter( zero = 0.0e+0, one = 1.0e+0 )
242 REAL TWENTY
243 parameter( twenty = 2.0e+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 REAL BQRA21, BRQA21, DDUM, DNORMA, DNORMB,
253 $ DSCALE,
254 $ DSUM, EPS, F, G, SA, SB, SCALE, SMLNUM,
255 $ THRESHA, THRESHB
256* ..
257* .. Local Arrays ..
258 INTEGER IWORK( LDST + 2 )
259 REAL AI( 2 ), AR( 2 ), BE( 2 ), IR( LDST, LDST ),
260 $ IRCOP( LDST, LDST ), LI( LDST, LDST ),
261 $ LICOP( LDST, LDST ), S( LDST, LDST ),
262 $ SCPY( LDST, LDST ), T( LDST, LDST ),
263 $ TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
264* ..
265* .. External Functions ..
266 REAL SLAMCH
267 EXTERNAL slamch
268* ..
269* .. External Subroutines ..
270 EXTERNAL sgemm, sgeqr2, sgerq2, slacpy, slagv2, slartg,
272 $ srot, sscal, stgsy2
273* ..
274* .. Intrinsic Functions ..
275 INTRINSIC abs, max, sqrt
276* ..
277* .. Executable Statements ..
278*
279 info = 0
280*
281* Quick return if possible
282*
283 IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
284 $ RETURN
285 IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
286 $ RETURN
287 m = n1 + n2
288 IF( lwork.LT.max( n*m, m*m*2 ) ) THEN
289 info = -16
290 work( 1 ) = max( n*m, m*m*2 )
291 RETURN
292 END IF
293*
294 weak = .false.
295 strong = .false.
296*
297* Make a local copy of selected block
298*
299 CALL slaset( 'Full', ldst, ldst, zero, zero, li, ldst )
300 CALL slaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
301 CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
302 CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
303*
304* Compute threshold for testing acceptance of swapping.
305*
306 eps = slamch( 'P' )
307 smlnum = slamch( 'S' ) / eps
308 dscale = zero
309 dsum = one
310 CALL slacpy( 'Full', m, m, s, ldst, work, m )
311 CALL slassq( m*m, work, 1, dscale, dsum )
312 dnorma = dscale*sqrt( dsum )
313 dscale = zero
314 dsum = one
315 CALL slacpy( 'Full', m, m, t, ldst, work, m )
316 CALL slassq( m*m, work, 1, dscale, dsum )
317 dnormb = dscale*sqrt( dsum )
318*
319* THRES has been changed from
320* THRESH = MAX( TEN*EPS*SA, SMLNUM )
321* to
322* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
323* on 04/01/10.
324* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
325* Jim Demmel and Guillaume Revy. See forum post 1783.
326*
327 thresha = max( twenty*eps*dnorma, smlnum )
328 threshb = max( twenty*eps*dnormb, smlnum )
329*
330 IF( m.EQ.2 ) THEN
331*
332* CASE 1: Swap 1-by-1 and 1-by-1 blocks.
333*
334* Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
335* using Givens rotations and perform the swap tentatively.
336*
337 f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
338 g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
339 sa = abs( s( 2, 2 ) ) * abs( t( 1, 1 ) )
340 sb = abs( s( 1, 1 ) ) * abs( t( 2, 2 ) )
341 CALL slartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
342 ir( 2, 1 ) = -ir( 1, 2 )
343 ir( 2, 2 ) = ir( 1, 1 )
344 CALL srot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
345 $ ir( 2, 1 ) )
346 CALL srot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
347 $ ir( 2, 1 ) )
348 IF( sa.GE.sb ) THEN
349 CALL slartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350 $ ddum )
351 ELSE
352 CALL slartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
353 $ ddum )
354 END IF
355 CALL srot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
356 $ li( 2, 1 ) )
357 CALL srot( 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 slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
378 $ m )
379 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
380 $ work, m )
381 CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
382 $ work( m*m+1 ), m )
383 dscale = zero
384 dsum = one
385 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
386 sa = dscale*sqrt( dsum )
387*
388 CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
389 $ m )
390 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
391 $ work, m )
392 CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
393 $ work( m*m+1 ), m )
394 dscale = zero
395 dsum = one
396 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
397 sb = dscale*sqrt( dsum )
398 strong = sa.LE.thresha .AND. sb.LE.threshb
399 IF( .NOT.strong )
400 $ GO TO 70
401 END IF
402*
403* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
404* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
405*
406 CALL srot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
407 $ ir( 2, 1 ) )
408 CALL srot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
409 $ ir( 2, 1 ) )
410 CALL srot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
411 $ li( 1, 1 ), li( 2, 1 ) )
412 CALL srot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
413 $ li( 1, 1 ), li( 2, 1 ) )
414*
415* Set N1-by-N2 (2,1) - blocks to ZERO.
416*
417 a( j1+1, j1 ) = zero
418 b( j1+1, j1 ) = zero
419*
420* Accumulate transformations into Q and Z if requested.
421*
422 IF( wantz )
423 $ CALL srot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
424 $ ir( 2, 1 ) )
425 IF( wantq )
426 $ CALL srot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
427 $ li( 2, 1 ) )
428*
429* Exit with INFO = 0 if swap was successfully performed.
430*
431 RETURN
432*
433 ELSE
434*
435* CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
436* and 2-by-2 blocks.
437*
438* Solve the generalized Sylvester equation
439* S11 * R - L * S22 = SCALE * S12
440* T11 * R - L * T22 = SCALE * T12
441* for R and L. Solutions in LI and IR.
442*
443 CALL slacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
444 CALL slacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
445 $ ir( n2+1, n1+1 ), ldst )
446 CALL stgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
447 $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
448 $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
449 $ linfo )
450 IF( linfo.NE.0 )
451 $ GO TO 70
452*
453* Compute orthogonal matrix QL:
454*
455* QL**T * LI = [ TL ]
456* [ 0 ]
457* where
458* LI = [ -L ]
459* [ SCALE * identity(N2) ]
460*
461 DO 10 i = 1, n2
462 CALL sscal( n1, -one, li( 1, i ), 1 )
463 li( n1+i, i ) = scale
464 10 CONTINUE
465 CALL sgeqr2( m, n2, li, ldst, taul, work, linfo )
466 IF( linfo.NE.0 )
467 $ GO TO 70
468 CALL sorg2r( m, m, n2, li, ldst, taul, work, linfo )
469 IF( linfo.NE.0 )
470 $ GO TO 70
471*
472* Compute orthogonal matrix RQ:
473*
474* IR * RQ**T = [ 0 TR],
475*
476* where IR = [ SCALE * identity(N1), R ]
477*
478 DO 20 i = 1, n1
479 ir( n2+i, i ) = scale
480 20 CONTINUE
481 CALL sgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
482 IF( linfo.NE.0 )
483 $ GO TO 70
484 CALL sorgr2( m, m, n1, ir, ldst, taur, work, linfo )
485 IF( linfo.NE.0 )
486 $ GO TO 70
487*
488* Perform the swapping tentatively:
489*
490 CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
491 $ work, m )
492 CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
493 $ ldst )
494 CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
495 $ work, m )
496 CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
497 $ ldst )
498 CALL slacpy( 'F', m, m, s, ldst, scpy, ldst )
499 CALL slacpy( 'F', m, m, t, ldst, tcpy, ldst )
500 CALL slacpy( 'F', m, m, ir, ldst, ircop, ldst )
501 CALL slacpy( 'F', m, m, li, ldst, licop, ldst )
502*
503* Triangularize the B-part by an RQ factorization.
504* Apply transformation (from left) to A-part, giving S.
505*
506 CALL sgerq2( m, m, t, ldst, taur, work, linfo )
507 IF( linfo.NE.0 )
508 $ GO TO 70
509 CALL sormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
510 $ linfo )
511 IF( linfo.NE.0 )
512 $ GO TO 70
513 CALL sormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst, work,
514 $ linfo )
515 IF( linfo.NE.0 )
516 $ GO TO 70
517*
518* Compute F-norm(S21) in BRQA21. (T21 is 0.)
519*
520 dscale = zero
521 dsum = one
522 DO 30 i = 1, n2
523 CALL slassq( n1, s( n2+1, i ), 1, dscale, dsum )
524 30 CONTINUE
525 brqa21 = dscale*sqrt( dsum )
526*
527* Triangularize the B-part by a QR factorization.
528* Apply transformation (from right) to A-part, giving S.
529*
530 CALL sgeqr2( m, m, tcpy, ldst, taul, work, linfo )
531 IF( linfo.NE.0 )
532 $ GO TO 70
533 CALL sorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
534 $ work, info )
535 CALL sorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop, ldst,
536 $ work, info )
537 IF( linfo.NE.0 )
538 $ GO TO 70
539*
540* Compute F-norm(S21) in BQRA21. (T21 is 0.)
541*
542 dscale = zero
543 dsum = one
544 DO 40 i = 1, n2
545 CALL slassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
546 40 CONTINUE
547 bqra21 = dscale*sqrt( dsum )
548*
549* Decide which method to use.
550* Weak stability test:
551* F-norm(S21) <= O(EPS * F-norm((S)))
552*
553 IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresha ) THEN
554 CALL slacpy( 'F', m, m, scpy, ldst, s, ldst )
555 CALL slacpy( 'F', m, m, tcpy, ldst, t, ldst )
556 CALL slacpy( 'F', m, m, ircop, ldst, ir, ldst )
557 CALL slacpy( 'F', m, m, licop, ldst, li, ldst )
558 ELSE IF( brqa21.GE.thresha ) THEN
559 GO TO 70
560 END IF
561*
562* Set lower triangle of B-part to zero
563*
564 CALL slaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
565*
566 IF( wands ) THEN
567*
568* Strong stability test:
569* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
570* and
571* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
572*
573 CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
574 $ m )
575 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
576 $ work, m )
577 CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
578 $ work( m*m+1 ), m )
579 dscale = zero
580 dsum = one
581 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
582 sa = dscale*sqrt( dsum )
583*
584 CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
585 $ m )
586 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
587 $ work, m )
588 CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
589 $ work( m*m+1 ), m )
590 dscale = zero
591 dsum = one
592 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
593 sb = dscale*sqrt( dsum )
594 strong = sa.LE.thresha .AND. sb.LE.threshb
595 IF( .NOT.strong )
596 $ GO TO 70
597*
598 END IF
599*
600* If the swap is accepted ("weakly" and "strongly"), apply the
601* transformations and set N1-by-N2 (2,1)-block to zero.
602*
603 CALL slaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
604*
605* copy back M-by-M diagonal block starting at index J1 of (A, B)
606*
607 CALL slacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
608 CALL slacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
609 CALL slaset( 'Full', ldst, ldst, zero, zero, t, ldst )
610*
611* Standardize existing 2-by-2 blocks.
612*
613 CALL slaset( 'Full', m, m, zero, zero, work, m )
614 work( 1 ) = one
615 t( 1, 1 ) = one
616 idum = lwork - m*m - 2
617 IF( n2.GT.1 ) THEN
618 CALL slagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
619 $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
620 work( m+1 ) = -work( 2 )
621 work( m+2 ) = work( 1 )
622 t( n2, n2 ) = t( 1, 1 )
623 t( 1, 2 ) = -t( 2, 1 )
624 END IF
625 work( m*m ) = one
626 t( m, m ) = one
627*
628 IF( n1.GT.1 ) THEN
629 CALL slagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
630 $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
631 $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
632 $ t( m, m-1 ) )
633 work( m*m ) = work( n2*m+n2+1 )
634 work( m*m-1 ) = -work( n2*m+n2+2 )
635 t( m, m ) = t( n2+1, n2+1 )
636 t( m-1, m ) = -t( m, m-1 )
637 END IF
638 CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
639 $ lda, zero, work( m*m+1 ), n2 )
640 CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
641 $ lda )
642 CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
643 $ ldb, zero, work( m*m+1 ), n2 )
644 CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
645 $ ldb )
646 CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
647 $ work( m*m+1 ), m )
648 CALL slacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
649 CALL sgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
650 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
651 CALL slacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
652 CALL sgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
653 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
654 CALL slacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
655 CALL sgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
656 $ work, m )
657 CALL slacpy( 'Full', m, m, work, m, ir, ldst )
658*
659* Accumulate transformations into Q and Z if requested.
660*
661 IF( wantq ) THEN
662 CALL sgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
663 $ ldst, zero, work, n )
664 CALL slacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
665*
666 END IF
667*
668 IF( wantz ) THEN
669 CALL sgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
670 $ ldst, zero, work, n )
671 CALL slacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
672*
673 END IF
674*
675* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
676* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
677*
678 i = j1 + m
679 IF( i.LE.n ) THEN
680 CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
681 $ a( j1, i ), lda, zero, work, m )
682 CALL slacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
683 CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
684 $ b( j1, i ), ldb, zero, work, m )
685 CALL slacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
686 END IF
687 i = j1 - 1
688 IF( i.GT.0 ) THEN
689 CALL sgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
690 $ ldst, zero, work, i )
691 CALL slacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
692 CALL sgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
693 $ ldst, zero, work, i )
694 CALL slacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
695 END IF
696*
697* Exit with INFO = 0 if swap was successfully performed.
698*
699 RETURN
700*
701 END IF
702*
703* Exit with INFO = 1 if swap was rejected.
704*
705 70 CONTINUE
706*
707 info = 1
708 RETURN
709*
710* End of STGEX2
711*
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
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:130
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:123
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
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:157
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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:110
subroutine slassq(n, x, incx, scale, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:124
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
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:274
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:114
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:114
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:159
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:159
Here is the call graph for this function:
Here is the caller graph for this function: