**Purpose**

To reorder the diagonal blocks of a principal submatrix of an upper quasi-triangular matrix A together with their eigenvalues by constructing an orthogonal similarity transformation UT. After reordering, the leading block of the selected submatrix of A has eigenvalues in a suitably defined domain of interest, usually related to stability/instability in a continuous- or discrete-time sense.

SUBROUTINE MB03QD( DICO, STDOM, JOBU, N, NLOW, NSUP, ALPHA, $ A, LDA, U, LDU, NDIM, DWORK, INFO ) C .. Scalar Arguments .. CHARACTER DICO, JOBU, STDOM INTEGER INFO, LDA, LDU, N, NDIM, NLOW, NSUP DOUBLE PRECISION ALPHA C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), DWORK(*), U(LDU,*)

**Mode Parameters**

DICO CHARACTER*1 Specifies the type of the spectrum separation to be performed as follows: = 'C': continuous-time sense; = 'D': discrete-time sense. STDOM CHARACTER*1 Specifies whether the domain of interest is of stability type (left part of complex plane or inside of a circle) or of instability type (right part of complex plane or outside of a circle) as follows: = 'S': stability type domain; = 'U': instability type domain. JOBU CHARACTER*1 Indicates how the performed orthogonal transformations UT are accumulated, as follows: = 'I': U is initialized to the unit matrix and the matrix UT is returned in U; = 'U': the given matrix U is updated and the matrix U*UT is returned in U.

N (input) INTEGER The order of the matrices A and U. N >= 1. NLOW, (input) INTEGER NSUP NLOW and NSUP specify the boundary indices for the rows and columns of the principal submatrix of A whose diagonal blocks are to be reordered. 1 <= NLOW <= NSUP <= N. ALPHA (input) DOUBLE PRECISION The boundary of the domain of interest for the eigenvalues of A. If DICO = 'C', ALPHA is the boundary value for the real parts of eigenvalues, while for DICO = 'D', ALPHA >= 0 represents the boundary value for the moduli of eigenvalues. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading N-by-N part of this array must contain a matrix in a real Schur form whose 1-by-1 and 2-by-2 diagonal blocks between positions NLOW and NSUP are to be reordered. On exit, the leading N-by-N part contains the ordered real Schur matrix UT' * A * UT with the elements below the first subdiagonal set to zero. The leading NDIM-by-NDIM part of the principal submatrix D = A(NLOW:NSUP,NLOW:NSUP) has eigenvalues in the domain of interest and the trailing part of this submatrix has eigenvalues outside the domain of interest. The domain of interest for lambda(D), the eigenvalues of D, is defined by the parameters ALPHA, DICO and STDOM as follows: For DICO = 'C': Real(lambda(D)) < ALPHA if STDOM = 'S'; Real(lambda(D)) > ALPHA if STDOM = 'U'. For DICO = 'D': Abs(lambda(D)) < ALPHA if STDOM = 'S'; Abs(lambda(D)) > ALPHA if STDOM = 'U'. LDA INTEGER The leading dimension of array A. LDA >= N. U (input/output) DOUBLE PRECISION array, dimension (LDU,N) On entry with JOBU = 'U', the leading N-by-N part of this array must contain a transformation matrix (e.g. from a previous call to this routine). On exit, if JOBU = 'U', the leading N-by-N part of this array contains the product of the input matrix U and the orthogonal matrix UT used to reorder the diagonal blocks of A. On exit, if JOBU = 'I', the leading N-by-N part of this array contains the matrix UT of the performed orthogonal transformations. Array U need not be set on entry if JOBU = 'I'. LDU INTEGER The leading dimension of array U. LDU >= N. NDIM (output) INTEGER The number of eigenvalues of the selected principal submatrix lying inside the domain of interest. If NLOW = 1, NDIM is also the dimension of the invariant subspace corresponding to the eigenvalues of the leading NDIM-by-NDIM submatrix. In this case, if U is the orthogonal transformation matrix used to compute and reorder the real Schur form of A, its first NDIM columns form an orthonormal basis for the above invariant subspace.

DWORK DOUBLE PRECISION array, dimension (N)

INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = 1: A(NLOW,NLOW-1) is nonzero, i.e. A(NLOW,NLOW) is not the leading element of a 1-by-1 or 2-by-2 diagonal block of A, or A(NSUP+1,NSUP) is nonzero, i.e. A(NSUP,NSUP) is not the bottom element of a 1-by-1 or 2-by-2 diagonal block of A; = 2: two adjacent blocks are too close to swap (the problem is very ill-conditioned).

Given an upper quasi-triangular matrix A with 1-by-1 or 2-by-2 diagonal blocks, the routine reorders its diagonal blocks along with its eigenvalues by performing an orthogonal similarity transformation UT' * A * UT. The column transformation UT is also performed on the given (initial) transformation U (resulted from a possible previous step or initialized as the identity matrix). After reordering, the eigenvalues inside the region specified by the parameters ALPHA, DICO and STDOM appear at the top of the selected diagonal block between positions NLOW and NSUP. In other words, lambda(A(NLOW:NSUP,NLOW:NSUP)) are ordered such that lambda(A(NLOW:NLOW+NDIM-1,NLOW:NLOW+NDIM-1)) are inside and lambda(A(NLOW+NDIM:NSUP,NLOW+NDIM:NSUP)) are outside the domain of interest. If NLOW = 1, the first NDIM columns of U*UT span the corresponding invariant subspace of A.

[1] Stewart, G.W. HQR3 and EXCHQZ: FORTRAN subroutines for calculating and ordering the eigenvalues of a real upper Hessenberg matrix. ACM TOMS, 2, pp. 275-280, 1976.

3 The algorithm requires less than 4*N operations.

None

**Program Text**

* MB03QD EXAMPLE PROGRAM TEXT * Copyright (c) 2002-2017 NICONET e.V. * * .. Parameters .. INTEGER NIN, NOUT PARAMETER ( NIN = 5, NOUT = 6 ) INTEGER NMAX PARAMETER ( NMAX = 10 ) INTEGER LDA, LDU PARAMETER ( LDA = NMAX, LDU = NMAX ) INTEGER LDWORK PARAMETER ( LDWORK = 3*NMAX ) * .. Local Scalars .. CHARACTER*1 DICO, JOBU, STDOM INTEGER I, INFO, J, N, NDIM, NLOW, NSUP DOUBLE PRECISION ALPHA * .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK), U(LDU,NMAX), $ WI(NMAX), WR(NMAX) LOGICAL BWORK(NMAX) * .. External Functions .. LOGICAL SELECT * .. External Subroutines .. EXTERNAL DGEES, MB03QD * .. Executable Statements .. * WRITE ( NOUT, FMT = 99999 ) * Skip the heading in the data file and read the data. READ ( NIN, FMT = '()' ) READ ( NIN, FMT = * ) N, NLOW, NSUP, ALPHA, DICO, STDOM, JOBU IF ( N.LT.0 .OR. N.GT.NMAX ) THEN WRITE ( NOUT, FMT = 99992 ) N ELSE READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N ) * Compute Schur form, eigenvalues and Schur vectors. CALL DGEES( 'Vectors', 'Not sorted', SELECT, N, A, LDA, NDIM, $ WR, WI, U, LDU, DWORK, LDWORK, BWORK, INFO ) IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99998 ) INFO ELSE * Block reordering. CALL MB03QD( DICO, STDOM, JOBU, N, NLOW, NSUP, ALPHA, $ A, LDA, U, LDU, NDIM, DWORK, INFO ) IF ( INFO.NE.0 ) THEN WRITE ( NOUT, FMT = 99997 ) INFO ELSE WRITE ( NOUT, FMT = 99996 ) NDIM WRITE ( NOUT, FMT = 99994 ) DO 10 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N ) 10 CONTINUE WRITE ( NOUT, FMT = 99993 ) DO 20 I = 1, N WRITE ( NOUT, FMT = 99995 ) ( U(I,J), J = 1,N ) 20 CONTINUE END IF END IF END IF * STOP * 99999 FORMAT (' MB03QD EXAMPLE PROGRAM RESULTS',/1X) 99998 FORMAT (' INFO on exit from DGEES = ',I2) 99997 FORMAT (' INFO on exit from MB03QD = ',I2) 99996 FORMAT (' The number of eigenvalues in the domain is ',I5) 99995 FORMAT (8X,20(1X,F8.4)) 99994 FORMAT (/' The ordered Schur form matrix is ') 99993 FORMAT (/' The transformation matrix is ') 99992 FORMAT (/' N is out of range.',/' N = ',I5) END

MB03QD EXAMPLE PROGRAM DATA 4 1 4 0.0 C S U -1.0 37.0 -12.0 -12.0 -1.0 -10.0 0.0 4.0 2.0 -4.0 7.0 -6.0 2.0 2.0 7.0 -9.0

MB03QD EXAMPLE PROGRAM RESULTS The number of eigenvalues in the domain is 4 The ordered Schur form matrix is -3.1300 -26.5066 27.2262 -16.2009 0.9070 -3.1300 13.6254 8.9206 0.0000 0.0000 -3.3700 0.3419 0.0000 0.0000 -1.7879 -3.3700 The transformation matrix is 0.9611 0.1784 0.2064 -0.0440 -0.1468 -0.2704 0.8116 -0.4965 -0.2224 0.7675 0.4555 0.3924 -0.0733 0.5531 -0.3018 -0.7730