LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2017
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 stgex2.f.

223 *
224 * -- LAPACK auxiliary routine (version 3.7.1) --
225 * -- LAPACK is a software package provided by Univ. of Tennessee, --
226 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227 * June 2017
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  REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
235  $ WORK( * ), Z( LDZ, * )
236 * ..
237 *
238 * =====================================================================
239 * Replaced various illegal calls to SCOPY by calls to SLASET, or by DO
240 * loops. Sven Hammarling, 1/5/02.
241 *
242 * .. Parameters ..
243  REAL ZERO, ONE
244  parameter( zero = 0.0e+0, one = 1.0e+0 )
245  REAL TWENTY
246  parameter( twenty = 2.0e+01 )
247  INTEGER LDST
248  parameter( ldst = 4 )
249  LOGICAL WANDS
250  parameter( wands = .true. )
251 * ..
252 * .. Local Scalars ..
253  LOGICAL STRONG, WEAK
254  INTEGER I, IDUM, LINFO, M
255  REAL 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  REAL 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  REAL SLAMCH
268  EXTERNAL slamch
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL sgemm, sgeqr2, sgerq2, slacpy, slagv2, slartg,
272  $ slaset, slassq, sorg2r, sorgr2, sorm2r, sormr2,
273  $ srot, sscal, stgsy2
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( n*m, m*m*2 ) ) THEN
290  info = -16
291  work( 1 ) = max( n*m, m*m*2 )
292  RETURN
293  END IF
294 *
295  weak = .false.
296  strong = .false.
297 *
298 * Make a local copy of selected block
299 *
300  CALL slaset( 'Full', ldst, ldst, zero, zero, li, ldst )
301  CALL slaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
302  CALL slacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
303  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
304 *
305 * Compute threshold for testing acceptance of swapping.
306 *
307  eps = slamch( 'P' )
308  smlnum = slamch( 'S' ) / eps
309  dscale = zero
310  dsum = one
311  CALL slacpy( 'Full', m, m, s, ldst, work, m )
312  CALL slassq( m*m, work, 1, dscale, dsum )
313  CALL slacpy( 'Full', m, m, t, ldst, work, m )
314  CALL slassq( 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 slartg( 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 srot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
342  $ ir( 2, 1 ) )
343  CALL srot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
344  $ ir( 2, 1 ) )
345  IF( sa.GE.sb ) THEN
346  CALL slartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
347  $ ddum )
348  ELSE
349  CALL slartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350  $ ddum )
351  END IF
352  CALL srot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
353  $ li( 2, 1 ) )
354  CALL srot( 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 slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
373  $ m )
374  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
375  $ work, m )
376  CALL sgemm( '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 slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
381 *
382  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
383  $ m )
384  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
385  $ work, m )
386  CALL sgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
387  $ work( m*m+1 ), m )
388  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389  ss = dscale*sqrt( dsum )
390  strong = ss.LE.thresh
391  IF( .NOT.strong )
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 srot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
399  $ ir( 2, 1 ) )
400  CALL srot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
401  $ ir( 2, 1 ) )
402  CALL srot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403  $ li( 1, 1 ), li( 2, 1 ) )
404  CALL srot( 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 srot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
416  $ ir( 2, 1 ) )
417  IF( wantq )
418  $ CALL srot( 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 slacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436  CALL slacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
437  $ ir( n2+1, n1+1 ), ldst )
438  CALL stgsy2( '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 sscal( n1, -one, li( 1, i ), 1 )
453  li( n1+i, i ) = scale
454  10 CONTINUE
455  CALL sgeqr2( m, n2, li, ldst, taul, work, linfo )
456  IF( linfo.NE.0 )
457  $ GO TO 70
458  CALL sorg2r( 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 sgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
472  IF( linfo.NE.0 )
473  $ GO TO 70
474  CALL sorgr2( 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 sgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
481  $ work, m )
482  CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
483  $ ldst )
484  CALL sgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
485  $ work, m )
486  CALL sgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
487  $ ldst )
488  CALL slacpy( 'F', m, m, s, ldst, scpy, ldst )
489  CALL slacpy( 'F', m, m, t, ldst, tcpy, ldst )
490  CALL slacpy( 'F', m, m, ir, ldst, ircop, ldst )
491  CALL slacpy( '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 sgerq2( m, m, t, ldst, taur, work, linfo )
497  IF( linfo.NE.0 )
498  $ GO TO 70
499  CALL sormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
500  $ linfo )
501  IF( linfo.NE.0 )
502  $ GO TO 70
503  CALL sormr2( '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 slassq( 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 sgeqr2( m, m, tcpy, ldst, taul, work, linfo )
521  IF( linfo.NE.0 )
522  $ GO TO 70
523  CALL sorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
524  $ work, info )
525  CALL sorm2r( '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 slassq( 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 slacpy( 'F', m, m, scpy, ldst, s, ldst )
545  CALL slacpy( 'F', m, m, tcpy, ldst, t, ldst )
546  CALL slacpy( 'F', m, m, ircop, ldst, ir, ldst )
547  CALL slacpy( '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 slaset( '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 slacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
562  $ m )
563  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
564  $ work, m )
565  CALL sgemm( '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 slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
570 *
571  CALL slacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
572  $ m )
573  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
574  $ work, m )
575  CALL sgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
576  $ work( m*m+1 ), m )
577  CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578  ss = dscale*sqrt( dsum )
579  strong = ( ss.LE.thresh )
580  IF( .NOT.strong )
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 slaset( '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 slacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
593  CALL slacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
594  CALL slaset( 'Full', ldst, ldst, zero, zero, t, ldst )
595 *
596 * Standardize existing 2-by-2 blocks.
597 *
598  CALL slaset( '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 slagv2( 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 slagv2( 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 sgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
624  $ lda, zero, work( m*m+1 ), n2 )
625  CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
626  $ lda )
627  CALL sgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
628  $ ldb, zero, work( m*m+1 ), n2 )
629  CALL slacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
630  $ ldb )
631  CALL sgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
632  $ work( m*m+1 ), m )
633  CALL slacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
634  CALL sgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
635  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
636  CALL slacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
637  CALL sgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
638  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
639  CALL slacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
640  CALL sgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
641  $ work, m )
642  CALL slacpy( 'Full', m, m, work, m, ir, ldst )
643 *
644 * Accumulate transformations into Q and Z if requested.
645 *
646  IF( wantq ) THEN
647  CALL sgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
648  $ ldst, zero, work, n )
649  CALL slacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
650 *
651  END IF
652 *
653  IF( wantz ) THEN
654  CALL sgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
655  $ ldst, zero, work, n )
656  CALL slacpy( '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 sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
666  $ a( j1, i ), lda, zero, work, m )
667  CALL slacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
668  CALL sgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
669  $ b( j1, i ), ldb, zero, work, m )
670  CALL slacpy( '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 sgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
675  $ ldst, zero, work, i )
676  CALL slacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
677  CALL sgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
678  $ ldst, zero, work, i )
679  CALL slacpy( '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 STGEX2
696 *
Here is the call graph for this function:
Here is the caller graph for this function:
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:132
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:94
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:105
slassq
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f:105
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
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:116
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:116
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:112
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
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:161
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:159
slartg
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
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:276
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:161
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:125