1 SUBROUTINE TG04BX( IP, IZ, A, LDA, E, B, C, D, PR, PI, ZR, ZI, 2 $ GAIN, IWORK ) 3C 4C WARNING 5C 6C This routine is a modified version of TB04BX. It is intended 7C for the Octave Control Systems Package and supports Descriptor 8C State-Space models. TG04BX is *NOT* part of SLICOT and the 9C authors from NICONET e.V. are *NOT* responsible for it. 10C See file DESCRIPTION for the current maintainer of the Octave 11C control package. 12C 13C SLICOT RELEASE 5.0. 14C 15C Copyright (c) 2002-2009 NICONET e.V. 16C 17C This program is free software: you can redistribute it and/or 18C modify it under the terms of the GNU General Public License as 19C published by the Free Software Foundation, either version 2 of 20C the License, or (at your option) any later version. 21C 22C This program is distributed in the hope that it will be useful, 23C but WITHOUT ANY WARRANTY; without even the implied warranty of 24C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25C GNU General Public License for more details. 26C 27C You should have received a copy of the GNU General Public License 28C along with this program. If not, see 29C <http://www.gnu.org/licenses/>. 30C 31C PURPOSE 32C 33C To compute the gain of a single-input single-output linear system, 34C given its state-space representation (A,b,c,d), and its poles and 35C zeros. The matrix A is assumed to be in an upper Hessenberg form. 36C The gain is computed using the formula 37C 38C -1 IP IZ 39C g = (c*( S0*E - A ) *b + d)*Prod( S0 - Pi )/Prod( S0 - Zi ) , 40C i=1 i=1 (1) 41C 42C where Pi, i = 1 : IP, and Zj, j = 1 : IZ, are the poles and zeros, 43C respectively, and S0 is a real scalar different from all poles and 44C zeros. 45C 46C ARGUMENTS 47C 48C Input/Output Parameters 49C 50C IP (input) INTEGER 51C The number of the system poles. IP >= 0. 52C 53C IZ (input) INTEGER 54C The number of the system zeros. IZ >= 0. 55C 56C A (input/output) DOUBLE PRECISION array, dimension (LDA,IP) 57C On entry, the leading IP-by-IP part of this array must 58C contain the state dynamics matrix A in an upper Hessenberg 59C form. The elements below the second diagonal are not 60C referenced. 61C On exit, the leading IP-by-IP upper Hessenberg part of 62C this array contains the LU factorization of the matrix 63C A - S0*I, as computed by SLICOT Library routine MB02SD. 64C 65C LDA INTEGER 66C The leading dimension of array A. LDA >= max(1,IP). 67C 68C B (input/output) DOUBLE PRECISION array, dimension (IP) 69C On entry, this array must contain the system input 70C vector b. 71C On exit, this array contains the solution of the linear 72C system ( A - S0*I )x = b . 73C 74C C (input) DOUBLE PRECISION array, dimension (IP) 75C This array must contain the system output vector c. 76C 77C D (input) DOUBLE PRECISION 78C The variable must contain the system feedthrough scalar d. 79C 80C PR (input) DOUBLE PRECISION array, dimension (IP) 81C This array must contain the real parts of the system 82C poles. Pairs of complex conjugate poles must be stored in 83C consecutive memory locations. 84C 85C PI (input) DOUBLE PRECISION array, dimension (IP) 86C This array must contain the imaginary parts of the system 87C poles. 88C 89C ZR (input) DOUBLE PRECISION array, dimension (IZ) 90C This array must contain the real parts of the system 91C zeros. Pairs of complex conjugate zeros must be stored in 92C consecutive memory locations. 93C 94C ZI (input) DOUBLE PRECISION array, dimension (IZ) 95C This array must contain the imaginary parts of the system 96C zeros. 97C 98C GAIN (output) DOUBLE PRECISION 99C The gain of the linear system (A,b,c,d), given by (1). 100C 101C Workspace 102C 103C IWORK INTEGER array, dimension (IP) 104C On exit, it contains the pivot indices; for 1 <= i <= IP, 105C row i of the matrix A - S0*I was interchanged with 106C row IWORK(i). 107C 108C METHOD 109C 110C The routine implements the method presented in [1]. A suitable 111C value of S0 is chosen based on the system poles and zeros. 112C Then, the LU factorization of the upper Hessenberg, nonsingular 113C matrix A - S0*I is computed and used to solve the linear system 114C in (1). 115C 116C REFERENCES 117C 118C [1] Varga, A. and Sima, V. 119C Numerically Stable Algorithm for Transfer Function Matrix 120C Evaluation. 121C Int. J. Control, vol. 33, nr. 6, pp. 1123-1133, 1981. 122C 123C NUMERICAL ASPECTS 124C 125C The algorithm is numerically stable in practice and requires 126C O(IP*IP) floating point operations. 127C 128C CONTRIBUTORS 129C 130C V. Sima, Research Institute for Informatics, Bucharest, May 2002. 131C Partly based on the BIMASC Library routine GAIN by A. Varga. 132C 133C REVISIONS 134C 135C 2011-02-08 (Lukas Reichlin) Modifications for Descriptor Systems. 136C 137C KEYWORDS 138C 139C Eigenvalue, state-space representation, transfer function, zeros. 140C 141C ****************************************************************** 142C 143C .. Parameters .. 144 DOUBLE PRECISION ZERO, ONE, TWO, P1, ONEP1 145 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 146 $ P1 = 0.1D0, ONEP1 = 1.1D0 ) 147C .. Scalar Arguments .. 148 DOUBLE PRECISION D, GAIN 149 INTEGER IP, IZ, LDA 150C .. Array Arguments .. 151 DOUBLE PRECISION A(LDA,*), E(LDA,*), B(*), C(*), PI(*), PR(*), 152 $ ZI(*), ZR(*) 153 INTEGER IWORK(*) 154C .. Local Scalars .. 155 INTEGER I, J, INFO 156 DOUBLE PRECISION S0, S 157C .. External Functions .. 158 DOUBLE PRECISION DDOT 159 EXTERNAL DDOT 160C .. External Subroutines .. 161C EXTERNAL MB02RD, MB02SD 162 EXTERNAL DGETRF, DGETRS 163C .. Intrinsic Functions .. 164 INTRINSIC ABS, MAX 165C .. 166C .. Executable Statements .. 167C 168C For efficiency, the input scalar parameters are not checked. 169C 170C Quick return if possible. 171C 172C IF( IP.EQ.0 ) THEN 173C GAIN = ZERO 174C RETURN 175C END IF 176C 177C Compute a suitable value for S0 . 178C 179 S0 = ZERO 180C 181 DO 10 I = 1, IP 182 S = ABS( PR(I) ) 183 IF ( PI(I).NE.ZERO ) 184 $ S = S + ABS( PI(I) ) 185 S0 = MAX( S0, S ) 186 10 CONTINUE 187C 188 DO 20 I = 1, IZ 189 S = ABS( ZR(I) ) 190 IF ( ZI(I).NE.ZERO ) 191 $ S = S + ABS( ZI(I) ) 192 S0 = MAX( S0, S ) 193 20 CONTINUE 194C 195 S0 = TWO*S0 + P1 196 IF ( S0.LE.ONE ) 197 $ S0 = ONEP1 198C 199C Form A - S0*E . 200C 201 DO 30 J = 1, LDA 202 DO 25 I = 1, LDA 203 A(I,J) = A(I,J) - S0*E(I,J) 204 25 CONTINUE 205 30 CONTINUE 206C 207C Compute the LU factorization of the matrix A - S0*E 208C (guaranteed to be nonsingular). 209C 210C CALL MB02SD( IP, A, LDA, IWORK, INFO ) 211C CALL MB02SD( LDA, A, LDA, IWORK, INFO )oo 212 CALL DGETRF( LDA, LDA, A, LDA, IWORK, INFO) 213C 214C Solve the linear system (A - S0*E)*x = b . 215C 216C CALL MB02RD( 'No Transpose', IP, 1, A, LDA, IWORK, B, IP, INFO ) 217C CALL MB02RD( 'No Transpose', IP, 1, A, LDA, IWORK, B, LDA, INFO ) 218C CALL MB02RD( 'No Transpose', LDA, 1, A, LDA, IWORK, B, LDA, INFO ) 219 CALL DGETRS( 'No Transpose', LDA, 1, A, LDA, IWORK, B, LDA, INFO ) 220C -1 221C Compute c*(S0*E - A) *b + d . 222C 223C GAIN = D - DDOT( IP, C, 1, B, 1 ) 224 GAIN = D - DDOT( LDA, C, 1, B, 1 ) 225C 226C Multiply by the products in terms of poles and zeros in (1). 227C 228 I = 1 229C 230C WHILE ( I <= IP ) DO 231C 232 40 IF ( I.LE.IP ) THEN 233 IF ( PI(I).EQ.ZERO ) THEN 234 GAIN = GAIN*( S0 - PR(I) ) 235 I = I + 1 236 ELSE 237 GAIN = GAIN*( S0*( S0 - TWO*PR(I) ) + PR(I)**2 + PI(I)**2 ) 238 I = I + 2 239 END IF 240 GO TO 40 241 END IF 242C 243C END WHILE 40 244C 245 I = 1 246C 247C WHILE ( I <= IZ ) DO 248C 249 50 IF ( I.LE.IZ ) THEN 250 IF ( ZI(I).EQ.ZERO ) THEN 251 GAIN = GAIN/( S0 - ZR(I) ) 252 I = I + 1 253 ELSE 254 GAIN = GAIN/( S0*( S0 - TWO*ZR(I) ) + ZR(I)**2 + ZI(I)**2 ) 255 I = I + 2 256 END IF 257 GO TO 50 258 END IF 259C 260C END WHILE 50 261C 262 RETURN 263C *** Last line of TG04BX *** 264 END 265