1 SUBROUTINE TB01ND( JOBU, UPLO, N, P, A, LDA, C, LDC, U, LDU, 2 $ DWORK, INFO ) 3C 4C WARNING 5C 6C This file alters the SLICOT routine TB01ND. 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 reduce the pair (A,C) to lower or upper observer Hessenberg 38C form using (and optionally accumulating) unitary state-space 39C transformations. 40C 41C ARGUMENTS 42C 43C Mode Parameters 44C 45C JOBU CHARACTER*1 46C Indicates whether the user wishes to accumulate in a 47C matrix U the unitary state-space transformations for 48C reducing the system, as follows: 49C = 'N': Do not form U; 50C = 'I': U is initialized to the unit matrix and the 51C unitary transformation matrix U is returned; 52C = 'U': The given matrix U is updated by the unitary 53C transformations used in the reduction. 54C 55C UPLO CHARACTER*1 56C Indicates whether the user wishes the pair (A,C) to be 57C reduced to upper or lower observer Hessenberg form as 58C follows: 59C = 'U': Upper observer Hessenberg form; 60C = 'L': Lower observer Hessenberg form. 61C 62C Input/Output Parameters 63C 64C N (input) INTEGER 65C The actual state dimension, i.e. the order of the 66C matrix A. N >= 0. 67C 68C P (input) INTEGER 69C The actual output dimension, i.e. the number of rows of 70C the matrix C. P >= 0. 71C 72C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 73C On entry, the leading N-by-N part of this array must 74C contain the state transition matrix A to be transformed. 75C On exit, the leading N-by-N part of this array contains 76C the transformed state transition matrix U' * A * U. 77C The annihilated elements are set to zero. 78C 79C LDA INTEGER 80C The leading dimension of array A. LDA >= MAX(1,N). 81C 82C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 83C On entry, the leading P-by-N part of this array must 84C contain the output matrix C to be transformed. 85C On exit, the leading P-by-N part of this array contains 86C the transformed output matrix C * U. 87C The annihilated elements are set to zero. 88C 89C LDC INTEGER 90C The leading dimension of array C. LDC >= MAX(1,P). 91C 92C U (input/output) DOUBLE PRECISION array, dimension (LDU,*) 93C On entry, if JOBU = 'U', then the leading N-by-N part of 94C this array must contain a given matrix U (e.g. from a 95C previous call to another SLICOT routine), and on exit, the 96C leading N-by-N part of this array contains the product of 97C the input matrix U and the state-space transformation 98C matrix which reduces the given pair to observer Hessenberg 99C form. 100C On exit, if JOBU = 'I', then the leading N-by-N part of 101C this array contains the matrix of accumulated unitary 102C similarity transformations which reduces the given pair 103C to observer Hessenberg form. 104C If JOBU = 'N', the array U is not referenced and can be 105C supplied as a dummy array (i.e. set parameter LDU = 1 and 106C declare this array to be U(1,1) in the calling program). 107C 108C LDU INTEGER 109C The leading dimension of array U. If JOBU = 'U' or 110C JOBU = 'I', LDU >= MAX(1,N); if JOBU = 'N', LDU >= 1. 111C 112C Workspace 113C 114C DWORK DOUBLE PRECISION array, dimension (MAX(N,P-1)) 115C 116C Error Indicator 117C 118C INFO INTEGER 119C = 0: successful exit; 120C < 0: if INFO = -i, the i-th argument had an illegal 121C value. 122C 123C METHOD 124C 125C The routine computes a unitary state-space transformation U, which 126C reduces the pair (A,C) to one of the following observer Hessenberg 127C forms: 128C 129C N 130C |* . . . . . . *| 131C |. .| 132C |. .| 133C |. .| N 134C |* .| 135C |U'AU| | . .| 136C |----| = | . .| 137C |CU | | * . . . *| 138C ------------------- 139C | * . . *| 140C | . .| P 141C | . .| 142C | *| 143C 144C if UPLO = 'U', or 145C 146C N 147C |* | 148C |. . | 149C |. . | P 150C |* . . * | 151C |CU | ------------------- 152C |----| = |* . . . * | 153C |U'AU| |. . | 154C |. . | 155C |. *| 156C |. .| N 157C |. .| 158C |. .| 159C |* . . . . . . *| 160C 161C if UPLO = 'L'. 162C 163C If P >= N, then the matrix CU is trapezoidal and U'AU is full. 164C 165C REFERENCES 166C 167C [1] Van Dooren, P. and Verhaegen, M.H.G. 168C On the use of unitary state-space transformations. 169C In : Contemporary Mathematics on Linear Algebra and its Role 170C in Systems Theory, 47, AMS, Providence, 1985. 171C 172C NUMERICAL ASPECTS 173C 174C The algorithm requires O((N + P) x N**2) operations and is 175C backward stable (see [1]). 176C 177C CONTRIBUTORS 178C 179C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Dec. 1996. 180C Supersedes Release 2.0 routine TB01BD by M. Vanbegin, and 181C P. Van Dooren, Philips Research Laboratory, Brussels, Belgium. 182C 183C REVISIONS 184C 185C February 1997. 186C 187C KEYWORDS 188C 189C Controllability, observer Hessenberg form, orthogonal 190C transformation, unitary transformation. 191C 192C ****************************************************************** 193C 194C .. Parameters .. 195 DOUBLE PRECISION ZERO, ONE 196 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 197C .. Scalar Arguments .. 198 INTEGER INFO, LDA, LDC, LDU, N, P 199 CHARACTER JOBU, UPLO 200C .. Array Arguments .. 201 DOUBLE PRECISION A(LDA,*), C(LDC,*), DWORK(*), U(LDU,*) 202C .. Local Scalars .. 203 LOGICAL LJOBA, LJOBI, LUPLO 204 INTEGER II, J, N1, NJ, P1, PAR1, PAR2, PAR3, PAR4, PAR5, 205 $ PAR6 206 DOUBLE PRECISION DZ 207C .. External Functions .. 208 LOGICAL LSAME 209 EXTERNAL LSAME 210C .. External Subroutines .. 211 EXTERNAL DLARFG, DLASET, ODLTZM, XERBLA 212C .. Intrinsic Functions .. 213 INTRINSIC MAX, MIN 214C .. Executable Statements .. 215C 216 INFO = 0 217 LUPLO = LSAME( UPLO, 'U' ) 218 LJOBI = LSAME( JOBU, 'I' ) 219 LJOBA = LJOBI.OR.LSAME( JOBU, 'U' ) 220C 221C Test the input scalar arguments. 222C 223 IF( .NOT.LJOBA .AND. .NOT.LSAME( JOBU, 'N' ) ) THEN 224 INFO = -1 225 ELSE IF( .NOT.LUPLO .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 226 INFO = -2 227 ELSE IF( N.LT.0 ) THEN 228 INFO = -3 229 ELSE IF( P.LT.0 ) THEN 230 INFO = -4 231 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 232 INFO = -6 233 ELSE IF( LDC.LT.MAX( 1, P ) ) THEN 234 INFO = -8 235 ELSE IF( .NOT.LJOBA .AND. LDU.LT.1 .OR. 236 $ LJOBA .AND. LDU.LT.MAX( 1, N ) ) THEN 237 INFO = -10 238 END IF 239C 240 IF ( INFO.NE.0 ) THEN 241C 242C Error return 243C 244 CALL XERBLA( 'TB01ND', -INFO ) 245 RETURN 246 END IF 247C 248C Quick return if possible. 249C 250 IF ( N.EQ.0 .OR. P.EQ.0 ) 251 $ RETURN 252C 253 P1 = P + 1 254 N1 = N - 1 255C 256 IF ( LJOBI ) THEN 257C 258C Initialize U to the identity matrix. 259C 260 CALL DLASET( 'Full', N, N, ZERO, ONE, U, LDU ) 261 END IF 262C 263C Perform transformations involving both C and A. 264C 265 DO 20 J = 1, MIN( P, N1 ) 266 NJ = N - J 267 IF ( LUPLO ) THEN 268 PAR1 = P - J + 1 269 PAR2 = NJ + 1 270 PAR3 = 1 271 PAR4 = P - J 272 PAR5 = NJ 273 ELSE 274 PAR1 = J 275 PAR2 = J 276 PAR3 = J + 1 277 PAR4 = P 278 PAR5 = N 279 END IF 280C 281 CALL DLARFG( NJ+1, C(PAR1,PAR2), C(PAR1,PAR3), LDC, DZ ) 282C 283C Update A. 284C 285 CALL ODLTZM( 'Left', NJ+1, N, C(PAR1,PAR3), LDC, DZ, A(PAR2,1), 286 $ A(PAR3,1), LDA, DWORK ) 287 CALL ODLTZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ, 288 $ A(1,PAR2), A(1,PAR3), LDA, DWORK ) 289C 290 IF ( LJOBA ) THEN 291C 292C Update U. 293C 294 CALL ODLTZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ, 295 $ U(1,PAR2), U(1,PAR3), LDU, DWORK ) 296 END IF 297C 298 IF ( J.NE.P ) THEN 299C 300C Update C. 301C 302 CALL ODLTZM( 'Right', PAR4-PAR3+1, NJ+1, C(PAR1,PAR3), LDC, 303 $ DZ, C(PAR3,PAR2), C(PAR3,PAR3), LDC, DWORK ) 304 END IF 305C 306 DO 10 II = PAR3, PAR5 307 C(PAR1,II) = ZERO 308 10 CONTINUE 309C 310 20 CONTINUE 311C 312 DO 40 J = P1, N1 313C 314C Perform next transformations only involving A. 315C 316 NJ = N - J 317 IF ( LUPLO ) THEN 318 PAR1 = N + P1 - J 319 PAR2 = NJ + 1 320 PAR3 = 1 321 PAR4 = NJ 322 PAR5 = 1 323 PAR6 = N + P - J 324 ELSE 325 PAR1 = J - P 326 PAR2 = J 327 PAR3 = J + 1 328 PAR4 = N 329 PAR5 = J - P + 1 330 PAR6 = N 331 END IF 332C 333 IF ( NJ.GT.0 ) THEN 334C 335 CALL DLARFG( NJ+1, A(PAR1,PAR2), A(PAR1,PAR3), LDA, DZ ) 336C 337C Update A. 338C 339 CALL ODLTZM( 'Left', NJ+1, N, A(PAR1,PAR3), LDA, DZ, 340 $ A(PAR2,1), A(PAR3,1), LDA, DWORK ) 341 CALL ODLTZM( 'Right', PAR6-PAR5+1, NJ+1, A(PAR1,PAR3), LDA, 342 $ DZ, A(PAR5,PAR2), A(PAR5,PAR3), LDA, DWORK ) 343C 344 IF ( LJOBA ) THEN 345C 346C Update U. 347C 348 CALL ODLTZM( 'Right', N, NJ+1, A(PAR1,PAR3), LDA, DZ, 349 $ U(1,PAR2), U(1,PAR3), LDU, DWORK ) 350 END IF 351C 352 DO 30 II = PAR3, PAR4 353 A(PAR1,II) = ZERO 354 30 CONTINUE 355C 356 END IF 357C 358 40 CONTINUE 359C 360 RETURN 361C *** Last line of TB01ND *** 362 END 363