Purpose
To compute the coefficients of a minimal polynomial basis
DK
K(s) = K(0) + K(1) * s + ... + K(DK) * s
for the right nullspace of the MP-by-NP polynomial matrix of
degree DP, given by
DP
P(s) = P(0) + P(1) * s + ... + P(DP) * s ,
which corresponds to solving the polynomial matrix equation
P(s) * K(s) = 0.
Specification
SUBROUTINE MC03ND( MP, NP, DP, P, LDP1, LDP2, DK, GAM, NULLSP,
$ LDNULL, KER, LDKER1, LDKER2, TOL, IWORK, DWORK,
$ LDWORK, INFO )
C .. Scalar Arguments ..
INTEGER DK, DP, INFO, LDKER1, LDKER2, LDNULL, LDP1,
$ LDP2, LDWORK, MP, NP
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER GAM(*), IWORK(*)
DOUBLE PRECISION DWORK(*), KER(LDKER1,LDKER2,*),
$ NULLSP(LDNULL,*), P(LDP1,LDP2,*)
Arguments
Input/Output Parameters
MP (input) INTEGER
The number of rows of the polynomial matrix P(s).
MP >= 0.
NP (input) INTEGER
The number of columns of the polynomial matrix P(s).
NP >= 0.
DP (input) INTEGER
The degree of the polynomial matrix P(s). DP >= 1.
P (input) DOUBLE PRECISION array, dimension (LDP1,LDP2,DP+1)
The leading MP-by-NP-by-(DP+1) part of this array must
contain the coefficients of the polynomial matrix P(s).
Specifically, P(i,j,k) must contain the (i,j)-th element
of P(k-1), which is the cofficient of s**(k-1) of P(s),
where i = 1,2,...,MP, j = 1,2,...,NP and k = 1,2,...,DP+1.
LDP1 INTEGER
The leading dimension of array P. LDP1 >= MAX(1,MP).
LDP2 INTEGER
The second dimension of array P. LDP2 >= MAX(1,NP).
DK (output) INTEGER
The degree of the minimal polynomial basis K(s) for the
right nullspace of P(s) unless DK = -1, in which case
there is no right nullspace.
GAM (output) INTEGER array, dimension (DP*MP+1)
The leading (DK+1) elements of this array contain
information about the ordering of the right nullspace
vectors stored in array NULLSP.
NULLSP (output) DOUBLE PRECISION array, dimension
(LDNULL,(DP*MP+1)*NP)
The leading NP-by-SUM(i*GAM(i)) part of this array
contains the right nullspace vectors of P(s) in condensed
form (as defined in METHOD), where i = 1,2,...,DK+1.
LDNULL INTEGER
The leading dimension of array NULLSP.
LDNULL >= MAX(1,NP).
KER (output) DOUBLE PRECISION array, dimension
(LDKER1,LDKER2,DP*MP+1)
The leading NP-by-nk-by-(DK+1) part of this array contains
the coefficients of the minimal polynomial basis K(s),
where nk = SUM(GAM(i)) and i = 1,2,...,DK+1. Specifically,
KER(i,j,m) contains the (i,j)-th element of K(m-1), which
is the coefficient of s**(m-1) of K(s), where i = 1,2,...,
NP, j = 1,2,...,nk and m = 1,2,...,DK+1.
LDKER1 INTEGER
The leading dimension of array KER. LDKER1 >= MAX(1,NP).
LDKER2 INTEGER
The second dimension of array KER. LDKER2 >= MAX(1,NP).
Tolerances
TOL DOUBLE PRECISION
A tolerance below which matrix elements are considered
to be zero. If the user sets TOL to be less than
10 * EPS * MAX( ||A|| , ||E|| ), then the tolerance is
F F
taken as 10 * EPS * MAX( ||A|| , ||E|| ), where EPS is the
F F
machine precision (see LAPACK Library Routine DLAMCH) and
A and E are matrices (as defined in METHOD).
Workspace
IWORK INTEGER array, dimension (m+2*MAX(n,m+1)+n),
where m = DP*MP and n = (DP-1)*MP + NP.
DWORK DOUBLE PRECISION array, dimension (LDWORK)
LDWORK The length of the array DWORK.
LDWORK >= m*n*n + 2*m*n + 2*n*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 taken during the
computations. This failure is not likely to occur.
The possible values are:
k, 1 <= k <= DK+1, the k-th diagonal submatrix had
not a full row rank;
DK+2, if incorrect dimensions of a full column
rank submatrix;
DK+3, if incorrect dimensions of a full row rank
submatrix.
Method
The computation of the right nullspace of the MP-by-NP polynomial
matrix P(s) of degree DP given by
DP-1 DP
P(s) = P(0) + P(1) * s + ... + P(DP-1) * s + P(DP) * s
is performed via the pencil s*E - A, associated with P(s), where
| I | | 0 -P(DP) |
| . | | I . . |
A = | . | and E = | . . . |. (1)
| . | | . 0 . |
| I | | I 0 -P(2) |
| P(0) | | I -P(1) |
The pencil s*E - A is transformed by unitary matrices Q and Z such
that
| sE(eps)-A(eps) | X | X |
|----------------|----------------|------------|
| 0 | sE(inf)-A(inf) | X |
Q'(s*E-A)Z = |=================================|============|.
| | |
| 0 | sE(r)-A(r) |
Since s*E(inf)-A(inf) and s*E(r)-A(r) have full column rank, the
minimal polynomial basis for the right nullspace of Q'(s*E-A)Z
(and consequently the basis for the right nullspace of s*E - A) is
completely determined by s*E(eps)-A(eps).
Let Veps(s) be a minimal polynomial basis for the right nullspace
of s*E(eps)-A(eps). Then
| Veps(s) |
V(s) = Z * |---------|
| 0 |
is a minimal polynomial basis for the right nullspace of s*E - A.
From the structure of s*E - A it can be shown that if V(s) is
partitioned as
| Vo(s) | (DP-1)*MP
V(s) = |------ |
| Ve(s) | NP
then the columns of Ve(s) form a minimal polynomial basis for the
right nullspace of P(s).
The vectors of Ve(s) are computed and stored in array NULLSP in
the following condensed form:
|| || | || | | || | |
|| U1,0 || U2,0 | U2,1 || U3,0 | U3,1 | U3,2 || U4,0 | ... |,
|| || | || | | || | |
where Ui,j is an NP-by-GAM(i) matrix which contains the i-th block
of columns of K(j), the j-th coefficient of the polynomial matrix
representation for the right nullspace
DK
K(s) = K(0) + K(1) * s + . . . + K(DK) * s .
The coefficients K(0), K(1), ..., K(DK) are NP-by-nk matrices
given by
K(0) = | U1,0 | U2,0 | U3,0 | . . . | U(DK+1,0) |
K(1) = | 0 | U2,1 | U3,1 | . . . | U(DK+1,1) |
K(2) = | 0 | 0 | U3,2 | . . . | U(DK+1,2) |
. . . . . . . . . .
K(DK) = | 0 | 0 | 0 | . . . | 0 | U(DK+1,DK)|.
Note that the degree of K(s) satisfies the inequality DK <=
DP * MIN(MP,NP) and that the dimension of K(s) satisfies the
inequality (NP-MP) <= nk <= NP.
References
[1] Beelen, Th.G.J.
New Algorithms for Computing the Kronecker structure of a
Pencil with Applications to Systems and Control Theory.
Ph.D.Thesis, Eindhoven University of Technology, 1987.
[2] Van Den Hurk, G.J.H.H.
New Algorithms for Solving Polynomial Matrix Problems.
Master's Thesis, Eindhoven University of Technology, 1987.
Numerical Aspects
The algorithm used by the routine involves the construction of a special block echelon form with pivots considered to be non-zero when they are larger than TOL. These pivots are then inverted in order to construct the columns of the kernel of the polynomial matrix. If TOL is chosen to be too small then these inversions may be sensitive whereas increasing TOL will make the inversions more robust but will affect the block echelon form (and hence the column degrees of the polynomial kernel). Furthermore, if the elements of the computed polynomial kernel are large relative to the polynomial matrix, then the user should consider trying several values of TOL.Further Comments
It also possible to compute a minimal polynomial basis for the right nullspace of a pencil, since a pencil is a polynomial matrix of degree 1. Thus for the pencil (s*E - A), the required input is P(1) = E and P(0) = -A. The routine can also be used to compute a minimal polynomial basis for the left nullspace of a polynomial matrix by simply transposing P(s).Example
Program Text
* MC03ND EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER DPMAX, MPMAX, NPMAX
* PARAMETER ( DPMAX = 5, MPMAX = 5, NPMAX = 5 )
PARAMETER ( DPMAX = 2, MPMAX = 5, NPMAX = 4 )
INTEGER LDP1, LDP2, LDNULL, LDKER1, LDKER2
PARAMETER ( LDP1 = MPMAX, LDP2 = NPMAX, LDNULL = NPMAX,
$ LDKER1 = NPMAX, LDKER2 = NPMAX )
INTEGER M, N
PARAMETER ( M = DPMAX*MPMAX, N = ( DPMAX-1 )*MPMAX+NPMAX )
INTEGER LIWORK, LDWORK
* PARAMETER ( LIWORK = 3*( N+M )+2,
PARAMETER ( LIWORK = M+2*MAX( N,M+1 )+N,
$ LDWORK = M*N**2+2*M*N+2*N**2 )
* .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER DK, DP, I, INFO, J, K, M1, MP, NK, NP
* .. Local Arrays ..
DOUBLE PRECISION DWORK(LDWORK), KER(LDKER1,LDKER2,M+1),
$ NULLSP(LDNULL,(M+1)*NPMAX), P(LDP1,LDP2,DPMAX+1)
INTEGER GAM(M+1), IWORK(LIWORK)
* .. External Subroutines ..
EXTERNAL MC03ND
* .. 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 = * ) MP, NP, DP, TOL
IF ( MP.LT.0 .OR. MP.GT.MPMAX ) THEN
WRITE ( NOUT, FMT = 99990 ) MP
ELSE IF ( NP.LT.0 .OR. NP.GT.NPMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) NP
ELSE IF ( DP.LE.0 .OR. DP.GT.DPMAX ) THEN
WRITE ( NOUT, FMT = 99992 ) DP
ELSE
DO 40 K = 1, DP + 1
DO 20 I = 1, MP
READ ( NIN, FMT = * ) ( P(I,J,K), J = 1,NP )
20 CONTINUE
40 CONTINUE
* Compute a minimal polynomial basis K(s) of the given P(s).
CALL MC03ND( MP, NP, DP, P, LDP1, LDP2, DK, GAM, NULLSP,
$ LDNULL, KER, LDKER1, LDKER2, TOL, IWORK, DWORK,
$ LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE IF ( DK.LT.0 ) THEN
WRITE ( NOUT, FMT = 99997 )
ELSE
NK = 0
M1 = 0
DO 60 I = 1, DK + 1
NK = NK + GAM(I)
M1 = M1 + GAM(I)*I
60 CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 80 I = 1, NP
WRITE ( NOUT, FMT = 99995 ) ( NULLSP(I,J), J = 1,M1 )
80 CONTINUE
WRITE ( NOUT, FMT = 99994 ) DK, ( I-1, I = 1,DK+1 )
DO 120 I = 1, NP
DO 100 J = 1, NK
WRITE ( NOUT, FMT = 99993 )
$ I, J, ( KER(I,J,K), K = 1,DK+1 )
100 CONTINUE
120 CONTINUE
END IF
END IF
STOP
*
99999 FORMAT (' MC03ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MC03ND = ',I2)
99997 FORMAT (' The polynomial matrix P(s) has no right nullspace')
99996 FORMAT (' The right nullspace vectors of P(s) are ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (/' The minimal polynomial basis K(s) (of degree ',I2,') ',
$ 'for the right nullspace is ',//' power of s ',
$ 20I8)
99993 FORMAT (/' element (',I2,',',I2,') is ',20(1X,F7.2))
99992 FORMAT (/' DP is out of range.',/' DP = ',I5)
99991 FORMAT (/' NP is out of range.',/' NP = ',I5)
99990 FORMAT (/' MP is out of range.',/' MP = ',I5)
END
Program Data
MC03ND EXAMPLE PROGRAM DATA 5 4 2 0.0 2.0 2.0 0.0 3.0 0.0 4.0 0.0 6.0 8.0 8.0 0.0 12.0 0.0 0.0 0.0 0.0 2.0 2.0 0.0 3.0 1.0 0.0 1.0 0.0 0.0 0.0 2.0 0.0 4.0 0.0 4.0 0.0 2.0 2.0 0.0 3.0 3.0 2.0 1.0 3.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0Program Results
MC03ND EXAMPLE PROGRAM RESULTS The right nullspace vectors of P(s) are 0.0000 0.0000 0.0000 -0.8321 0.0000 0.1538 0.0000 -1.0000 0.0000 0.5547 0.0000 0.2308 The minimal polynomial basis K(s) (of degree 1) for the right nullspace is power of s 0 1 element ( 1, 1) is 0.00 0.00 element ( 1, 2) is 0.00 0.00 element ( 2, 1) is -0.83 0.00 element ( 2, 2) is 0.00 0.15 element ( 3, 1) is 0.00 0.00 element ( 3, 2) is -1.00 0.00 element ( 4, 1) is 0.55 0.00 element ( 4, 2) is 0.00 0.23