Purpose
To compute for the descriptor system (A-lambda E,B,C)
the orthogonal transformation matrices Q and Z such that the
transformed system (Q'*A*Z-lambda Q'*E*Z, Q'*B, C*Z) is
in a SVD-like coordinate form with
( A11 A12 ) ( Er 0 )
Q'*A*Z = ( ) , Q'*E*Z = ( ) ,
( A21 A22 ) ( 0 0 )
where Er is an upper triangular invertible matrix.
Optionally, the A22 matrix can be further reduced to the form
( Ar X )
A22 = ( ) ,
( 0 0 )
with Ar an upper triangular invertible matrix, and X either a full
or a zero matrix.
The left and/or right orthogonal transformations performed
to reduce E and A22 can be optionally accumulated.
Specification
SUBROUTINE TG01FD( COMPQ, COMPZ, JOBA, L, N, M, P, A, LDA, E, LDE,
$ B, LDB, C, LDC, Q, LDQ, Z, LDZ, RANKE, RNKA22,
$ TOL, IWORK, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
CHARACTER COMPQ, COMPZ, JOBA
INTEGER INFO, L, LDA, LDB, LDC, LDE, LDQ, LDWORK,
$ LDZ, M, N, P, RANKE, RNKA22
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IWORK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
$ DWORK( * ), E( LDE, * ), Q( LDQ, * ),
$ Z( LDZ, * )
Arguments
Mode Parameters
COMPQ CHARACTER*1
= 'N': do not compute Q;
= 'I': Q is initialized to the unit matrix, and the
orthogonal matrix Q is returned;
= 'U': Q must contain an orthogonal matrix Q1 on entry,
and the product Q1*Q is returned.
COMPZ CHARACTER*1
= 'N': do not compute Z;
= 'I': Z is initialized to the unit matrix, and the
orthogonal matrix Z is returned;
= 'U': Z must contain an orthogonal matrix Z1 on entry,
and the product Z1*Z is returned.
JOBA CHARACTER*1
= 'N': do not reduce A22;
= 'R': reduce A22 to a SVD-like upper triangular form.
= 'T': reduce A22 to an upper trapezoidal form.
Input/Output Parameters
L (input) INTEGER
The number of rows of matrices A, B, and E. L >= 0.
N (input) INTEGER
The number of columns of matrices A, E, and C. N >= 0.
M (input) INTEGER
The number of columns of matrix B. M >= 0.
P (input) INTEGER
The number of rows of matrix C. P >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading L-by-N part of this array must
contain the state dynamics matrix A.
On exit, the leading L-by-N part of this array contains
the transformed matrix Q'*A*Z. If JOBA = 'T', this matrix
is in the form
( A11 * * )
Q'*A*Z = ( * Ar X ) ,
( * 0 0 )
where A11 is a RANKE-by-RANKE matrix and Ar is a
RNKA22-by-RNKA22 invertible upper triangular matrix.
If JOBA = 'R' then A has the above form with X = 0.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,L).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading L-by-N part of this array must
contain the descriptor matrix E.
On exit, the leading L-by-N part of this array contains
the transformed matrix Q'*E*Z.
( Er 0 )
Q'*E*Z = ( ) ,
( 0 0 )
where Er is a RANKE-by-RANKE upper triangular invertible
matrix.
LDE INTEGER
The leading dimension of array E. LDE >= MAX(1,L).
B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
On entry, the leading L-by-M part of this array must
contain the input/state matrix B.
On exit, the leading L-by-M part of this array contains
the transformed matrix Q'*B.
LDB INTEGER
The leading dimension of array B.
LDB >= MAX(1,L) if M > 0 or LDB >= 1 if M = 0.
C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
On entry, the leading P-by-N part of this array must
contain the state/output matrix C.
On exit, the leading P-by-N part of this array contains
the transformed matrix C*Z.
LDC INTEGER
The leading dimension of array C. LDC >= MAX(1,P).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,L)
If COMPQ = 'N': Q is not referenced.
If COMPQ = 'I': on entry, Q need not be set;
on exit, the leading L-by-L part of this
array contains the orthogonal matrix Q,
where Q' is the product of Householder
transformations which are applied to A,
E, and B on the left.
If COMPQ = 'U': on entry, the leading L-by-L part of this
array must contain an orthogonal matrix
Q1;
on exit, the leading L-by-L part of this
array contains the orthogonal matrix
Q1*Q.
LDQ INTEGER
The leading dimension of array Q.
LDQ >= 1, if COMPQ = 'N';
LDQ >= MAX(1,L), if COMPQ = 'U' or 'I'.
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
If COMPZ = 'N': Z is not referenced.
If COMPZ = 'I': on entry, Z need not be set;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix Z,
which is the product of Householder
transformations applied to A, E, and C
on the right.
If COMPZ = 'U': on entry, the leading N-by-N part of this
array must contain an orthogonal matrix
Z1;
on exit, the leading N-by-N part of this
array contains the orthogonal matrix
Z1*Z.
LDZ INTEGER
The leading dimension of array Z.
LDZ >= 1, if COMPZ = 'N';
LDZ >= MAX(1,N), if COMPZ = 'U' or 'I'.
RANKE (output) INTEGER
The estimated rank of matrix E, and thus also the order
of the invertible upper triangular submatrix Er.
RNKA22 (output) INTEGER
If JOBA = 'R' or 'T', then RNKA22 is the estimated rank of
matrix A22, and thus also the order of the invertible
upper triangular submatrix Ar.
If JOBA = 'N', then RNKA22 is not referenced.
Tolerances
TOL DOUBLE PRECISION
The tolerance to be used in determining the rank of E
and of A22. If the user sets TOL > 0, then the given
value of TOL is used as a lower bound for the
reciprocal condition numbers of leading submatrices
of R or R22 in the QR decompositions E * P = Q * R of E
or A22 * P22 = Q22 * R22 of A22.
A submatrix whose estimated condition number is less than
1/TOL is considered to be of full rank. If the user sets
TOL <= 0, then an implicitly computed, default tolerance,
defined by TOLDEF = L*N*EPS, is used instead, where
EPS is the machine precision (see LAPACK Library routine
DLAMCH). TOL < 1.
Workspace
IWORK INTEGER array, dimension (N)
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 >= MAX( 1, N+P, MIN(L,N)+MAX(3*N-1,M,L) ).
For optimal 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
The routine computes a truncated QR factorization with column
pivoting of E, in the form
( E11 E12 )
E * P = Q * ( )
( 0 E22 )
and finds the largest RANKE-by-RANKE leading submatrix E11 whose
estimated condition number is less than 1/TOL. RANKE defines thus
the rank of matrix E. Further E22, being negligible, is set to
zero, and an orthogonal matrix Y is determined such that
( E11 E12 ) = ( Er 0 ) * Y .
The overal transformation matrix Z results as Z = P * Y' and the
resulting transformed matrices Q'*A*Z and Q'*E*Z have the form
( Er 0 ) ( A11 A12 )
E <- Q'* E * Z = ( ) , A <- Q' * A * Z = ( ) ,
( 0 0 ) ( A21 A22 )
where Er is an upper triangular invertible matrix.
If JOBA = 'R' the same reduction is performed on A22 to obtain it
in the form
( Ar 0 )
A22 = ( ) ,
( 0 0 )
with Ar an upper triangular invertible matrix.
If JOBA = 'T' then A22 is row compressed using the QR
factorization with column pivoting to the form
( Ar X )
A22 = ( )
( 0 0 )
with Ar an upper triangular invertible matrix.
The transformations are also applied to the rest of system
matrices
B <- Q' * B, C <- C * Z.
Numerical Aspects
The algorithm is numerically backward stable and requires 0( L*L*N ) floating point operations.Further Comments
NoneExample
Program Text
* TG01FD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER LMAX, NMAX, MMAX, PMAX
PARAMETER ( LMAX = 20, NMAX = 20, MMAX = 20, PMAX = 20 )
INTEGER LDA, LDB, LDC, LDE, LDQ, LDZ
PARAMETER ( LDA = LMAX, LDB = LMAX, LDC = PMAX,
$ LDE = LMAX, LDQ = LMAX, LDZ = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( 1, PMAX,
$ MIN(LMAX,NMAX)+MAX( 3*NMAX, MMAX, LMAX ) ) )
* .. Local Scalars ..
CHARACTER*1 COMPQ, COMPZ, JOBA
INTEGER I, INFO, J, L, M, N, P, RANKE, RNKA22
DOUBLE PRECISION TOL
* .. Local Arrays ..
INTEGER IWORK(NMAX)
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ DWORK(LDWORK), E(LDE,NMAX), Q(LDQ,LMAX),
$ Z(LDZ,NMAX)
* .. External Subroutines ..
EXTERNAL TG01FD
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) L, N, M, P, TOL
COMPQ = 'I'
COMPZ = 'I'
JOBA = 'R'
IF ( L.LT.0 .OR. L.GT.LMAX ) THEN
WRITE ( NOUT, FMT = 99989 ) L
ELSE
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99988 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,L )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,L )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99987 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,L )
IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
WRITE ( NOUT, FMT = 99986 ) P
ELSE
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
* Find the transformed descriptor system
* (A-lambda E,B,C).
CALL TG01FD( COMPQ, COMPZ, JOBA, L, N, M, P, A, LDA,
$ E, LDE, B, LDB, C, LDC, Q, LDQ, Z, LDZ,
$ RANKE, RNKA22, TOL, IWORK, DWORK, LDWORK,
$ INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99994 ) RANKE, RNKA22
WRITE ( NOUT, FMT = 99997 )
DO 10 I = 1, L
WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
10 CONTINUE
WRITE ( NOUT, FMT = 99996 )
DO 20 I = 1, L
WRITE ( NOUT, FMT = 99995 ) ( E(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 30 I = 1, L
WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,M )
30 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 40 I = 1, P
WRITE ( NOUT, FMT = 99995 ) ( C(I,J), J = 1,N )
40 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 50 I = 1, L
WRITE ( NOUT, FMT = 99995 ) ( Q(I,J), J = 1,L )
50 CONTINUE
WRITE ( NOUT, FMT = 99990 )
DO 60 I = 1, N
WRITE ( NOUT, FMT = 99995 ) ( Z(I,J), J = 1,N )
60 CONTINUE
END IF
END IF
END IF
END IF
END IF
STOP
*
99999 FORMAT (' TG01FD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TG01FD = ',I2)
99997 FORMAT (/' The transformed state dynamics matrix Q''*A*Z is ')
99996 FORMAT (/' The transformed descriptor matrix Q''*E*Z is ')
99995 FORMAT (20(1X,F8.4))
99994 FORMAT (' Rank of matrix E =', I5/
$ ' Rank of matrix A22 =', I5)
99993 FORMAT (/' The transformed input/state matrix Q''*B is ')
99992 FORMAT (/' The transformed state/output matrix C*Z is ')
99991 FORMAT (/' The left transformation matrix Q is ')
99990 FORMAT (/' The right transformation matrix Z is ')
99989 FORMAT (/' L is out of range.',/' L = ',I5)
99988 FORMAT (/' N is out of range.',/' N = ',I5)
99987 FORMAT (/' M is out of range.',/' M = ',I5)
99986 FORMAT (/' P is out of range.',/' P = ',I5)
END
Program Data
TG01FD EXAMPLE PROGRAM DATA
4 4 2 2 0.0
-1 0 0 3
0 0 1 2
1 1 0 4
0 0 0 0
1 2 0 0
0 1 0 1
3 9 6 3
0 0 2 0
1 0
0 0
0 1
1 1
-1 0 1 0
0 1 -1 1
Program Results
TG01FD EXAMPLE PROGRAM RESULTS Rank of matrix E = 3 Rank of matrix A22 = 1 The transformed state dynamics matrix Q'*A*Z is 2.0278 0.1078 3.9062 -2.1571 -0.0980 0.2544 1.6053 -0.1269 0.2713 0.7760 -0.3692 -0.4853 0.0690 -0.5669 -2.1974 0.3086 The transformed descriptor matrix Q'*E*Z is 10.1587 5.8230 1.3021 0.0000 0.0000 -2.4684 -0.1896 0.0000 0.0000 0.0000 1.0338 0.0000 0.0000 0.0000 0.0000 0.0000 The transformed input/state matrix Q'*B is -0.2157 -0.9705 0.3015 0.9516 0.7595 0.0991 1.1339 0.3780 The transformed state/output matrix C*Z is 0.3651 -1.0000 -0.4472 -0.8165 -1.0954 1.0000 -0.8944 0.0000 The left transformation matrix Q is -0.2157 -0.5088 0.6109 0.5669 -0.1078 -0.2544 -0.7760 0.5669 -0.9705 0.1413 -0.0495 -0.1890 0.0000 0.8102 0.1486 0.5669 The right transformation matrix Z is -0.3651 0.0000 0.4472 0.8165 -0.9129 0.0000 0.0000 -0.4082 0.0000 -1.0000 0.0000 0.0000 -0.1826 0.0000 -0.8944 0.4082