LAPACK 3.12.1
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 216 of file dtgexc.f.

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