 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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.```
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
```   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: