1 SUBROUTINE SB01BY( N, M, S, P, A, B, F, TOL, DWORK, INFO ) 2C 3C WARNING 4C 5C This file alters the SLICOT routine SB01BY. It is intended 6C for use from the Octave Control Systems Package and modifies 7C the original SLICOT implementation. This file itself is *NOT* 8C part of SLICOT and the authors from NICONET e.V. are *NOT* 9C responsible for it. See file DESCRIPTION for the current 10C maintainer of the Octave control package. 11C 12C In altered implementation the deprecated LAPACK routine DLATZM 13C is replaced by ODLTZM. 14C 15C 16C SLICOT RELEASE 5.0. 17C 18C Copyright (c) 2002-2010 NICONET e.V. 19C 20C This program is free software: you can redistribute it and/or 21C modify it under the terms of the GNU General Public License as 22C published by the Free Software Foundation, either version 2 of 23C the License, or (at your option) any later version. 24C 25C This program is distributed in the hope that it will be useful, 26C but WITHOUT ANY WARRANTY; without even the implied warranty of 27C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 28C GNU General Public License for more details. 29C 30C You should have received a copy of the GNU General Public License 31C along with this program. If not, see 32C <http://www.gnu.org/licenses/>. 33C 34C PURPOSE 35C 36C To solve an N-by-N pole placement problem for the simple cases 37C N = 1 or N = 2: given the N-by-N matrix A and N-by-M matrix B, 38C construct an M-by-N matrix F such that A + B*F has prescribed 39C eigenvalues. These eigenvalues are specified by their sum S and 40C product P (if N = 2). The resulting F has minimum Frobenius norm. 41C 42C ARGUMENTS 43C 44C Input/Output Parameters 45C 46C N (input) INTEGER 47C The order of the matrix A and also the number of rows of 48C the matrix B and the number of columns of the matrix F. 49C N is either 1, if a single real eigenvalue is prescribed 50C or 2, if a complex conjugate pair or a set of two real 51C eigenvalues are prescribed. 52C 53C M (input) INTEGER 54C The number of columns of the matrix B and also the number 55C of rows of the matrix F. M >= 1. 56C 57C S (input) DOUBLE PRECISION 58C The sum of the prescribed eigenvalues if N = 2 or the 59C value of prescribed eigenvalue if N = 1. 60C 61C P (input) DOUBLE PRECISION 62C The product of the prescribed eigenvalues if N = 2. 63C Not referenced if N = 1. 64C 65C A (input/output) DOUBLE PRECISION array, dimension (N,N) 66C On entry, this array must contain the N-by-N state 67C dynamics matrix whose eigenvalues have to be moved to 68C prescribed locations. 69C On exit, this array contains no useful information. 70C 71C B (input/output) DOUBLE PRECISION array, dimension (N,M) 72C On entry, this array must contain the N-by-M input/state 73C matrix B. 74C On exit, this array contains no useful information. 75C 76C F (output) DOUBLE PRECISION array, dimension (M,N) 77C The state feedback matrix F which assigns one pole or two 78C poles of the closed-loop matrix A + B*F. 79C If N = 2 and the pair (A,B) is not controllable 80C (INFO = 1), then F(1,1) and F(1,2) contain the elements of 81C an orthogonal rotation which can be used to remove the 82C uncontrollable part of the pair (A,B). 83C 84C Tolerances 85C 86C TOL DOUBLE PRECISION 87C The absolute tolerance level below which the elements of A 88C and B are considered zero (used for controllability test). 89C 90C Workspace 91C 92C DWORK DOUBLE PRECISION array, dimension (M) 93C 94C Error Indicator 95C 96C INFO INTEGER 97C = 0: successful exit; 98C = 1: if uncontrollability of the pair (A,B) is detected. 99C 100C CONTRIBUTOR 101C 102C A. Varga, German Aerospace Center, 103C DLR Oberpfaffenhofen, July 1998. 104C Based on the RASP routine SB01BY. 105C 106C REVISIONS 107C 108C Nov. 1998, V. Sima, Research Institute for Informatics, Bucharest. 109C Dec. 1998, V. Sima, Katholieke Univ. Leuven, Leuven. 110C May 2003, A. Varga, German Aerospace Center. 111C 112C KEYWORDS 113C 114C Eigenvalue, eigenvalue assignment, feedback control, pole 115C placement, state-space model. 116C 117C ****************************************************************** 118C 119C .. Parameters .. 120 DOUBLE PRECISION FOUR, ONE, THREE, TWO, ZERO 121 PARAMETER ( FOUR = 4.0D0, ONE = 1.0D0, THREE = 3.0D0, 122 $ TWO = 2.0D0, ZERO = 0.0D0 ) 123C .. Scalar Arguments .. 124 INTEGER INFO, M, N 125 DOUBLE PRECISION P, S, TOL 126C .. Array Arguments .. 127 DOUBLE PRECISION A(N,*), B(N,*), DWORK(*), F(M,*) 128C .. Local Scalars .. 129 INTEGER IR, J 130 DOUBLE PRECISION ABSR, B1, B2, B21, C, C0, C1, C11, C12, C21, 131 $ C22, C3, C4, CS, CU, CV, DC0, DC2, DC3, DIFFR, 132 $ R, RN, S12, S21, SIG, SN, SU, SV, TAU1, TAU2, 133 $ WI, WI1, WR, WR1, X, Y, Z 134C .. External Functions .. 135 DOUBLE PRECISION DLAMC3, DLAMCH 136 EXTERNAL DLAMC3, DLAMCH 137C .. External Subroutines .. 138 EXTERNAL DLANV2, DLARFG, DLASET, DLASV2, ODLTZM, DROT 139C .. Intrinsic Functions .. 140 INTRINSIC ABS, MIN 141C .. Executable Statements .. 142C 143C For efficiency reasons, the parameters are not checked. 144C 145 INFO = 0 146 IF( N.EQ.1 ) THEN 147C 148C The case N = 1. 149C 150 IF( M.GT.1 ) 151 $ CALL DLARFG( M, B(1,1), B(1,2), N, TAU1 ) 152 B1 = B(1,1) 153 IF( ABS( B1 ).LE.TOL ) THEN 154C 155C The pair (A,B) is uncontrollable. 156C 157 INFO = 1 158 RETURN 159 END IF 160C 161 F(1,1) = ( S - A(1,1) )/B1 162 IF( M.GT.1 ) THEN 163 CALL DLASET( 'Full', M-1, 1, ZERO, ZERO, F(2,1), M ) 164 CALL ODLTZM( 'Left', M, N, B(1,2), N, TAU1, F(1,1), F(2,1), 165 $ M, DWORK ) 166 END IF 167 RETURN 168 END IF 169C 170C In the sequel N = 2. 171C 172C Compute the singular value decomposition of B in the form 173C 174C ( V 0 ) ( B1 0 ) 175C B = U*( G1 0 )*( )*H2*H1 , G1 = ( ), 176C ( 0 I ) ( 0 B2 ) 177C 178C ( CU SU ) ( CV SV ) 179C where U = ( ) and V = ( ) are orthogonal 180C (-SU CU ) (-SV CV ) 181C 182C rotations and H1 and H2 are elementary Householder reflectors. 183C ABS(B1) and ABS(B2) are the singular values of matrix B, 184C with ABS(B1) >= ABS(B2). 185C 186C Reduce first B to the lower bidiagonal form ( B1 0 ... 0 ). 187C ( B21 B2 ... 0 ) 188 IF( M.EQ.1 ) THEN 189C 190C Initialization for the case M = 1; no reduction required. 191C 192 B1 = B(1,1) 193 B21 = B(2,1) 194 B2 = ZERO 195 ELSE 196C 197C Postmultiply B with elementary Householder reflectors H1 198C and H2. 199C 200 CALL DLARFG( M, B(1,1), B(1,2), N, TAU1 ) 201 CALL ODLTZM( 'Right', N-1, M, B(1,2), N, TAU1, B(2,1), B(2,2), 202 $ N, DWORK ) 203 B1 = B(1,1) 204 B21 = B(2,1) 205 IF( M.GT.2 ) 206 $ CALL DLARFG( M-1, B(2,2), B(2,3), N, TAU2 ) 207 B2 = B(2,2) 208 END IF 209C 210C Reduce B to a diagonal form by premultiplying and postmultiplying 211C it with orthogonal rotations U and V, respectively, and order the 212C diagonal elements to have decreasing magnitudes. 213C Note: B2 has been set to zero if M = 1. Thus in the following 214C computations the case M = 1 need not to be distinguished. 215C Note also that LAPACK routine DLASV2 assumes an upper triangular 216C matrix, so the results should be adapted. 217C 218 CALL DLASV2( B1, B21, B2, X, Y, SU, CU, SV, CV ) 219 SU = -SU 220 B1 = Y 221 B2 = X 222C 223C Compute A1 = U'*A*U. 224C 225 CALL DROT( 2, A(2,1), 2, A(1,1), 2, CU, SU ) 226 CALL DROT( 2, A(1,2), 1, A(1,1), 1, CU, SU ) 227C 228C Compute the rank of B and check the controllability of the 229C pair (A,B). 230C 231 IR = 0 232 IF( ABS( B2 ).GT.TOL ) IR = IR + 1 233 IF( ABS( B1 ).GT.TOL ) IR = IR + 1 234 IF( IR.EQ.0 .OR. ( IR.EQ.1 .AND. ABS( A(2,1) ).LE.TOL ) ) THEN 235 F(1,1) = CU 236 F(1,2) = -SU 237C 238C The pair (A,B) is uncontrollable. 239C 240 INFO = 1 241 RETURN 242 END IF 243C 244C Compute F1 which assigns N poles for the reduced pair (A1,G1). 245C 246 X = DLAMC3( B1, B2 ) 247 IF( X.EQ.B1 ) THEN 248C 249C Rank one G1. 250C 251 F(1,1) = ( S - ( A(1,1) + A(2,2) ) )/B1 252 F(1,2) = -( A(2,2)*( A(2,2) - S ) + A(2,1)*A(1,2) + P )/ 253 $ A(2,1)/B1 254 IF( M.GT.1 ) THEN 255 F(2,1) = ZERO 256 F(2,2) = ZERO 257 END IF 258 ELSE 259C 260C Rank two G1. 261C 262 Z = ( S - ( A(1,1) + A(2,2) ) )/( B1*B1 + B2*B2 ) 263 F(1,1) = B1*Z 264 F(2,2) = B2*Z 265C 266C Compute an approximation for the minimum norm parameter 267C selection. 268C 269 X = A(1,1) + B1*F(1,1) 270 C = X*( S - X ) - P 271 IF( C.GE.ZERO ) THEN 272 SIG = ONE 273 ELSE 274 SIG = -ONE 275 END IF 276 S12 = B1/B2 277 S21 = B2/B1 278 C11 = ZERO 279 C12 = ONE 280 C21 = SIG*S12*C 281 C22 = A(1,2) - SIG*S12*A(2,1) 282 CALL DLANV2( C11, C12, C21, C22, WR, WI, WR1, WI1, CS, SN ) 283 IF( ABS( WR - A(1,2) ).GT.ABS( WR1 - A(1,2) ) ) THEN 284 R = WR1 285 ELSE 286 R = WR 287 END IF 288C 289C Perform Newton iteration to solve the equation for minimum. 290C 291 C0 = -C*C 292 C1 = C*A(2,1) 293 C4 = S21*S21 294 C3 = -C4*A(1,2) 295 DC0 = C1 296 DC2 = THREE*C3 297 DC3 = FOUR*C4 298C 299 DO 10 J = 1, 10 300 X = C0 + R*( C1 + R*R*( C3 + R*C4 ) ) 301 Y = DC0 + R*R*( DC2 + R*DC3 ) 302 IF( Y.EQ.ZERO ) GO TO 20 303 RN = R - X/Y 304 ABSR = ABS( R ) 305 DIFFR = ABS( R - RN ) 306 Z = DLAMC3( ABSR, DIFFR ) 307 IF( Z.EQ.ABSR ) 308 $ GO TO 20 309 R = RN 310 10 CONTINUE 311C 312 20 CONTINUE 313 IF( R.EQ.ZERO ) R = DLAMCH( 'Epsilon' ) 314 F(1,2) = ( R - A(1,2) )/B1 315 F(2,1) = ( C/R - A(2,1) )/B2 316 END IF 317C 318C Back-transform F1. Compute first F1*U'. 319C 320 CALL DROT( MIN( M, 2 ), F(1,1), 1, F(1,2), 1, CU, SU ) 321 IF( M.EQ.1 ) 322 $ RETURN 323C 324C Compute V'*F1. 325C 326 CALL DROT( 2, F(2,1), M, F(1,1), M, CV, SV ) 327C 328C ( F1 ) 329C Form F = ( ) . 330C ( 0 ) 331C 332 IF( M.GT.N ) 333 $ CALL DLASET( 'Full', M-N, N, ZERO, ZERO, F(N+1,1), M ) 334C 335C Compute H1*H2*F. 336C 337 IF( M.GT.2 ) 338 $ CALL ODLTZM( 'Left', M-1, N, B(2,3), N, TAU2, F(2,1), F(3,1), 339 $ M, DWORK ) 340 CALL ODLTZM( 'Left', M, N, B(1,2), N, TAU1, F(1,1), F(2,1), M, 341 $ DWORK ) 342C 343 RETURN 344C *** Last line of SB01BY *** 345 END 346