Purpose
To compute orthogonal transformations Q and Z such that the
transformed pencil Q'(sE-A)Z is in upper block triangular form,
where E is an M-by-N matrix in column echelon form (see SLICOT
Library routine MB04UD) and A is an M-by-N matrix.
If MODE = 'B', then the matrices A and E are transformed into the
following generalized Schur form by unitary transformations Q1
and Z1 :
| sE(eps,inf)-A(eps,inf) | X |
Q1'(sE-A)Z1 = |------------------------|------------|. (1)
| O | sE(r)-A(r) |
The pencil sE(eps,inf)-A(eps,inf) is in staircase form, and it
contains all Kronecker column indices and infinite elementary
divisors of the pencil sE-A. The pencil sE(r)-A(r) contains all
Kronecker row indices and elementary divisors of sE-A.
Note: X is a pencil.
If MODE = 'T', then the submatrices having full row and column
rank in the pencil sE(eps,inf)-A(eps,inf) in (1) are
triangularized by applying unitary transformations Q2 and Z2 to
Q1'*(sE-A)*Z1.
If MODE = 'S', then the pencil sE(eps,inf)-A(eps,inf) in (1) is
separated into sE(eps)-A(eps) and sE(inf)-A(inf) by applying
unitary transformations Q3 and Z3 to Q2'*Q1'*(sE-A)*Z1*Z2.
This gives
| sE(eps)-A(eps) | X | X |
|----------------|----------------|------------|
| O | sE(inf)-A(inf) | X |
Q'(sE-A)Z =|=================================|============| (2)
| | |
| O | sE(r)-A(r) |
where Q = Q1*Q2*Q3 and Z = Z1*Z2*Z3.
Note: the pencil sE(r)-A(r) is not reduced further.
Specification
SUBROUTINE MB04VD( MODE, JOBQ, JOBZ, M, N, RANKE, A, LDA, E, LDE,
$ Q, LDQ, Z, LDZ, ISTAIR, NBLCKS, NBLCKI, IMUK,
$ INUK, IMUK0, MNEI, TOL, IWORK, INFO )
C .. Scalar Arguments ..
CHARACTER JOBQ, JOBZ, MODE
INTEGER INFO, LDA, LDE, LDQ, LDZ, M, N, NBLCKI, NBLCKS,
$ RANKE
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IMUK(*), IMUK0(*), INUK(*), ISTAIR(*), IWORK(*),
$ MNEI(*)
DOUBLE PRECISION A(LDA,*), E(LDE,*), Q(LDQ,*), Z(LDZ,*)
Arguments
Mode Parameters
MODE CHARACTER*1
Specifies the desired structure of the transformed
pencil Q'(sE-A)Z to be computed as follows:
= 'B': Basic reduction given by (1);
= 'T': Further reduction of (1) to triangular form;
= 'S': Further separation of sE(eps,inf)-A(eps,inf)
in (1) into the two pencils in (2).
JOBQ CHARACTER*1
Indicates whether the user wishes to accumulate in a
matrix Q the orthogonal row transformations, as follows:
= 'N': Do not form Q;
= 'I': Q is initialized to the unit matrix and the
orthogonal transformation matrix Q is returned;
= 'U': The given matrix Q is updated by the orthogonal
row transformations used in the reduction.
JOBZ CHARACTER*1
Indicates whether the user wishes to accumulate in a
matrix Z the orthogonal column transformations, as
follows:
= 'N': Do not form Z;
= 'I': Z is initialized to the unit matrix and the
orthogonal transformation matrix Z is returned;
= 'U': The given matrix Z is updated by the orthogonal
transformations used in the reduction.
Input/Output Parameters
M (input) INTEGER
The number of rows in the matrices A, E and the order of
the matrix Q. M >= 0.
N (input) INTEGER
The number of columns in the matrices A, E and the order
of the matrix Z. N >= 0.
RANKE (input) INTEGER
The rank of the matrix E in column echelon form.
RANKE >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading M-by-N part of this array must
contain the matrix to be row compressed.
On exit, the leading M-by-N part of this array contains
the matrix that has been row compressed while keeping
matrix E in column echelon form.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,M).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading M-by-N part of this array must
contain the matrix in column echelon form to be
transformed equivalent to matrix A.
On exit, the leading M-by-N part of this array contains
the matrix that has been transformed equivalent to matrix
A.
LDE INTEGER
The leading dimension of array E. LDE >= MAX(1,M).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,*)
On entry, if JOBQ = 'U', then the leading M-by-M part of
this array must contain a given matrix Q (e.g. from a
previous call to another SLICOT routine), and on exit, the
leading M-by-M part of this array contains the product of
the input matrix Q and the row transformation matrix used
to transform the rows of matrices A and E.
On exit, if JOBQ = 'I', then the leading M-by-M part of
this array contains the matrix of accumulated orthogonal
row transformations performed.
If JOBQ = 'N', the array Q is not referenced and can be
supplied as a dummy array (i.e. set parameter LDQ = 1 and
declare this array to be Q(1,1) in the calling program).
LDQ INTEGER
The leading dimension of array Q. If JOBQ = 'U' or
JOBQ = 'I', LDQ >= MAX(1,M); if JOBQ = 'N', LDQ >= 1.
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,*)
On entry, if JOBZ = 'U', then the leading N-by-N part of
this array must contain a given matrix Z (e.g. from a
previous call to another SLICOT routine), and on exit, the
leading N-by-N part of this array contains the product of
the input matrix Z and the column transformation matrix
used to transform the columns of matrices A and E.
On exit, if JOBZ = 'I', then the leading N-by-N part of
this array contains the matrix of accumulated orthogonal
column transformations performed.
If JOBZ = 'N', the array Z is not referenced and can be
supplied as a dummy array (i.e. set parameter LDZ = 1 and
declare this array to be Z(1,1) in the calling program).
LDZ INTEGER
The leading dimension of array Z. If JOBZ = 'U' or
JOBZ = 'I', LDZ >= MAX(1,N); if JOBZ = 'N', LDZ >= 1.
ISTAIR (input/output) INTEGER array, dimension (M)
On entry, this array must contain information on the
column echelon form of the unitary transformed matrix E.
Specifically, ISTAIR(i) must be set to +j if the first
non-zero element E(i,j) is a corner point and -j
otherwise, for i = 1,2,...,M.
On exit, this array contains no useful information.
NBLCKS (output) INTEGER
The number of submatrices having full row rank greater
than or equal to 0 detected in matrix A in the pencil
sE(x)-A(x),
where x = eps,inf if MODE = 'B' or 'T',
or x = eps if MODE = 'S'.
NBLCKI (output) INTEGER
If MODE = 'S', the number of diagonal submatrices in the
pencil sE(inf)-A(inf). If MODE = 'B' or 'T' then
NBLCKI = 0.
IMUK (output) INTEGER array, dimension (MAX(N,M+1))
The leading NBLCKS elements of this array contain the
column dimensions mu(1),...,mu(NBLCKS) of the submatrices
having full column rank in the pencil sE(x)-A(x),
where x = eps,inf if MODE = 'B' or 'T',
or x = eps if MODE = 'S'.
INUK (output) INTEGER array, dimension (MAX(N,M+1))
The leading NBLCKS elements of this array contain the
row dimensions nu(1),...,nu(NBLCKS) of the submatrices
having full row rank in the pencil sE(x)-A(x),
where x = eps,inf if MODE = 'B' or 'T',
or x = eps if MODE = 'S'.
IMUK0 (output) INTEGER array, dimension (limuk0),
where limuk0 = N if MODE = 'S' and 1, otherwise.
If MODE = 'S', then the leading NBLCKI elements of this
array contain the dimensions mu0(1),...,mu0(NBLCKI)
of the square diagonal submatrices in the pencil
sE(inf)-A(inf).
Otherwise, IMUK0 is not referenced and can be supplied
as a dummy array.
MNEI (output) INTEGER array, dimension (3)
If MODE = 'B' or 'T' then
MNEI(1) contains the row dimension of
sE(eps,inf)-A(eps,inf);
MNEI(2) contains the column dimension of
sE(eps,inf)-A(eps,inf);
MNEI(3) = 0.
If MODE = 'S', then
MNEI(1) contains the row dimension of sE(eps)-A(eps);
MNEI(2) contains the column dimension of sE(eps)-A(eps);
MNEI(3) contains the order of the regular pencil
sE(inf)-A(inf).
Tolerances
TOL DOUBLE PRECISION
A tolerance below which matrix elements are considered
to be zero. If the user sets TOL to be less than (or
equal to) zero then the tolerance is taken as
EPS * MAX( ABS(A(I,J)), ABS(E(I,J)) ), where EPS is the
machine precision (see LAPACK Library routine DLAMCH),
I = 1,2,...,M and J = 1,2,...,N.
Workspace
IWORK INTEGER array, dimension (N)Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value.
> 0: if incorrect rank decisions were revealed during the
triangularization phase. This failure is not likely
to occur. The possible values are:
= 1: if incorrect dimensions of a full column rank
submatrix;
= 2: if incorrect dimensions of a full row rank
submatrix.
Method
Let sE - A be an arbitrary pencil. Prior to calling the routine, this pencil must be transformed into a pencil with E in column echelon form. This may be accomplished by calling the SLICOT Library routine MB04UD. Depending on the value of MODE, submatrices of A and E are then reduced to one of the forms described above. Further details can be found in [1].References
[1] Beelen, Th. and Van Dooren, P.
An improved algorithm for the computation of Kronecker's
canonical form of a singular pencil.
Linear Algebra and Applications, 105, pp. 9-65, 1988.
Numerical Aspects
It is shown in [1] that the algorithm is numerically backward stable. The operations count is proportional to (MAX(M,N))**3.Further Comments
The difference mu(k)-nu(k), for k = 1,2,...,NBLCKS, is the number of elementary Kronecker blocks of size k x (k+1). If MODE = 'B' or 'T' on entry, then the difference nu(k)-mu(k+1), for k = 1,2,...,NBLCKS, is the number of infinite elementary divisors of degree k (with mu(NBLCKS+1) = 0). If MODE = 'S' on entry, then the difference mu0(k)-mu0(k+1), for k = 1,2,...,NBLCKI, is the number of infinite elementary divisors of degree k (with mu0(NBLCKI+1) = 0). In the pencil sE(r)-A(r), the pencils sE(f)-A(f) and sE(eta)-A(eta) can be separated by pertransposing the pencil sE(r)-A(r) and calling the routine with MODE set to 'B'. The result has got to be pertransposed again. (For more details see [1]).Example
Program Text
* MB04VD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER MMAX, NMAX
PARAMETER ( MMAX = 20, NMAX = 20 )
INTEGER LDA, LDE, LDQ, LDZ
PARAMETER ( LDA = MMAX, LDE = MMAX, LDQ = MMAX,
$ LDZ = NMAX )
INTEGER LINUK
PARAMETER ( LINUK = MAX( NMAX,MMAX+1 ) )
* PARAMETER ( LINUK = NMAX+MMAX+1 )
INTEGER LIWORK
PARAMETER ( LIWORK = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( NMAX,MMAX ) )
* PARAMETER ( LDWORK = NMAX+MMAX )
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
* .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER I, INFO, J, M, N, NBLCKI, NBLCKS, RANKE
CHARACTER*1 JOBQ, JOBZ, MODE
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK), E(LDE,NMAX),
$ Q(LDQ,MMAX), Z(LDZ,NMAX)
INTEGER IMUK(LINUK), IMUK0(NMAX), INUK(LINUK),
$ ISTAIR(MMAX), IWORK(LIWORK), MNEI(3)
C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL MB04UD, MB04VD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) M, N, TOL, MODE
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99984 ) M
ELSE IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99983 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,M )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,M )
JOBQ = 'I'
JOBZ = 'I'
* Reduce E to column echelon form and compute Q'*A*Z.
CALL MB04UD( JOBQ, JOBZ, M, N, A, LDA, E, LDE, Q, LDQ, Z, LDZ,
$ RANKE, ISTAIR, TOL, DWORK, INFO )
JOBQ = 'U'
JOBZ = 'U'
*
IF ( INFO.EQ.0 ) THEN
* Compute a unitary transformed pencil Q'*(s*E-A)*Z.
CALL MB04VD( MODE, JOBQ, JOBZ, M, N, RANKE, A, LDA, E, LDE,
$ Q, LDQ, Z, LDZ, ISTAIR, NBLCKS, NBLCKI, IMUK,
$ INUK, IMUK0, MNEI, TOL, IWORK, INFO )
*
IF ( INFO.EQ.0 ) THEN
WRITE ( NOUT, FMT = 99996 )
WRITE ( NOUT, FMT = 99995 )
DO 140 I = 1, M
WRITE ( NOUT, FMT = 99994 ) ( Q(I,J), J = 1,M )
140 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 160 I = 1, M
WRITE ( NOUT, FMT = 99994 ) ( E(I,J), J = 1,N )
160 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 180 I = 1, M
WRITE ( NOUT, FMT = 99994 ) ( A(I,J), J = 1,N )
180 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 200 I = 1, N
WRITE ( NOUT, FMT = 99994 ) ( Z(I,J), J = 1,N )
200 CONTINUE
WRITE ( NOUT, FMT = 99990 ) NBLCKS
IF ( .NOT. LSAME( MODE, 'S' ) ) THEN
WRITE ( NOUT, FMT = 99989 ) ( IMUK(I), I = 1,NBLCKS )
WRITE ( NOUT, FMT = 99988 ) ( INUK(I), I = 1,NBLCKS )
ELSE
WRITE ( NOUT, FMT = 99987 ) ( IMUK(I), I = 1,NBLCKS )
WRITE ( NOUT, FMT = 99986 ) ( INUK(I), I = 1,NBLCKS )
WRITE ( NOUT, FMT = 99982 ) ( IMUK0(I), I = 1,NBLCKI )
WRITE ( NOUT, FMT = 99985 ) ( MNEI(I), I = 1,3 )
END IF
ELSE
WRITE ( NOUT, FMT = 99998 ) INFO
END IF
ELSE
WRITE ( NOUT, FMT = 99997 ) INFO
END IF
END IF
STOP
*
99999 FORMAT (' MB04VD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04VD = ',I2)
99997 FORMAT (' INFO on exit from MB04UD = ',I2)
99996 FORMAT (' The unitary transformed pencil is Q''*(s*E-A)*Z, where',
$ /)
99995 FORMAT (' Matrix Q',/)
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' Matrix E',/)
99992 FORMAT (/' Matrix A',/)
99991 FORMAT (/' Matrix Z',/)
99990 FORMAT (/' The number of submatrices having full row rank detect',
$ 'ed in matrix A = ',I3)
99989 FORMAT (/' The column dimensions of the submatrices having full ',
$ 'column rank in the pencil',/' sE(eps,inf) - A(eps,inf) a',
$ 're',/20(1X,I5))
99988 FORMAT (/' The row dimensions of the submatrices having full row',
$ ' rank in the pencil',/' sE(eps,inf) - A(eps,inf) are',
$ /20(1X,I5))
99987 FORMAT (/' The column dimensions of the submatrices having full ',
$ 'column rank in the pencil',/' sE(eps) - A(eps) are',
$ /20(1X,I5))
99986 FORMAT (/' The row dimensions of the submatrices having full row',
$ ' rank in the pencil',/' sE(eps) - A(eps) are',/20(1X,I5))
99985 FORMAT (/' MNEI is ',/20(1X,I5))
99984 FORMAT (/' M is out of range.',/' M = ',I5)
99983 FORMAT (/' N is out of range.',/' N = ',I5)
99982 FORMAT (/' The orders of the diagonal submatrices in the pencil ',
$ 'sE(inf) - A(inf) are',/20(1X,I5))
END
Program Data
MB04VD EXAMPLE PROGRAM DATA 2 4 0.0 S 1.0 0.0 -1.0 0.0 1.0 1.0 0.0 -1.0 0.0 -1.0 0.0 0.0 0.0 -1.0 0.0 0.0Program Results
MB04VD EXAMPLE PROGRAM RESULTS
The unitary transformed pencil is Q'*(s*E-A)*Z, where
Matrix Q
0.7071 -0.7071
0.7071 0.7071
Matrix E
0.0000 0.0000 -1.1547 0.8165
0.0000 0.0000 0.0000 0.0000
Matrix A
0.0000 1.7321 0.5774 -0.4082
0.0000 0.0000 0.0000 -1.2247
Matrix Z
0.5774 0.8165 0.0000 0.0000
0.0000 0.0000 0.8165 -0.5774
0.5774 -0.4082 -0.4082 -0.5774
0.5774 -0.4082 0.4082 0.5774
The number of submatrices having full row rank detected in matrix A = 2
The column dimensions of the submatrices having full column rank in the pencil
sE(eps) - A(eps) are
2 1
The row dimensions of the submatrices having full row rank in the pencil
sE(eps) - A(eps) are
1 0
The orders of the diagonal submatrices in the pencil sE(inf) - A(inf) are
1
MNEI is
1 3 1