LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine stgexc ( 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  IFST,
integer  ILST,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

STGEXC

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

Purpose: 
 STGEXC reorders the generalized real Schur decomposition of a real
 matrix pair (A,B) using an orthogonal equivalence transformation

                (A, B) = Q * (A, B) * Z**T,

 so that the diagonal block of (A, B) with row index IFST is moved
 to row ILST.

 (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 generalized real Schur canonical
          form.
          On exit, the updated matrix A, again in generalized
          real Schur canonical form.
[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 generalized real Schur canonical
          form (A,B).
          On exit, the updated matrix B, again in generalized
          real Schur canonical form (A,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 (LDZ,N)
          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
          On exit, the updated matrix Q.
          If WANTQ = .FALSE., Q is not referenced.
[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.
          If WANTZ = .FALSE., Z is not referenced.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z. LDZ >= 1.
          If WANTZ = .TRUE., LDZ >= N.
[in,out]IFST
          IFST is INTEGER
[in,out]ILST
          ILST is INTEGER
          Specify the reordering of the diagonal blocks of (A, B).
          The block with row index IFST is moved to row ILST, by a
          sequence of swapping between adjacent blocks.
          On exit, if IFST pointed on entry to the second row of
          a 2-by-2 block, it is changed to point to the first row;
          ILST always points to the first row of the block in its
          final position (which may differ from its input value by
          +1 or -1). 1 <= IFST, ILST <= N.
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
           =0:  successful exit.
           <0:  if INFO = -i, the i-th argument had an illegal value.
           =1:  The transformed matrix pair (A, B) would be too far
                from generalized Schur form; the problem is ill-
                conditioned. (A, B) may have been partially reordered,
                and ILST points to the first row of the current
                position of the block being moved.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
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.

Definition at line 222 of file stgexc.f.

222 *
223 * -- LAPACK computational routine (version 3.4.0) --
224 * -- LAPACK is a software package provided by Univ. of Tennessee, --
225 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226 * November 2011
227 *
228 * .. Scalar Arguments ..
229  LOGICAL wantq, wantz
230  INTEGER ifst, ilst, info, lda, ldb, ldq, ldz, lwork, n
231 * ..
232 * .. Array Arguments ..
233  REAL a( lda, * ), b( ldb, * ), q( ldq, * ),
234  $ work( * ), z( ldz, * )
235 * ..
236 *
237 * =====================================================================
238 *
239 * .. Parameters ..
240  REAL zero
241  parameter ( zero = 0.0e+0 )
242 * ..
243 * .. Local Scalars ..
244  LOGICAL lquery
245  INTEGER here, lwmin, nbf, nbl, nbnext
246 * ..
247 * .. External Subroutines ..
248  EXTERNAL stgex2, xerbla
249 * ..
250 * .. Intrinsic Functions ..
251  INTRINSIC max
252 * ..
253 * .. Executable Statements ..
254 *
255 * Decode and test input arguments.
256 *
257  info = 0
258  lquery = ( lwork.EQ.-1 )
259  IF( n.LT.0 ) THEN
260  info = -3
261  ELSE IF( lda.LT.max( 1, n ) ) THEN
262  info = -5
263  ELSE IF( ldb.LT.max( 1, n ) ) THEN
264  info = -7
265  ELSE IF( ldq.LT.1 .OR. wantq .AND. ( ldq.LT.max( 1, n ) ) ) THEN
266  info = -9
267  ELSE IF( ldz.LT.1 .OR. wantz .AND. ( ldz.LT.max( 1, n ) ) ) THEN
268  info = -11
269  ELSE IF( ifst.LT.1 .OR. ifst.GT.n ) THEN
270  info = -12
271  ELSE IF( ilst.LT.1 .OR. ilst.GT.n ) THEN
272  info = -13
273  END IF
274 *
275  IF( info.EQ.0 ) THEN
276  IF( n.LE.1 ) THEN
277  lwmin = 1
278  ELSE
279  lwmin = 4*n + 16
280  END IF
281  work(1) = lwmin
282 *
283  IF (lwork.LT.lwmin .AND. .NOT.lquery) THEN
284  info = -15
285  END IF
286  END IF
287 *
288  IF( info.NE.0 ) THEN
289  CALL xerbla( 'STGEXC', -info )
290  RETURN
291  ELSE IF( lquery ) THEN
292  RETURN
293  END IF
294 *
295 * Quick return if possible
296 *
297  IF( n.LE.1 )
298  $ RETURN
299 *
300 * Determine the first row of the specified block and find out
301 * if it is 1-by-1 or 2-by-2.
302 *
303  IF( ifst.GT.1 ) THEN
304  IF( a( ifst, ifst-1 ).NE.zero )
305  $ ifst = ifst - 1
306  END IF
307  nbf = 1
308  IF( ifst.LT.n ) THEN
309  IF( a( ifst+1, ifst ).NE.zero )
310  $ nbf = 2
311  END IF
312 *
313 * Determine the first row of the final block
314 * and find out if it is 1-by-1 or 2-by-2.
315 *
316  IF( ilst.GT.1 ) THEN
317  IF( a( ilst, ilst-1 ).NE.zero )
318  $ ilst = ilst - 1
319  END IF
320  nbl = 1
321  IF( ilst.LT.n ) THEN
322  IF( a( ilst+1, ilst ).NE.zero )
323  $ nbl = 2
324  END IF
325  IF( ifst.EQ.ilst )
326  $ RETURN
327 *
328  IF( ifst.LT.ilst ) THEN
329 *
330 * Update ILST.
331 *
332  IF( nbf.EQ.2 .AND. nbl.EQ.1 )
333  $ ilst = ilst - 1
334  IF( nbf.EQ.1 .AND. nbl.EQ.2 )
335  $ ilst = ilst + 1
336 *
337  here = ifst
338 *
339  10 CONTINUE
340 *
341 * Swap with next one below.
342 *
343  IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
344 *
345 * Current block either 1-by-1 or 2-by-2.
346 *
347  nbnext = 1
348  IF( here+nbf+1.LE.n ) THEN
349  IF( a( here+nbf+1, here+nbf ).NE.zero )
350  $ nbnext = 2
351  END IF
352  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
353  $ ldz, here, nbf, nbnext, work, lwork, info )
354  IF( info.NE.0 ) THEN
355  ilst = here
356  RETURN
357  END IF
358  here = here + nbnext
359 *
360 * Test if 2-by-2 block breaks into two 1-by-1 blocks.
361 *
362  IF( nbf.EQ.2 ) THEN
363  IF( a( here+1, here ).EQ.zero )
364  $ nbf = 3
365  END IF
366 *
367  ELSE
368 *
369 * Current block consists of two 1-by-1 blocks, each of which
370 * must be swapped individually.
371 *
372  nbnext = 1
373  IF( here+3.LE.n ) THEN
374  IF( a( here+3, here+2 ).NE.zero )
375  $ nbnext = 2
376  END IF
377  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
378  $ ldz, here+1, 1, nbnext, work, lwork, info )
379  IF( info.NE.0 ) THEN
380  ilst = here
381  RETURN
382  END IF
383  IF( nbnext.EQ.1 ) THEN
384 *
385 * Swap two 1-by-1 blocks.
386 *
387  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
388  $ ldz, here, 1, 1, work, lwork, info )
389  IF( info.NE.0 ) THEN
390  ilst = here
391  RETURN
392  END IF
393  here = here + 1
394 *
395  ELSE
396 *
397 * Recompute NBNEXT in case of 2-by-2 split.
398 *
399  IF( a( here+2, here+1 ).EQ.zero )
400  $ nbnext = 1
401  IF( nbnext.EQ.2 ) THEN
402 *
403 * 2-by-2 block did not split.
404 *
405  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
406  $ z, ldz, here, 1, nbnext, work, lwork,
407  $ info )
408  IF( info.NE.0 ) THEN
409  ilst = here
410  RETURN
411  END IF
412  here = here + 2
413  ELSE
414 *
415 * 2-by-2 block did split.
416 *
417  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
418  $ z, ldz, here, 1, 1, work, lwork, info )
419  IF( info.NE.0 ) THEN
420  ilst = here
421  RETURN
422  END IF
423  here = here + 1
424  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
425  $ z, ldz, here, 1, 1, work, lwork, info )
426  IF( info.NE.0 ) THEN
427  ilst = here
428  RETURN
429  END IF
430  here = here + 1
431  END IF
432 *
433  END IF
434  END IF
435  IF( here.LT.ilst )
436  $ GO TO 10
437  ELSE
438  here = ifst
439 *
440  20 CONTINUE
441 *
442 * Swap with next one below.
443 *
444  IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
445 *
446 * Current block either 1-by-1 or 2-by-2.
447 *
448  nbnext = 1
449  IF( here.GE.3 ) THEN
450  IF( a( here-1, here-2 ).NE.zero )
451  $ nbnext = 2
452  END IF
453  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
454  $ ldz, here-nbnext, nbnext, nbf, work, lwork,
455  $ info )
456  IF( info.NE.0 ) THEN
457  ilst = here
458  RETURN
459  END IF
460  here = here - nbnext
461 *
462 * Test if 2-by-2 block breaks into two 1-by-1 blocks.
463 *
464  IF( nbf.EQ.2 ) THEN
465  IF( a( here+1, here ).EQ.zero )
466  $ nbf = 3
467  END IF
468 *
469  ELSE
470 *
471 * Current block consists of two 1-by-1 blocks, each of which
472 * must be swapped individually.
473 *
474  nbnext = 1
475  IF( here.GE.3 ) THEN
476  IF( a( here-1, here-2 ).NE.zero )
477  $ nbnext = 2
478  END IF
479  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
480  $ ldz, here-nbnext, nbnext, 1, work, lwork,
481  $ info )
482  IF( info.NE.0 ) THEN
483  ilst = here
484  RETURN
485  END IF
486  IF( nbnext.EQ.1 ) THEN
487 *
488 * Swap two 1-by-1 blocks.
489 *
490  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
491  $ ldz, here, nbnext, 1, work, lwork, info )
492  IF( info.NE.0 ) THEN
493  ilst = here
494  RETURN
495  END IF
496  here = here - 1
497  ELSE
498 *
499 * Recompute NBNEXT in case of 2-by-2 split.
500 *
501  IF( a( here, here-1 ).EQ.zero )
502  $ nbnext = 1
503  IF( nbnext.EQ.2 ) THEN
504 *
505 * 2-by-2 block did not split.
506 *
507  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
508  $ z, ldz, here-1, 2, 1, work, lwork, info )
509  IF( info.NE.0 ) THEN
510  ilst = here
511  RETURN
512  END IF
513  here = here - 2
514  ELSE
515 *
516 * 2-by-2 block did split.
517 *
518  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
519  $ z, ldz, here, 1, 1, work, lwork, info )
520  IF( info.NE.0 ) THEN
521  ilst = here
522  RETURN
523  END IF
524  here = here - 1
525  CALL stgex2( wantq, wantz, n, a, lda, b, ldb, q, ldq,
526  $ z, ldz, here, 1, 1, work, lwork, info )
527  IF( info.NE.0 ) THEN
528  ilst = here
529  RETURN
530  END IF
531  here = here - 1
532  END IF
533  END IF
534  END IF
535  IF( here.GT.ilst )
536  $ GO TO 20
537  END IF
538  ilst = here
539  work( 1 ) = lwmin
540  RETURN
541 *
542 * End of STGEXC
543 *
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
stgex2
subroutine stgex2(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO)
STGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition: stgex2.f:223

Here is the call graph for this function:

Here is the caller graph for this function: