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

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