LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ stgexc()

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 (LDQ,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.
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 218 of file stgexc.f.

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