Purpose
To reduce the system state matrix A to an upper Hessenberg form by using an orthogonal similarity transformation A <-- U'*A*U and to apply the transformation to the matrices B and C: B <-- U'*B and C <-- C*U.Specification
SUBROUTINE TB01WX( COMPU, N, M, P, A, LDA, B, LDB, C, LDC, U, LDU,
$ DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER COMPU
INTEGER INFO, LDA, LDB, LDC, LDU, LDWORK, M, N, P
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), U(LDU,*)
Arguments
Mode Parameters
COMPU CHARACTER*1
= 'N': do not compute U;
= 'I': U is initialized to the unit matrix, and the
orthogonal matrix U is returned;
= 'U': U must contain an orthogonal matrix U1 on entry,
and the product U1*U is returned.
Input/Output Parameters
N (input) INTEGER
The order of the original state-space representation,
i.e., the order of the matrix A. N >= 0.
M (input) INTEGER
The number of system inputs, or of columns of B. M >= 0.
P (input) INTEGER
The number of system outputs, or of rows of C. P >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the original state dynamics matrix A.
On exit, the leading N-by-N part of this array contains
the matrix U' * A * U in Hessenberg form. The elements
below the first subdiagonal are set to zero.
LDA INTEGER
The leading dimension of the array A. LDA >= MAX(1,N).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
On entry, the leading N-by-M part of this array must
contain the input matrix B.
On exit, the leading N-by-M part of this array contains
the transformed input matrix U' * B.
LDB INTEGER
The leading dimension of the array B. LDB >= MAX(1,N).
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading P-by-N part of this array must
contain the output matrix C.
On exit, the leading P-by-N part of this array contains
the transformed output matrix C * U.
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,P).
U (input/output) DOUBLE PRECISION array, dimension (LDU,*)
On entry, if COMPU = 'U', the leading N-by-N part of this
array must contain the given matrix U1. Otherwise, this
array need not be set on input.
On exit, if COMPU <> 'N', the leading N-by-N part of this
array contains the orthogonal transformation matrix used
to reduce A to the Hessenberg form (U1*U if COMPU = 'U').
If COMPU = 'N', this array is not referenced.
LDU INTEGER
The leading dimension of the array U.
LDU >= 1, if COMPU = 'N';
LDU >= max(1,N), if COMPU <> 'N'.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= 1, and if N > 0,
LDWORK >= N - 1 + MAX(N,M,P).
For optimum performance LDWORK should be larger.
If LDWORK = -1, then a workspace query is assumed; the
routine only calculates the optimal size of the DWORK
array, returns this value as the first entry of the DWORK
array, and no error message related to LDWORK is issued by
XERBLA.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value.
Method
Matrix A is reduced to the Hessenberg form using an orthogonal similarity transformation A <- U'*A*U. Then, the transformation is applied to the matrices B and C: B <-- U'*B and C <-- C*U.Numerical Aspects
3 2
The algorithm requires about 5N /3 + N (M+P) floating point
3
operations, if COMPU = 'N'. Otherwise, 2N /3 additional operations
are needed.
Further Comments
NoneExample
Program Text
* TB01WX EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX, PMAX
PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER LDA, LDB, LDC, LDU
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
$ LDU = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = NMAX - 1 + MAX( NMAX, MMAX, PMAX ) )
* .. Local Scalars ..
CHARACTER COMPU
INTEGER I, INFO, J, M, N, P
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ DWORK(LDWORK), U(LDU,NMAX)
* .. External Subroutines ..
EXTERNAL TB01WX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M, P, COMPU
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99990 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1, N )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99988 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
* Find the transformed ssr for (A,B,C).
CALL TB01WX( COMPU, N, M, P, A, LDA, B, LDB, C, LDC, U,
$ LDU, DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 40 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
40 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 60 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,N )
60 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 70 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( U(I,J), J = 1,N )
70 CONTINUE
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TB01WX EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB01WX = ',I2)
99996 FORMAT (/' The transformed state dynamics matrix U''*A*U is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT ( ' (',F8.4,', ',F8.4,' )')
99993 FORMAT (/' The transformed input/state matrix U''*B is ')
99992 FORMAT (/' The transformed state/output matrix C*U is ')
99991 FORMAT (/' The similarity transformation matrix U is ')
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' M is out of range.',/' M = ',I5)
99988 FORMAT (/' P is out of range.',/' P = ',I5)
END
Program Data
TB01WX EXAMPLE PROGRAM DATA (Continuous system)
5 2 3 I
-0.04165 4.9200 -4.9200 0 0
-1.387944 -3.3300 0 0 0
0.5450 0 0 -0.5450 0
0 0 4.9200 -0.04165 4.9200
0 0 0 -1.387944 -3.3300
0 0
3.3300 0
0 0
0 0
0 3.3300
1 0 0 0 0
0 0 1 0 0
0 0 0 1 0
Program Results
TB01WX EXAMPLE PROGRAM RESULTS The transformed state dynamics matrix U'*A*U is -0.0416 -6.3778 1.4826 -1.9856 1.2630 1.4911 -2.8851 -0.4353 0.8984 -0.5714 0.0000 -2.1254 1.6804 -4.9686 -1.7731 0.0000 0.0000 2.1880 -3.3545 -2.6069 0.0000 0.0000 0.0000 0.7554 -2.1424 The transformed input/state matrix U'*B is 0.0000 0.0000 -3.0996 0.0000 -0.6488 0.0000 0.8689 1.7872 -0.5527 2.8098 The transformed state/output matrix C*U is 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.3655 -0.4962 0.6645 -0.4227 0.0000 0.0000 -0.8461 -0.4498 0.2861 The similarity transformation matrix U is 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.9308 -0.1948 0.2609 -0.1660 0.0000 0.3655 -0.4962 0.6645 -0.4227 0.0000 0.0000 -0.8461 -0.4498 0.2861 0.0000 0.0000 0.0000 0.5367 0.8438