Purpose
To compute the inverse (Ai,Bi,Ci,Di) of a given system (A,B,C,D).Specification
SUBROUTINE AB07ND( N, M, A, LDA, B, LDB, C, LDC, D, LDD, RCOND,
$ IWORK, DWORK, LDWORK, INFO )
C .. Scalar Arguments ..
DOUBLE PRECISION RCOND
INTEGER INFO, LDA, LDB, LDC, LDD, LDWORK, M, N
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
$ DWORK(*)
INTEGER IWORK(*)
Arguments
Input/Output Parameters
N (input) INTEGER
The order of the state matrix A. N >= 0.
M (input) INTEGER
The number of system inputs and outputs. M >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading N-by-N part of this array must
contain the state matrix A of the original system.
On exit, the leading N-by-N part of this array contains
the state matrix Ai of the inverse system.
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 of the original system.
On exit, the leading N-by-M part of this array contains
the input matrix Bi of the inverse system.
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 M-by-N part of this array must
contain the output matrix C of the original system.
On exit, the leading M-by-N part of this array contains
the output matrix Ci of the inverse system.
LDC INTEGER
The leading dimension of the array C. LDC >= MAX(1,M).
D (input/output) DOUBLE PRECISION array, dimension (LDD,M)
On entry, the leading M-by-M part of this array must
contain the feedthrough matrix D of the original system.
On exit, the leading M-by-M part of this array contains
the feedthrough matrix Di of the inverse system.
LDD INTEGER
The leading dimension of the array D. LDD >= MAX(1,M).
RCOND (output) DOUBLE PRECISION
The estimated reciprocal condition number of the
feedthrough matrix D of the original system.
Workspace
IWORK INTEGER array, dimension (2*M)
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0 or M+1, DWORK(1) returns the optimal
value of LDWORK.
LDWORK INTEGER
The length of the array DWORK. LDWORK >= MAX(1,4*M).
For good 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;
= i: the matrix D is exactly singular; the (i,i) diagonal
element is zero, i <= M; RCOND was set to zero;
= M+1: the matrix D is numerically singular, i.e., RCOND
is less than the relative machine precision, EPS
(see LAPACK Library routine DLAMCH). The
calculations have been completed, but the results
could be very inaccurate.
Method
The matrices of the inverse system are computed with the formulas:
-1 -1 -1 -1
Ai = A - B*D *C, Bi = -B*D , Ci = D *C, Di = D .
Numerical Aspects
The accuracy depends mainly on the condition number of the matrix D to be inverted. The estimated reciprocal condition number is returned in RCOND.Further Comments
NoneExample
Program Text
* AB07ND EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX, MMAX
PARAMETER ( NMAX = 20, MMAX = 20 )
INTEGER LDA, LDB, LDC, LDD
PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = MMAX,
$ LDD = MMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = 4*MMAX )
* .. Local Scalars ..
INTEGER I, INFO, J, M, N
DOUBLE PRECISION RCOND
* .. Local Arrays ..
INTEGER IWORK(2*MMAX)
DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
$ D(LDD,MMAX), DWORK(LDWORK)
* .. External Subroutines ..
EXTERNAL AB07ND
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read in the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, M
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 )
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99991 ) M
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,M )
READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,M )
* Find the inverse of the ssr (A,B,C,D).
CALL AB07ND( N, M, A, LDA, B, LDB, C, LDC, D, LDD, RCOND,
$ IWORK, DWORK, LDWORK, INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
WRITE ( NOUT, FMT = 99997 )
DO 20 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( A(I,J), J = 1,N )
20 CONTINUE
WRITE ( NOUT, FMT = 99995 )
DO 40 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,M )
40 CONTINUE
WRITE ( NOUT, FMT = 99994 )
DO 60 I = 1, M
WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,N )
60 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 80 I = 1, M
WRITE ( NOUT, FMT = 99996 ) ( D(I,J), J = 1,M )
80 CONTINUE
END IF
END IF
END IF
STOP
*
99999 FORMAT (' AB07ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from AB07ND = ',I2)
99997 FORMAT (' The state dynamics matrix of the inverse system is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' The input/state matrix of the inverse system is ')
99994 FORMAT (/' The state/output matrix of the inverse system is ')
99993 FORMAT (/' The feedthrough matrix of the inverse system is ')
99992 FORMAT (/' N is out of range.',/' N = ',I5)
99991 FORMAT (/' M is out of range.',/' M = ',I5)
END
Program Data
AB07ND EXAMPLE PROGRAM DATA 3 2 1.0 2.0 0.0 4.0 -1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 -1.0 0.0 0.0 1.0 4.0 0.0 0.0 1.0Program Results
AB07ND EXAMPLE PROGRAM RESULTS The state dynamics matrix of the inverse system is 1.0000 1.7500 0.2500 4.0000 -1.0000 -1.0000 0.0000 -0.2500 1.2500 The input/state matrix of the inverse system is -0.2500 0.0000 0.0000 -1.0000 -0.2500 0.0000 The state/output matrix of the inverse system is 0.0000 0.2500 -0.2500 0.0000 0.0000 1.0000 The feedthrough matrix of the inverse system is 0.2500 0.0000 0.0000 1.0000