1 SUBROUTINE AB08NX( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD, 2 $ NINFZ, INFZ, KRONL, MU, NU, NKROL, TOL, IWORK, 3 $ DWORK, LDWORK, INFO ) 4C 5C WARNING 6C 7C This file alters the SLICOT routine AB08NX. It is intended 8C for use from the Octave Control Systems Package and modifies 9C the original SLICOT implementation. This file itself is *NOT* 10C part of SLICOT and the authors from NICONET e.V. are *NOT* 11C responsible for it. See file DESCRIPTION for the current 12C maintainer of the Octave control package. 13C 14C In altered implementation the deprecated LAPACK routine DLATZM 15C is replaced by ODLTZM. 16C 17C 18C SLICOT RELEASE 5.0. 19C 20C Copyright (c) 2002-2010 NICONET e.V. 21C 22C This program is free software: you can redistribute it and/or 23C modify it under the terms of the GNU General Public License as 24C published by the Free Software Foundation, either version 2 of 25C the License, or (at your option) any later version. 26C 27C This program is distributed in the hope that it will be useful, 28C but WITHOUT ANY WARRANTY; without even the implied warranty of 29C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 30C GNU General Public License for more details. 31C 32C You should have received a copy of the GNU General Public License 33C along with this program. If not, see 34C <http://www.gnu.org/licenses/>. 35C 36C PURPOSE 37C 38C To extract from the (N+P)-by-(M+N) system 39C ( B A ) 40C ( D C ) 41C an (NU+MU)-by-(M+NU) "reduced" system 42C ( B' A') 43C ( D' C') 44C having the same transmission zeros but with D' of full row rank. 45C 46C ARGUMENTS 47C 48C Input/Output Parameters 49C 50C N (input) INTEGER 51C The number of state variables. N >= 0. 52C 53C M (input) INTEGER 54C The number of system inputs. M >= 0. 55C 56C P (input) INTEGER 57C The number of system outputs. P >= 0. 58C 59C RO (input/output) INTEGER 60C On entry, 61C = P for the original system; 62C = MAX(P-M, 0) for the pertransposed system. 63C On exit, RO contains the last computed rank. 64C 65C SIGMA (input/output) INTEGER 66C On entry, 67C = 0 for the original system; 68C = M for the pertransposed system. 69C On exit, SIGMA contains the last computed value sigma in 70C the algorithm. 71C 72C SVLMAX (input) DOUBLE PRECISION 73C During each reduction step, the rank-revealing QR 74C factorization of a matrix stops when the estimated minimum 75C singular value is smaller than TOL * MAX(SVLMAX,EMSV), 76C where EMSV is the estimated maximum singular value. 77C SVLMAX >= 0. 78C 79C ABCD (input/output) DOUBLE PRECISION array, dimension 80C (LDABCD,M+N) 81C On entry, the leading (N+P)-by-(M+N) part of this array 82C must contain the compound input matrix of the system. 83C On exit, the leading (NU+MU)-by-(M+NU) part of this array 84C contains the reduced compound input matrix of the system. 85C 86C LDABCD INTEGER 87C The leading dimension of array ABCD. 88C LDABCD >= MAX(1,N+P). 89C 90C NINFZ (input/output) INTEGER 91C On entry, the currently computed number of infinite zeros. 92C It should be initialized to zero on the first call. 93C NINFZ >= 0. 94C On exit, the number of infinite zeros. 95C 96C INFZ (input/output) INTEGER array, dimension (N) 97C On entry, INFZ(i) must contain the current number of 98C infinite zeros of degree i, where i = 1,2,...,N, found in 99C the previous call(s) of the routine. It should be 100C initialized to zero on the first call. 101C On exit, INFZ(i) contains the number of infinite zeros of 102C degree i, where i = 1,2,...,N. 103C 104C KRONL (input/output) INTEGER array, dimension (N+1) 105C On entry, this array must contain the currently computed 106C left Kronecker (row) indices found in the previous call(s) 107C of the routine. It should be initialized to zero on the 108C first call. 109C On exit, the leading NKROL elements of this array contain 110C the left Kronecker (row) indices. 111C 112C MU (output) INTEGER 113C The normal rank of the transfer function matrix of the 114C original system. 115C 116C NU (output) INTEGER 117C The dimension of the reduced system matrix and the number 118C of (finite) invariant zeros if D' is invertible. 119C 120C NKROL (output) INTEGER 121C The number of left Kronecker indices. 122C 123C Tolerances 124C 125C TOL DOUBLE PRECISION 126C A tolerance used in rank decisions to determine the 127C effective rank, which is defined as the order of the 128C largest leading (or trailing) triangular submatrix in the 129C QR (or RQ) factorization with column (or row) pivoting 130C whose estimated condition number is less than 1/TOL. 131C NOTE that when SVLMAX > 0, the estimated ranks could be 132C less than those defined above (see SVLMAX). 133C 134C Workspace 135C 136C IWORK INTEGER array, dimension (MAX(M,P)) 137C 138C DWORK DOUBLE PRECISION array, dimension (LDWORK) 139C On exit, if INFO = 0, DWORK(1) returns the optimal value 140C of LDWORK. 141C 142C LDWORK INTEGER 143C The length of the array DWORK. 144C LDWORK >= MAX( 1, MIN(P,M) + MAX(3*M-1,N), 145C MIN(P,N) + MAX(3*P-1,N+P,N+M) ). 146C For optimum performance LDWORK should be larger. 147C 148C If LDWORK = -1, then a workspace query is assumed; 149C the routine only calculates the optimal size of the 150C DWORK array, returns this value as the first entry of 151C the DWORK array, and no error message related to LDWORK 152C is issued by XERBLA. 153C 154C Error Indicator 155C 156C INFO INTEGER 157C = 0: successful exit; 158C < 0: if INFO = -i, the i-th argument had an illegal 159C value. 160C 161C REFERENCES 162C 163C [1] Svaricek, F. 164C Computation of the Structural Invariants of Linear 165C Multivariable Systems with an Extended Version of 166C the Program ZEROS. 167C System & Control Letters, 6, pp. 261-266, 1985. 168C 169C [2] Emami-Naeini, A. and Van Dooren, P. 170C Computation of Zeros of Linear Multivariable Systems. 171C Automatica, 18, pp. 415-430, 1982. 172C 173C NUMERICAL ASPECTS 174C 175C The algorithm is backward stable. 176C 177C CONTRIBUTOR 178C 179C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996. 180C Supersedes Release 2.0 routine AB08BZ by F. Svaricek. 181C 182C REVISIONS 183C 184C V. Sima, Oct. 1997; Feb. 1998, Jan. 2009, Apr. 2009. 185C A. Varga, May 1999; May 2001. 186C 187C KEYWORDS 188C 189C Generalized eigenvalue problem, Kronecker indices, multivariable 190C system, orthogonal transformation, structural invariant. 191C 192C ****************************************************************** 193C 194C .. Parameters .. 195 DOUBLE PRECISION ZERO 196 PARAMETER ( ZERO = 0.0D0 ) 197C .. Scalar Arguments .. 198 INTEGER INFO, LDABCD, LDWORK, M, MU, N, NINFZ, NKROL, 199 $ NU, P, RO, SIGMA 200 DOUBLE PRECISION SVLMAX, TOL 201C .. Array Arguments .. 202 INTEGER INFZ(*), IWORK(*), KRONL(*) 203 DOUBLE PRECISION ABCD(LDABCD,*), DWORK(*) 204C .. Local Scalars .. 205 LOGICAL LQUERY 206 INTEGER I1, IK, IROW, ITAU, IZ, JWORK, MM1, MNTAU, MNU, 207 $ MPM, NB, NP, RANK, RO1, TAU, WRKOPT 208 DOUBLE PRECISION T 209C .. Local Arrays .. 210 DOUBLE PRECISION SVAL(3) 211C .. External Functions .. 212 INTEGER ILAENV 213 EXTERNAL ILAENV 214C .. External Subroutines .. 215 EXTERNAL DLAPMT, DLARFG, DLASET, ODLTZM, DORMQR, DORMRQ, 216 $ MB03OY, MB03PY, XERBLA 217C .. Intrinsic Functions .. 218 INTRINSIC INT, MAX, MIN 219C .. Executable Statements .. 220C 221 NP = N + P 222 MPM = MIN( P, M ) 223 INFO = 0 224 LQUERY = ( LDWORK.EQ.-1 ) 225C 226C Test the input scalar arguments. 227C 228 IF( N.LT.0 ) THEN 229 INFO = -1 230 ELSE IF( M.LT.0 ) THEN 231 INFO = -2 232 ELSE IF( P.LT.0 ) THEN 233 INFO = -3 234 ELSE IF( RO.NE.P .AND. RO.NE.MAX( P-M, 0 ) ) THEN 235 INFO = -4 236 ELSE IF( SIGMA.NE.0 .AND. SIGMA.NE.M ) THEN 237 INFO = -5 238 ELSE IF( SVLMAX.LT.ZERO ) THEN 239 INFO = -6 240 ELSE IF( LDABCD.LT.MAX( 1, NP ) ) THEN 241 INFO = -8 242 ELSE IF( NINFZ.LT.0 ) THEN 243 INFO = -9 244 ELSE 245 JWORK = MAX( 1, MPM + MAX( 3*M - 1, N ), 246 $ MIN( P, N ) + MAX( 3*P - 1, NP, N+M ) ) 247 IF( LQUERY ) THEN 248 IF( M.GT.0 ) THEN 249 NB = MIN( 64, ILAENV( 1, 'DORMQR', 'LT', P, N, MPM, 250 $ -1 ) ) 251 WRKOPT = MAX( JWORK, MPM + MAX( 1, N )*NB ) 252 ELSE 253 WRKOPT = JWORK 254 END IF 255 NB = MIN( 64, ILAENV( 1, 'DORMRQ', 'RT', NP, N, MIN( P, N ), 256 $ -1 ) ) 257 WRKOPT = MAX( WRKOPT, MIN( P, N ) + MAX( 1, NP )*NB ) 258 NB = MIN( 64, ILAENV( 1, 'DORMRQ', 'LN', N, M+N, 259 $ MIN( P, N ), -1 ) ) 260 WRKOPT = MAX( WRKOPT, MIN( P, N ) + MAX( 1, M+N )*NB ) 261 ELSE IF( LDWORK.LT.JWORK ) THEN 262 INFO = -18 263 END IF 264 END IF 265C 266 IF ( INFO.NE.0 ) THEN 267C 268C Error return. 269C 270 CALL XERBLA( 'AB08NX', -INFO ) 271 RETURN 272 ELSE IF( LQUERY ) THEN 273 DWORK(1) = WRKOPT 274 RETURN 275 END IF 276C 277 MU = P 278 NU = N 279C 280 IZ = 0 281 IK = 1 282 MM1 = M + 1 283 ITAU = 1 284 NKROL = 0 285 WRKOPT = 1 286C 287C Main reduction loop: 288C 289C M NU M NU 290C NU [ B A ] NU [ B A ] 291C MU [ D C ] --> SIGMA [ RD C1 ] (SIGMA = rank(D) = 292C TAU [ 0 C2 ] row size of RD) 293C 294C M NU-RO RO 295C NU-RO [ B1 A11 A12 ] 296C --> RO [ B2 A21 A22 ] (RO = rank(C2) = 297C SIGMA [ RD C11 C12 ] col size of LC) 298C TAU [ 0 0 LC ] 299C 300C M NU-RO 301C NU-RO [ B1 A11 ] NU := NU - RO 302C [----------] MU := RO + SIGMA 303C --> RO [ B2 A21 ] D := [B2;RD] 304C SIGMA [ RD C11 ] C := [A21;C11] 305C 306 20 IF ( MU.EQ.0 ) 307 $ GO TO 80 308C 309C (Note: Comments in the code beginning "Workspace:" describe the 310C minimal amount of real workspace needed at that point in the 311C code, as well as the preferred amount for good performance.) 312C 313 RO1 = RO 314 MNU = M + NU 315 IF ( M.GT.0 ) THEN 316 IF ( SIGMA.NE.0 ) THEN 317 IROW = NU + 1 318C 319C Compress rows of D. First exploit triangular shape. 320C Workspace: need M+N-1. 321C 322 DO 40 I1 = 1, SIGMA 323 CALL DLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, T ) 324 CALL ODLTZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, T, 325 $ ABCD(IROW,I1+1), ABCD(IROW+1,I1+1), LDABCD, 326 $ DWORK ) 327 IROW = IROW + 1 328 40 CONTINUE 329 CALL DLASET( 'Lower', RO+SIGMA-1, SIGMA, ZERO, ZERO, 330 $ ABCD(NU+2,1), LDABCD ) 331 END IF 332C 333C Continue with Householder with column pivoting. 334C 335C The rank of D is the number of (estimated) singular values 336C that are greater than TOL * MAX(SVLMAX,EMSV). This number 337C includes the singular values of the first SIGMA columns. 338C Integer workspace: need M; 339C Workspace: need min(RO1,M) + 3*M - 1. RO1 <= P. 340C 341 IF ( SIGMA.LT.M ) THEN 342 JWORK = ITAU + MIN( RO1, M ) 343 I1 = SIGMA + 1 344 IROW = NU + I1 345 CALL MB03OY( RO1, M-SIGMA, ABCD(IROW,I1), LDABCD, TOL, 346 $ SVLMAX, RANK, SVAL, IWORK, DWORK(ITAU), 347 $ DWORK(JWORK), INFO ) 348 WRKOPT = MAX( WRKOPT, JWORK + 3*M - 2 ) 349C 350C Apply the column permutations to matrices B and part of D. 351C 352 CALL DLAPMT( .TRUE., NU+SIGMA, M-SIGMA, ABCD(1,I1), LDABCD, 353 $ IWORK ) 354C 355 IF ( RANK.GT.0 ) THEN 356C 357C Apply the Householder transformations to the submatrix C. 358C Workspace: need min(RO1,M) + NU; 359C prefer min(RO1,M) + NU*NB. 360C 361 CALL DORMQR( 'Left', 'Transpose', RO1, NU, RANK, 362 $ ABCD(IROW,I1), LDABCD, DWORK(ITAU), 363 $ ABCD(IROW,MM1), LDABCD, DWORK(JWORK), 364 $ LDWORK-JWORK+1, INFO ) 365 WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 ) 366 IF ( RO1.GT.1 ) 367 $ CALL DLASET( 'Lower', RO1-1, MIN( RO1-1, RANK ), ZERO, 368 $ ZERO, ABCD(IROW+1,I1), LDABCD ) 369 RO1 = RO1 - RANK 370 END IF 371 END IF 372 END IF 373C 374 TAU = RO1 375 SIGMA = MU - TAU 376C 377C Determination of the orders of the infinite zeros. 378C 379 IF ( IZ.GT.0 ) THEN 380 INFZ(IZ) = INFZ(IZ) + RO - TAU 381 NINFZ = NINFZ + IZ*( RO - TAU ) 382 END IF 383 IF ( RO1.EQ.0 ) 384 $ GO TO 80 385 IZ = IZ + 1 386C 387 IF ( NU.LE.0 ) THEN 388 MU = SIGMA 389 NU = 0 390 RO = 0 391 ELSE 392C 393C Compress the columns of C2 using RQ factorization with row 394C pivoting, P * C2 = R * Q. 395C 396 I1 = NU + SIGMA + 1 397 MNTAU = MIN( TAU, NU ) 398 JWORK = ITAU + MNTAU 399C 400C The rank of C2 is the number of (estimated) singular values 401C greater than TOL * MAX(SVLMAX,EMSV). 402C Integer Workspace: need TAU; 403C Workspace: need min(TAU,NU) + 3*TAU - 1. 404C 405 CALL MB03PY( TAU, NU, ABCD(I1,MM1), LDABCD, TOL, SVLMAX, RANK, 406 $ SVAL, IWORK, DWORK(ITAU), DWORK(JWORK), INFO ) 407 WRKOPT = MAX( WRKOPT, JWORK + 3*TAU - 1 ) 408 IF ( RANK.GT.0 ) THEN 409 IROW = I1 + TAU - RANK 410C 411C Apply Q' to the first NU columns of [A; C1] from the right. 412C Workspace: need min(TAU,NU) + NU + SIGMA; SIGMA <= P; 413C prefer min(TAU,NU) + (NU + SIGMA)*NB. 414C 415 CALL DORMRQ( 'Right', 'Transpose', I1-1, NU, RANK, 416 $ ABCD(IROW,MM1), LDABCD, DWORK(MNTAU-RANK+1), 417 $ ABCD(1,MM1), LDABCD, DWORK(JWORK), 418 $ LDWORK-JWORK+1, INFO ) 419 WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 ) 420C 421C Apply Q to the first NU rows and M + NU columns of [ B A ] 422C from the left. 423C Workspace: need min(TAU,NU) + M + NU; 424C prefer min(TAU,NU) + (M + NU)*NB. 425C 426 CALL DORMRQ( 'Left', 'NoTranspose', NU, MNU, RANK, 427 $ ABCD(IROW,MM1), LDABCD, DWORK(MNTAU-RANK+1), 428 $ ABCD, LDABCD, DWORK(JWORK), LDWORK-JWORK+1, 429 $ INFO ) 430 WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 ) 431C 432 CALL DLASET( 'Full', RANK, NU-RANK, ZERO, ZERO, 433 $ ABCD(IROW,MM1), LDABCD ) 434 IF ( RANK.GT.1 ) 435 $ CALL DLASET( 'Lower', RANK-1, RANK-1, ZERO, ZERO, 436 $ ABCD(IROW+1,MM1+NU-RANK), LDABCD ) 437 END IF 438C 439 RO = RANK 440 END IF 441C 442C Determine the left Kronecker indices (row indices). 443C 444 KRONL(IK) = KRONL(IK) + TAU - RO 445 NKROL = NKROL + KRONL(IK) 446 IK = IK + 1 447C 448C C and D are updated to [A21 ; C11] and [B2 ; RD]. 449C 450 NU = NU - RO 451 MU = SIGMA + RO 452 IF ( RO.NE.0 ) 453 $ GO TO 20 454C 455 80 CONTINUE 456 DWORK(1) = WRKOPT 457 RETURN 458C *** Last line of AB08NX *** 459 END 460