1 SUBROUTINE SB10QD( N, M, NP, NCON, NMEAS, GAMMA, A, LDA, B, LDB, 2 $ C, LDC, D, LDD, F, LDF, H, LDH, X, LDX, Y, LDY, 3 $ XYCOND, IWORK, DWORK, LDWORK, BWORK, INFO ) 4C 5C RELEASE 4.0, WGS COPYRIGHT 1999. 6C 7C PURPOSE 8C 9C To compute the state feedback and the output injection 10C matrices for an H-infinity (sub)optimal n-state controller, 11C using Glover's and Doyle's 1988 formulas, for the system 12C 13C | A | B1 B2 | | A | B | 14C P = |----|---------| = |---|---| 15C | C1 | D11 D12 | | C | D | 16C | C2 | D21 D22 | 17C 18C and for a given value of gamma, where B2 has as column size the 19C number of control inputs (NCON) and C2 has as row size the number 20C of measurements (NMEAS) being provided to the controller. 21C 22C It is assumed that 23C 24C (A1) (A,B2) is stabilizable and (C2,A) is detectable, 25C 26C (A2) D12 is full column rank with D12 = | 0 | and D21 is 27C | I | 28C full row rank with D21 = | 0 I | as obtained by the 29C subroutine SB10PD, 30C 31C (A3) | A-j*omega*I B2 | has full column rank for all omega, 32C | C1 D12 | 33C 34C 35C (A4) | A-j*omega*I B1 | has full row rank for all omega. 36C | C2 D21 | 37C 38C 39C ARGUMENTS 40C 41C Input/Output Parameters 42C 43C N (input) INTEGER 44C The order of the system. N >= 0. 45C 46C M (input) INTEGER 47C The column size of the matrix B. M >= 0. 48C 49C NP (input) INTEGER 50C The row size of the matrix C. NP >= 0. 51C 52C NCON (input) INTEGER 53C The number of control inputs (M2). M >= NCON >= 0, 54C NP-NMEAS >= NCON. 55C 56C NMEAS (input) INTEGER 57C The number of measurements (NP2). NP >= NMEAS >= 0, 58C M-NCON >= NMEAS. 59C 60C GAMMA (input) DOUBLE PRECISION 61C The value of gamma. It is assumed that gamma is 62C sufficiently large so that the controller is admissible. 63C GAMMA >= 0. 64C 65C A (input) DOUBLE PRECISION array, dimension (LDA,N) 66C The leading N-by-N part of this array must contain the 67C system state matrix A. 68C 69C LDA INTEGER 70C The leading dimension of the array A. LDA >= max(1,N). 71C 72C B (input) DOUBLE PRECISION array, dimension (LDB,M) 73C The leading N-by-M part of this array must contain the 74C system input matrix B. 75C 76C LDB INTEGER 77C The leading dimension of the array B. LDB >= max(1,N). 78C 79C C (input) DOUBLE PRECISION array, dimension (LDC,N) 80C The leading NP-by-N part of this array must contain the 81C system output matrix C. 82C 83C LDC INTEGER 84C The leading dimension of the array C. LDC >= max(1,NP). 85C 86C D (input) DOUBLE PRECISION array, dimension (LDD,M) 87C The leading NP-by-M part of this array must contain the 88C system input/output matrix D. 89C 90C LDD INTEGER 91C The leading dimension of the array D. LDD >= max(1,NP). 92C 93C F (output) DOUBLE PRECISION array, dimension (LDF,N) 94C The leading M-by-N part of this array contains the state 95C feedback matrix F. 96C 97C LDF INTEGER 98C The leading dimension of the array F. LDF >= max(1,M). 99C 100C H (output) DOUBLE PRECISION array, dimension (LDH,NP) 101C The leading N-by-NP part of this array contains the output 102C injection matrix H. 103C 104C LDH INTEGER 105C The leading dimension of the array H. LDH >= max(1,N). 106C 107C X (output) DOUBLE PRECISION array, dimension (LDX,N) 108C The leading N-by-N part of this array contains the matrix 109C X, solution of the X-Riccati equation. 110C 111C LDX INTEGER 112C The leading dimension of the array X. LDX >= max(1,N). 113C 114C Y (output) DOUBLE PRECISION array, dimension (LDY,N) 115C The leading N-by-N part of this array contains the matrix 116C Y, solution of the Y-Riccati equation. 117C 118C LDY INTEGER 119C The leading dimension of the array Y. LDY >= max(1,N). 120C 121C XYCOND (output) DOUBLE PRECISION array, dimension (2) 122C XYCOND(1) contains an estimate of the reciprocal condition 123C number of the X-Riccati equation; 124C XYCOND(2) contains an estimate of the reciprocal condition 125C number of the Y-Riccati equation. 126C 127C Workspace 128C 129C IWORK INTEGER array, dimension max(2*max(N,M-NCON,NP-NMEAS),N*N) 130C 131C DWORK DOUBLE PRECISION array, dimension (LDWORK) 132C On exit, if INFO = 0, DWORK(1) contains the optimal 133C LDWORK. 134C 135C LDWORK INTEGER 136C The dimension of the array DWORK. 137C LDWORK >= max(1,M*M + max(2*M1,3*N*N + 138C max(N*M,10*N*N+12*N+5)), 139C NP*NP + max(2*NP1,3*N*N + 140C max(N*NP,10*N*N+12*N+5))), 141C where M1 = M - M2 and NP1 = NP - NP2. 142C For good performance, LDWORK must generally be larger. 143C Denoting Q = MAX(M1,M2,NP1,NP2), an upper bound is 144C max(1,4*Q*Q+max(2*Q,3*N*N + max(2*N*Q,10*N*N+12*N+5))). 145C 146C BWORK LOGICAL array, dimension (2*N) 147C 148C Error Indicator 149C 150C INFO INTEGER 151C = 0: successful exit; 152C < 0: if INFO = -i, the i-th argument had an illegal 153C value; 154C = 1: if the controller is not admissible (too small value 155C of gamma); 156C = 2: if the X-Riccati equation was not solved 157C successfully (the controller is not admissible or 158C there are numerical difficulties); 159C = 3: if the Y-Riccati equation was not solved 160C successfully (the controller is not admissible or 161C there are numerical difficulties). 162C 163C METHOD 164C 165C The routine implements the Glover's and Doyle's formulas [1],[2] 166C modified as described in [3]. The X- and Y-Riccati equations 167C are solved with condition and accuracy estimates [4]. 168C 169C REFERENCES 170C 171C [1] Glover, K. and Doyle, J.C. 172C State-space formulae for all stabilizing controllers that 173C satisfy an Hinf norm bound and relations to risk sensitivity. 174C Systems and Control Letters, vol. 11, pp. 167-172, 1988. 175C 176C [2] Balas, G.J., Doyle, J.C., Glover, K., Packard, A., and 177C Smith, R. 178C mu-Analysis and Synthesis Toolbox. 179C The MathWorks Inc., Natick, Mass., 1995. 180C 181C [3] Petkov, P.Hr., Gu, D.W., and Konstantinov, M.M. 182C Fortran 77 routines for Hinf and H2 design of continuous-time 183C linear control systems. 184C Rep. 98-14, Department of Engineering, Leicester University, 185C Leicester, U.K., 1998. 186C 187C [4] Petkov, P.Hr., Konstantinov, M.M., and Mehrmann, V. 188C DGRSVX and DMSRIC: Fortan 77 subroutines for solving 189C continuous-time matrix algebraic Riccati equations with 190C condition and accuracy estimates. 191C Preprint SFB393/98-16, Fak. f. Mathematik, Tech. Univ. 192C Chemnitz, May 1998. 193C 194C NUMERICAL ASPECTS 195C 196C The precision of the solution of the matrix Riccati equations 197C can be controlled by the values of the condition numbers 198C XYCOND(1) and XYCOND(2) of these equations. 199C 200C FURTHER COMMENTS 201C 202C The Riccati equations are solved by the Schur approach 203C implementing condition and accuracy estimates. 204C 205C CONTRIBUTORS 206C 207C P.Hr. Petkov, D.W. Gu and M.M. Konstantinov, October 1998. 208C 209C REVISIONS 210C 211C V. Sima, Research Institute for Informatics, Bucharest, May 1999, 212C Sept. 1999. 213C 214C KEYWORDS 215C 216C Algebraic Riccati equation, H-infinity optimal control, robust 217C control. 218C 219C ****************************************************************** 220C 221C .. Parameters .. 222 DOUBLE PRECISION ZERO, ONE 223 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 224C 225C .. Scalar Arguments .. 226 INTEGER INFO, LDA, LDB, LDC, LDD, LDF, LDH, LDWORK, 227 $ LDX, LDY, M, N, NCON, NMEAS, NP 228 DOUBLE PRECISION GAMMA 229C .. 230C .. Array Arguments .. 231 INTEGER IWORK( * ) 232 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), 233 $ D( LDD, * ), DWORK( * ), F( LDF, * ), 234 $ H( LDH, * ), X( LDX, * ), XYCOND( 2 ), 235 $ Y( LDY, * ) 236 LOGICAL BWORK( * ) 237C 238C .. 239C .. Local Scalars .. 240 INTEGER INFO2, IW2, IWA, IWG, IWI, IWQ, IWR, IWRK, IWS, 241 $ IWT, IWV, LWAMAX, M1, M2, MINWRK, N2, ND1, ND2, 242 $ NN, NP1, NP2 243 DOUBLE PRECISION ANORM, EPS, FERR, RCOND, SEP 244C .. 245C .. External Functions .. 246C 247 DOUBLE PRECISION DLAMCH, DLANSY 248 EXTERNAL DLAMCH, DLANSY 249C .. 250C .. External Subroutines .. 251 EXTERNAL DGEMM, DLACPY, DLASET, DSYCON, DSYMM, DSYRK, 252 $ DSYTRF, DSYTRI, MB01RU, MB01RX, SB02RD, XERBLA 253C .. 254C .. Intrinsic Functions .. 255 INTRINSIC DBLE, INT, MAX 256C .. 257C .. Executable Statements .. 258C 259C Decode and Test input parameters. 260C 261 M1 = M - NCON 262 M2 = NCON 263 NP1 = NP - NMEAS 264 NP2 = NMEAS 265 NN = N*N 266C 267 INFO = 0 268 IF( N.LT.0 ) THEN 269 INFO = -1 270 ELSE IF( M.LT.0 ) THEN 271 INFO = -2 272 ELSE IF( NP.LT.0 ) THEN 273 INFO = -3 274 ELSE IF( NCON.LT.0 .OR. M1.LT.0 .OR. M2.GT.NP1 ) THEN 275 INFO = -4 276 ELSE IF( NMEAS.LT.0 .OR. NP1.LT.0 .OR. NP2.GT.M1 ) THEN 277 INFO = -5 278 ELSE IF( GAMMA.LT.ZERO ) THEN 279 INFO = -6 280 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 281 INFO = -8 282 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 283 INFO = -10 284 ELSE IF( LDC.LT.MAX( 1, NP ) ) THEN 285 INFO = -12 286 ELSE IF( LDD.LT.MAX( 1, NP ) ) THEN 287 INFO = -14 288 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 289 INFO = -16 290 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN 291 INFO = -18 292 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 293 INFO = -20 294 ELSE IF( LDY.LT.MAX( 1, N ) ) THEN 295 INFO = -22 296 ELSE 297C 298C Compute workspace. 299C 300 MINWRK = MAX( 1, M*M + MAX( 2*M1, 3*NN + 301 $ MAX( N*M, 10*NN + 12*N + 5 ) ), 302 $ NP*NP + MAX( 2*NP1, 3*NN + 303 $ MAX( N*NP, 10*NN + 12*N + 5 ) ) ) 304 IF( LDWORK.LT.MINWRK ) 305 $ INFO = -26 306 END IF 307 IF( INFO.NE.0 ) THEN 308 CALL XERBLA( 'SB10QD', -INFO ) 309 RETURN 310 END IF 311C 312C Quick return if possible. 313C 314 IF( N.EQ.0 .OR. M.EQ.0 .OR. NP.EQ.0 .OR. M1.EQ.0 .OR. M2.EQ.0 315 $ .OR. NP1.EQ.0 .OR. NP2.EQ.0 ) THEN 316 XYCOND( 1 ) = ONE 317 XYCOND( 2 ) = ONE 318 DWORK( 1 ) = ONE 319 RETURN 320 END IF 321 ND1 = NP1 - M2 322 ND2 = M1 - NP2 323 N2 = 2*N 324C 325C Get the machine precision. 326C 327 EPS = DLAMCH( 'Epsilon' ) 328C 329C Workspace usage. 330C 331 IWA = M*M + 1 332 IWQ = IWA + NN 333 IWG = IWQ + NN 334 IW2 = IWG + NN 335C 336C Compute |D1111'||D1111 D1112| - gamma^2*Im1 . 337C |D1112'| 338C 339 CALL DLASET( 'L', M1, M1, ZERO, -GAMMA*GAMMA, DWORK, M ) 340 IF( ND1.GT.0 ) 341 $ CALL DSYRK( 'L', 'T', M1, ND1, ONE, D, LDD, ONE, DWORK, M ) 342C 343C Compute inv(|D1111'|*|D1111 D1112| - gamma^2*Im1) . 344C |D1112'| 345C 346 IWRK = IWA 347 ANORM = DLANSY( 'I', 'L', M1, DWORK, M, DWORK( IWRK ) ) 348 CALL DSYTRF( 'L', M1, DWORK, M, IWORK, DWORK( IWRK ), 349 $ LDWORK-IWRK+1, INFO2 ) 350 IF( INFO2.GT.0 ) THEN 351 INFO = 1 352 RETURN 353 END IF 354C 355 LWAMAX = INT( DWORK( IWRK ) ) + IWRK - 1 356 CALL DSYCON( 'L', M1, DWORK, M, IWORK, ANORM, RCOND, 357 $ DWORK( IWRK ), IWORK( M1+1 ), INFO2 ) 358 IF( RCOND.LT.EPS ) THEN 359 INFO = 1 360 RETURN 361 END IF 362C 363C Compute inv(R) block by block. 364C 365 CALL DSYTRI( 'L', M1, DWORK, M, IWORK, DWORK( IWRK ), INFO2 ) 366C 367C Compute -|D1121 D1122|*inv(|D1111'|*|D1111 D1112| - gamma^2*Im1) . 368C |D1112'| 369C 370 CALL DSYMM( 'R', 'L', M2, M1, -ONE, DWORK, M, D( ND1+1, 1 ), LDD, 371 $ ZERO, DWORK( M1+1 ), M ) 372C 373C Compute |D1121 D1122|*inv(|D1111'|*|D1111 D1112| - 374C |D1112'| 375C 376C gamma^2*Im1)*|D1121'| + Im2 . 377C |D1122'| 378C 379 CALL DLASET( 'Lower', M2, M2, ZERO, ONE, DWORK( M1*(M+1)+1 ), M ) 380 CALL MB01RX( 'Right', 'Lower', 'Transpose', M2, M1, ONE, -ONE, 381 $ DWORK( M1*(M+1)+1 ), M, D( ND1+1, 1 ), LDD, 382 $ DWORK( M1+1 ), M, INFO2 ) 383C 384C Compute D11'*C1 . 385C 386 CALL DGEMM( 'T', 'N', M1, N, NP1, ONE, D, LDD, C, LDC, ZERO, 387 $ DWORK( IW2 ), M ) 388C 389C Compute D1D'*C1 . 390C 391 CALL DLACPY( 'Full', M2, N, C( ND1+1, 1 ), LDC, DWORK( IW2+M1 ), 392 $ M ) 393C 394C Compute inv(R)*D1D'*C1 in F . 395C 396 CALL DSYMM( 'L', 'L', M, N, ONE, DWORK, M, DWORK( IW2 ), M, ZERO, 397 $ F, LDF ) 398C 399C Compute Ax = A - B*inv(R)*D1D'*C1 . 400C 401 CALL DLACPY( 'Full', N, N, A, LDA, DWORK( IWA ), N ) 402 CALL DGEMM( 'N', 'N', N, N, M, -ONE, B, LDB, F, LDF, ONE, 403 $ DWORK( IWA ), N ) 404C 405C Compute Cx = C1'*C1 - C1'*D1D*inv(R)*D1D'*C1 . 406C 407 IF( ND1.EQ.0 ) THEN 408 CALL DLASET( 'L', N, N, ZERO, ZERO, DWORK( IWQ ), N ) 409 ELSE 410 CALL DSYRK( 'L', 'T', N, NP1, ONE, C, LDC, ZERO, 411 $ DWORK( IWQ ), N ) 412 CALL MB01RX( 'Left', 'Lower', 'Transpose', N, M, ONE, -ONE, 413 $ DWORK( IWQ ), N, DWORK( IW2 ), M, F, LDF, INFO2 ) 414 END IF 415C 416C Compute Dx = B*inv(R)*B' . 417C 418 IWRK = IW2 419 CALL MB01RU( 'Lower', 'NoTranspose', N, M, ZERO, ONE, 420 $ DWORK( IWG ), N, B, LDB, DWORK, M, DWORK( IWRK ), 421 $ M*N, INFO2 ) 422C 423C Solution of the Riccati equation Ax'*X + X*Ax + Cx - X*Dx*X = 0 . 424C Workspace: need M*M + 13*N*N + 12*N + 5; 425C prefer larger. 426C 427 IWT = IW2 428 IWV = IWT + NN 429 IWR = IWV + NN 430 IWI = IWR + N2 431 IWS = IWI + N2 432 IWRK = IWS + 4*NN 433C 434 CALL SB02RD( 'All', 'Continuous', 'NotUsed', 'NoTranspose', 435 $ 'Lower', 'GeneralScaling', 'Stable', 'NotFactored', 436 $ 'Original', N, DWORK( IWA ), N, DWORK( IWT ), N, 437 $ DWORK( IWV ), N, DWORK( IWG ), N, DWORK( IWQ ), N, 438 $ X, LDX, SEP, XYCOND( 1 ), FERR, DWORK( IWR ), 439 $ DWORK( IWI ), DWORK( IWS ), N2, IWORK, DWORK( IWRK ), 440 $ LDWORK-IWRK+1, BWORK, INFO2 ) 441 IF( INFO2.GT.0 ) THEN 442 INFO = 2 443 RETURN 444 END IF 445C 446 LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX ) 447C 448C Compute F = -inv(R)*|D1D'*C1 + B'*X| . 449C 450 IWRK = IW2 451 CALL DGEMM( 'T', 'N', M, N, N, ONE, B, LDB, X, LDX, ZERO, 452 $ DWORK( IWRK ), M ) 453 CALL DSYMM( 'L', 'L', M, N, -ONE, DWORK, M, DWORK( IWRK ), M, 454 $ -ONE, F, LDF ) 455C 456C Workspace usage. 457C 458 IWA = NP*NP + 1 459 IWQ = IWA + NN 460 IWG = IWQ + NN 461 IW2 = IWG + NN 462C 463C Compute |D1111|*|D1111' D1121'| - gamma^2*Inp1 . 464C |D1121| 465C 466 CALL DLASET( 'U', NP1, NP1, ZERO, -GAMMA*GAMMA, DWORK, NP ) 467 IF( ND2.GT.0 ) 468 $ CALL DSYRK( 'U', 'N', NP1, ND2, ONE, D, LDD, ONE, DWORK, NP ) 469C 470C Compute inv(|D1111|*|D1111' D1121'| - gamma^2*Inp1) . 471C |D1121| 472C 473 IWRK = IWA 474 ANORM = DLANSY( 'I', 'U', NP1, DWORK, NP, DWORK( IWRK ) ) 475 CALL DSYTRF( 'U', NP1, DWORK, NP, IWORK, DWORK( IWRK ), 476 $ LDWORK-IWRK+1, INFO2 ) 477 IF( INFO2.GT.0 ) THEN 478 INFO = 1 479 RETURN 480 END IF 481C 482 LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX ) 483 CALL DSYCON( 'U', NP1, DWORK, NP, IWORK, ANORM, RCOND, 484 $ DWORK( IWRK ), IWORK( NP1+1 ), INFO2 ) 485 IF( RCOND.LT.EPS ) THEN 486 INFO = 1 487 RETURN 488 END IF 489C 490C Compute inv(RT) . 491C 492 CALL DSYTRI( 'U', NP1, DWORK, NP, IWORK, DWORK( IWRK ), INFO2 ) 493C 494C Compute -inv(|D1111||D1111' D1121'| - gamma^2*Inp1)*|D1112| . 495C |D1121| |D1122| 496C 497 CALL DSYMM( 'L', 'U', NP1, NP2, -ONE, DWORK, NP, D( 1, ND2+1 ), 498 $ LDD, ZERO, DWORK( NP1*NP+1 ), NP ) 499C 500C Compute [D1112' D1122']*inv(|D1111||D1111' D1121'| - 501C |D1121| 502C 503C gamma^2*Inp1)*|D1112| + Inp2 . 504C |D1122| 505C 506 CALL DLASET( 'Full', NP2, NP2, ZERO, ONE, DWORK( NP1*(NP+1)+1 ), 507 $ NP ) 508 CALL MB01RX( 'Left', 'Upper', 'Transpose', NP2, NP1, ONE, -ONE, 509 $ DWORK( NP1*(NP+1)+1 ), NP, D( 1, ND2+1 ), LDD, 510 $ DWORK( NP1*NP+1 ), NP, INFO2 ) 511C 512C Compute B1*D11' . 513C 514 CALL DGEMM( 'N', 'T', N, NP1, M1, ONE, B, LDB, D, LDD, ZERO, 515 $ DWORK( IW2 ), N ) 516C 517C Compute B1*DD1' . 518C 519 CALL DLACPY( 'Full', N, NP2, B( 1, ND2+1 ), LDB, 520 $ DWORK( IW2+NP1*N ), N ) 521C 522C Compute B1*DD1'*inv(RT) in H . 523C 524 CALL DSYMM( 'R', 'U', N, NP, ONE, DWORK, NP, DWORK( IW2 ), N, 525 $ ZERO, H, LDH ) 526C 527C Compute Ay = A - B1*DD1'*inv(RT)*C . 528C 529 CALL DLACPY( 'Full', N, N, A, LDA, DWORK( IWA ), N ) 530 CALL DGEMM( 'N', 'N', N, N, NP, -ONE, H, LDH, C, LDC, ONE, 531 $ DWORK( IWA ), N ) 532C 533C Compute Cy = B1*B1' - B1*DD1'*inv(RT)*DD1*B1' . 534C 535 IF( ND2.EQ.0 ) THEN 536 CALL DLASET( 'U', N, N, ZERO, ZERO, DWORK( IWQ ), N ) 537 ELSE 538 CALL DSYRK( 'U', 'N', N, M1, ONE, B, LDB, ZERO, DWORK( IWQ ), 539 $ N ) 540 CALL MB01RX( 'Right', 'Upper', 'Transpose', N, NP, ONE, -ONE, 541 $ DWORK( IWQ ), N, H, LDH, DWORK( IW2 ), N, INFO2 ) 542 END IF 543C 544C Compute Dy = C'*inv(RT)*C . 545C 546 IWRK = IW2 547 CALL MB01RU( 'Upper', 'Transpose', N, NP, ZERO, ONE, DWORK( IWG ), 548 $ N, C, LDC, DWORK, NP, DWORK( IWRK), N*NP, INFO2 ) 549C 550C Solution of the Riccati equation Ay*Y + Y*Ay' + Cy - Y*Dy*Y = 0 . 551C Workspace: need NP*NP + 13*N*N + 12*N + 5; 552C prefer larger. 553C 554 IWT = IW2 555 IWV = IWT + NN 556 IWR = IWV + NN 557 IWI = IWR + N2 558 IWS = IWI + N2 559 IWRK = IWS + 4*NN 560C 561 CALL SB02RD( 'All', 'Continuous', 'NotUsed', 'Transpose', 562 $ 'Upper', 'GeneralScaling', 'Stable', 'NotFactored', 563 $ 'Original', N, DWORK( IWA ), N, DWORK( IWT ), N, 564 $ DWORK( IWV ), N, DWORK( IWG ), N, DWORK( IWQ ), N, 565 $ Y, LDY, SEP, XYCOND( 2 ), FERR, DWORK( IWR ), 566 $ DWORK( IWI ), DWORK( IWS ), N2, IWORK, DWORK( IWRK ), 567 $ LDWORK-IWRK+1, BWORK, INFO2 ) 568 IF( INFO2.GT.0 ) THEN 569 INFO = 3 570 RETURN 571 END IF 572C 573 LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX ) 574C 575C Compute H = -|B1*DD1' + Y*C'|*inv(RT) . 576C 577 IWRK = IW2 578 CALL DGEMM( 'N', 'T', N, NP, N, ONE, Y, LDY, C, LDC, ZERO, 579 $ DWORK( IWRK ), N ) 580 CALL DSYMM( 'R', 'U', N, NP, -ONE, DWORK, NP, DWORK( IWRK ), N, 581 $ -ONE, H, LDH ) 582C 583 DWORK( 1 ) = DBLE( LWAMAX ) 584 RETURN 585C *** Last line of SB10QD *** 586 END 587