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

◆ dtgexc()

subroutine dtgexc ( logical  wantq,
logical  wantz,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision, dimension( ldq, * )  q,
integer  ldq,
double precision, dimension( ldz, * )  z,
integer  ldz,
integer  ifst,
integer  ilst,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DTGEXC

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

Purpose:
 DTGEXC 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 DGGES), 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dtgexc.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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
231 $ WORK( * ), Z( LDZ, * )
232* ..
233*
234* =====================================================================
235*
236* .. Parameters ..
237 DOUBLE PRECISION ZERO
238 parameter( zero = 0.0d+0 )
239* ..
240* .. Local Scalars ..
241 LOGICAL LQUERY
242 INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
243* ..
244* .. External Subroutines ..
245 EXTERNAL dtgex2, 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( 'DTGEXC', -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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 dtgex2( 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 DTGEXC
540*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dtgex2(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, n1, n2, work, lwork, info)
DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equ...
Definition dtgex2.f:221
Here is the call graph for this function:
Here is the caller graph for this function: