1 SUBROUTINE SB06ND( N, M, KMAX, A, LDA, B, LDB, KSTAIR, U, LDU, F, 2 $ LDF, DWORK, INFO ) 3C 4C WARNING 5C 6C This file alters the SLICOT routine SB06ND. It is intended 7C for use from the Octave Control Systems Package and modifies 8C the original SLICOT implementation. This file itself is *NOT* 9C part of SLICOT and the authors from NICONET e.V. are *NOT* 10C responsible for it. See file DESCRIPTION for the current 11C maintainer of the Octave control package. 12C 13C In altered implementation the deprecated LAPACK routine DLATZM 14C is replaced by ODLTZM. 15C 16C 17C SLICOT RELEASE 5.0. 18C 19C Copyright (c) 2002-2010 NICONET e.V. 20C 21C This program is free software: you can redistribute it and/or 22C modify it under the terms of the GNU General Public License as 23C published by the Free Software Foundation, either version 2 of 24C the License, or (at your option) any later version. 25C 26C This program is distributed in the hope that it will be useful, 27C but WITHOUT ANY WARRANTY; without even the implied warranty of 28C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 29C GNU General Public License for more details. 30C 31C You should have received a copy of the GNU General Public License 32C along with this program. If not, see 33C <http://www.gnu.org/licenses/>. 34C 35C PURPOSE 36C 37C To construct the minimum norm feedback matrix F to perform 38C "deadbeat control" on a (A,B)-pair of a state-space model (which 39C must be preliminarily reduced to upper "staircase" form using 40C SLICOT Library routine AB01OD) such that the matrix R = A + BFU' 41C is nilpotent. 42C (The transformation matrix U reduces R to upper Schur form with 43C zero blocks on its diagonal (of dimension KSTAIR(i)) and 44C therefore contains bases for the i-th controllable subspaces, 45C where i = 1,...,KMAX). 46C 47C ARGUMENTS 48C 49C Input/Output Parameters 50C 51C N (input) INTEGER 52C The actual state dimension, i.e. the order of the 53C matrix A. N >= 0. 54C 55C M (input) INTEGER 56C The actual input dimension. M >= 0. 57C 58C KMAX (input) INTEGER 59C The number of "stairs" in the staircase form as produced 60C by SLICOT Library routine AB01OD. 0 <= KMAX <= N. 61C 62C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 63C On entry, the leading N-by-N part of this array must 64C contain the transformed state-space matrix of the 65C (A,B)-pair with triangular stairs, as produced by SLICOT 66C Library routine AB01OD (with option STAGES = 'A'). 67C On exit, the leading N-by-N part of this array contains 68C the matrix U'AU + U'BF. 69C 70C LDA INTEGER 71C The leading dimension of array A. LDA >= MAX(1,N). 72C 73C B (input/output) DOUBLE PRECISION array, dimension (LDB,M) 74C On entry, the leading N-by-M part of this array must 75C contain the transformed triangular input matrix of the 76C (A,B)-pair as produced by SLICOT Library routine AB01OD 77C (with option STAGES = 'A'). 78C On exit, the leading N-by-M part of this array contains 79C the matrix U'B. 80C 81C LDB INTEGER 82C The leading dimension of array B. LDB >= MAX(1,N). 83C 84C KSTAIR (input) INTEGER array, dimension (KMAX) 85C The leading KMAX elements of this array must contain the 86C dimensions of each "stair" as produced by SLICOT Library 87C routine AB01OD. 88C 89C U (input/output) DOUBLE PRECISION array, dimension (LDU,N) 90C On entry, the leading N-by-N part of this array must 91C contain either a transformation matrix (e.g. from a 92C previous call to other SLICOT routine) or be initialised 93C as the identity matrix. 94C On exit, the leading N-by-N part of this array contains 95C the product of the input matrix U and the state-space 96C transformation matrix which reduces A + BFU' to real 97C Schur form. 98C 99C LDU INTEGER 100C The leading dimension of array U. LDU >= MAX(1,N). 101C 102C F (output) DOUBLE PRECISION array, dimension (LDF,N) 103C The leading M-by-N part of this array contains the 104C deadbeat feedback matrix F. 105C 106C LDF INTEGER 107C The leading dimension of array F. LDF >= MAX(1,M). 108C 109C Workspace 110C 111C DWORK DOUBLE PRECISION array, dimension (2*N) 112C 113C Error Indicator 114C 115C INFO INTEGER 116C = 0: successful exit; 117C < 0: if INFO = -i, the i-th argument had an illegal 118C value. 119C 120C METHOD 121C 122C Starting from the (A,B)-pair in "staircase form" with "triangular" 123C stairs, dimensions KSTAIR(i+1) x KSTAIR(i), (described by the 124C vector KSTAIR): 125C 126C | B | A * . . . * | 127C | 1| 11 . . | 128C | | A A . . | 129C | | 21 22 . . | 130C | | . . . | 131C [ B | A ] = | | . . * | 132C | | . . | 133C | 0 | 0 | 134C | | A A | 135C | | r,r-1 rr | 136C 137C where the i-th diagonal block of A has dimension KSTAIR(i), for 138C i = 1,2,...,r, the feedback matrix F is constructed recursively in 139C r steps (where the number of "stairs" r is given by KMAX). In each 140C step a unitary state-space transformation U and a part of F are 141C updated in order to achieve the final form: 142C 143C | 0 A * . . . * | 144C | 12 . . | 145C | . . | 146C | 0 A . . | 147C | 23 . . | 148C | . . | 149C [ U'AU + U'BF ] = | . . * | . 150C | . . | 151C | | 152C | A | 153C | r-1,r| 154C | | 155C | 0 | 156C 157C 158C REFERENCES 159C 160C [1] Van Dooren, P. 161C Deadbeat control: a special inverse eigenvalue problem. 162C BIT, 24, pp. 681-699, 1984. 163C 164C NUMERICAL ASPECTS 165C 166C The algorithm requires O((N + M) * N**2) operations and is mixed 167C numerical stable (see [1]). 168C 169C CONTRIBUTORS 170C 171C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997. 172C Supersedes Release 2.0 routine SB06BD by M. Vanbegin, and 173C P. Van Dooren, Philips Research Laboratory, Brussels, Belgium. 174C 175C REVISIONS 176C 177C 1997, December 10; 2003, September 27. 178C 179C KEYWORDS 180C 181C Canonical form, deadbeat control, eigenvalue assignment, feedback 182C control, orthogonal transformation, real Schur form, staircase 183C form. 184C 185C ****************************************************************** 186C 187C .. Parameters .. 188 DOUBLE PRECISION ZERO, ONE 189 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 190C .. Scalar Arguments .. 191 INTEGER INFO, KMAX, LDA, LDB, LDF, LDU, M, N 192C .. Array Arguments .. 193 INTEGER KSTAIR(*) 194 DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*), F(LDF,*), U(LDU,*) 195C .. Local Scalars .. 196 INTEGER J, J0, JCUR, JKCUR, JMKCUR, KCUR, KK, KMIN, 197 $ KSTEP, MKCUR, NCONT 198C .. External Subroutines .. 199 EXTERNAL DCOPY, DGEMM, DLACPY, DLARFG, DLASET, ODLTZM, 200 $ DTRSM, XERBLA 201C .. Executable Statements .. 202C 203 INFO = 0 204C 205C Test the input scalar arguments. 206C 207 IF( N.LT.0 ) THEN 208 INFO = -1 209 ELSE IF( M.LT.0 ) THEN 210 INFO = -2 211 ELSE IF( KMAX.LT.0 .OR. KMAX.GT.N ) THEN 212 INFO = -3 213 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 214 INFO = -5 215 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 216 INFO = -7 217 ELSE IF( LDU.LT.MAX( 1, N ) ) THEN 218 INFO = -10 219 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 220 INFO = -12 221 ELSE 222 NCONT = 0 223C 224 DO 10 KK = 1, KMAX 225 NCONT = NCONT + KSTAIR(KK) 226 10 CONTINUE 227C 228 IF( NCONT.GT.N ) 229 $ INFO = -8 230 END IF 231C 232 IF ( INFO.NE.0 ) THEN 233C 234C Error return. 235C 236 CALL XERBLA( 'SB06ND', -INFO ) 237 RETURN 238 END IF 239C 240C Quick return if possible. 241C 242 IF ( N.EQ.0 .OR. M.EQ.0 ) 243 $ RETURN 244C 245 DO 120 KMIN = 1, KMAX 246 JCUR = NCONT 247 KSTEP = KMAX - KMIN 248C 249C Triangularize bottom part of A (if KSTEP > 0). 250C 251 DO 40 KK = KMAX, KMAX - KSTEP + 1, -1 252 KCUR = KSTAIR(KK) 253C 254C Construct Ukk and store in Fkk. 255C 256 DO 20 J = 1, KCUR 257 JMKCUR = JCUR - KCUR 258 CALL DCOPY( KCUR, A(JCUR,JMKCUR), LDA, F(1,JCUR), 1 ) 259 CALL DLARFG( KCUR+1, A(JCUR,JCUR), F(1,JCUR), 1, 260 $ DWORK(JCUR) ) 261 CALL DLASET( 'Full', 1, KCUR, ZERO, ZERO, A(JCUR,JMKCUR), 262 $ LDA ) 263C 264C Backmultiply A and U with Ukk. 265C 266 CALL ODLTZM( 'Right', JCUR-1, KCUR+1, F(1,JCUR), 1, 267 $ DWORK(JCUR), A(1,JCUR), A(1,JMKCUR), LDA, 268 $ DWORK ) 269C 270 CALL ODLTZM( 'Right', N, KCUR+1, F(1,JCUR), 1, 271 $ DWORK(JCUR), U(1,JCUR), U(1,JMKCUR), LDU, 272 $ DWORK(N+1) ) 273 JCUR = JCUR - 1 274 20 CONTINUE 275C 276 40 CONTINUE 277C 278C Eliminate diagonal block Aii by feedback Fi. 279C 280 KCUR = KSTAIR(KMIN) 281 J0 = JCUR - KCUR + 1 282 MKCUR = M - KCUR + 1 283C 284C Solve for Fi and add B x Fi to A. 285C 286 CALL DLACPY( 'Full', KCUR, KCUR, A(J0,J0), LDA, F(MKCUR,J0), 287 $ LDF ) 288 CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', KCUR, 289 $ KCUR, -ONE, B(J0,MKCUR), LDB, F(MKCUR,J0), LDF ) 290 IF ( J0.GT.1 ) 291 $ CALL DGEMM( 'No transpose', 'No transpose', J0-1, KCUR, 292 $ KCUR, ONE, B(1,MKCUR), LDB, F(MKCUR,J0), LDF, 293 $ ONE, A(1,J0), LDA ) 294 CALL DLASET( 'Full', KCUR, KCUR, ZERO, ZERO, A(J0,J0), LDA ) 295 CALL DLASET( 'Full', M-KCUR, KCUR, ZERO, ZERO, F(1,J0), LDF ) 296C 297 IF ( KSTEP.NE.0 ) THEN 298 JKCUR = NCONT 299C 300C Premultiply A with Ukk. 301C 302 DO 80 KK = KMAX, KMAX - KSTEP + 1, -1 303 KCUR = KSTAIR(KK) 304 JCUR = JKCUR - KCUR 305C 306 DO 60 J = 1, KCUR 307 CALL ODLTZM( 'Left', KCUR+1, N-JCUR+1, F(1,JKCUR), 1, 308 $ DWORK(JKCUR), A(JKCUR,JCUR), 309 $ A(JCUR,JCUR), LDA, DWORK(N+1) ) 310 JCUR = JCUR - 1 311 JKCUR = JKCUR - 1 312 60 CONTINUE 313C 314 80 CONTINUE 315C 316C Premultiply B with Ukk. 317C 318 JCUR = JCUR + KCUR 319 JKCUR = JCUR + KCUR 320C 321 DO 100 J = M, M - KCUR + 1, -1 322 CALL ODLTZM( 'Left', KCUR+1, M-J+1, F(1,JKCUR), 1, 323 $ DWORK(JKCUR), B(JKCUR,J), B(JCUR,J), LDB, 324 $ DWORK(N+1) ) 325 JCUR = JCUR - 1 326 JKCUR = JKCUR - 1 327 100 CONTINUE 328C 329 END IF 330 120 CONTINUE 331C 332 IF ( NCONT.NE.N ) 333 $ CALL DLASET( 'Full', M, N-NCONT, ZERO, ZERO, F(1,NCONT+1), 334 $ LDF ) 335C 336 RETURN 337C *** Last line of SB06ND *** 338 END 339