1 SUBROUTINE SB04OD( REDUCE, TRANS, JOBD, M, N, A, LDA, B, LDB, C, 2 $ LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, P, 3 $ LDP, Q, LDQ, U, LDU, V, LDV, IWORK, DWORK, 4 $ LDWORK, INFO ) 5C 6C WARNING 7C 8C This file alters the SLICOT routine SB04OD. It is intended 9C for use from the Octave Control Systems Package and modifies 10C the original SLICOT implementation. This file itself is *NOT* 11C part of SLICOT and the authors from NICONET e.V. are *NOT* 12C responsible for it. See file DESCRIPTION for the current 13C maintainer of the Octave control package. 14C 15C In altered implementation the deprecated LAPACK routine DGEGS 16C is replaced by DGGES. 17C 18C 19C SLICOT RELEASE 5.0. 20C 21C Copyright (c) 2002-2010 NICONET e.V. 22C 23C This program is free software: you can redistribute it and/or 24C modify it under the terms of the GNU General Public License as 25C published by the Free Software Foundation, either version 2 of 26C the License, or (at your option) any later version. 27C 28C This program is distributed in the hope that it will be useful, 29C but WITHOUT ANY WARRANTY; without even the implied warranty of 30C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 31C GNU General Public License for more details. 32C 33C You should have received a copy of the GNU General Public License 34C along with this program. If not, see 35C <http://www.gnu.org/licenses/>. 36C 37C PURPOSE 38C 39C To solve for R and L one of the generalized Sylvester equations 40C 41C A * R - L * B = scale * C ) 42C ) (1) 43C D * R - L * E = scale * F ) 44C 45C or 46C 47C A' * R + D' * L = scale * C ) 48C ) (2) 49C R * B' + L * E' = scale * (-F) ) 50C 51C where A and D are M-by-M matrices, B and E are N-by-N matrices and 52C C, F, R and L are M-by-N matrices. 53C 54C The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an 55C output scaling factor chosen to avoid overflow. 56C 57C The routine also optionally computes a Dif estimate, which 58C measures the separation of the spectrum of the matrix pair (A,D) 59C from the spectrum of the matrix pair (B,E), Dif[(A,D),(B,E)]. 60C 61C ARGUMENTS 62C 63C MODE PARAMETERS 64C 65C REDUCE CHARACTER*1 66C Indicates whether the matrix pairs (A,D) and/or (B,E) are 67C to be reduced to generalized Schur form as follows: 68C = 'R': The matrix pairs (A,D) and (B,E) are to be reduced 69C to generalized (real) Schur canonical form; 70C = 'A': The matrix pair (A,D) only is to be reduced 71C to generalized (real) Schur canonical form, 72C and the matrix pair (B,E) already is in this form; 73C = 'B': The matrix pair (B,E) only is to be reduced 74C to generalized (real) Schur canonical form, 75C and the matrix pair (A,D) already is in this form; 76C = 'N': The matrix pairs (A,D) and (B,E) are already in 77C generalized (real) Schur canonical form, as 78C produced by LAPACK routine DGEES. 79C 80C TRANS CHARACTER*1 81C Indicates which of the equations, (1) or (2), is to be 82C solved as follows: 83C = 'N': The generalized Sylvester equation (1) is to be 84C solved; 85C = 'T': The "transposed" generalized Sylvester equation 86C (2) is to be solved. 87C 88C JOBD CHARACTER*1 89C Indicates whether the Dif estimator is to be computed as 90C follows: 91C = '1': Only the one-norm-based Dif estimate is computed 92C and stored in DIF; 93C = '2': Only the Frobenius norm-based Dif estimate is 94C computed and stored in DIF; 95C = 'D': The equation (1) is solved and the one-norm-based 96C Dif estimate is computed and stored in DIF; 97C = 'F': The equation (1) is solved and the Frobenius norm- 98C based Dif estimate is computed and stored in DIF; 99C = 'N': The Dif estimator is not required and hence DIF is 100C not referenced. (Solve either (1) or (2) only.) 101C JOBD is not referenced if TRANS = 'T'. 102C 103C Input/Output Parameters 104C 105C M (input) INTEGER 106C The order of the matrices A and D and the number of rows 107C of the matrices C, F, R and L. M >= 0. 108C 109C N (input) INTEGER 110C The order of the matrices B and E and the number of 111C columns of the matrices C, F, R and L. N >= 0. 112C 113C A (input/output) DOUBLE PRECISION array, dimension (LDA,M) 114C On entry, the leading M-by-M part of this array must 115C contain the coefficient matrix A of the equation; A must 116C be in upper quasi-triangular form if REDUCE = 'B' or 'N'. 117C On exit, the leading M-by-M part of this array contains 118C the upper quasi-triangular form of A. 119C 120C LDA INTEGER 121C The leading dimension of array A. LDA >= MAX(1,M). 122C 123C B (input/output) DOUBLE PRECISION array, dimension (LDB,N) 124C On entry, the leading N-by-N part of this array must 125C contain the coefficient matrix B of the equation; B must 126C be in upper quasi-triangular form if REDUCE = 'A' or 'N'. 127C On exit, the leading N-by-N part of this array contains 128C the upper quasi-triangular form of B. 129C 130C LDB INTEGER 131C The leading dimension of array B. LDB >= MAX(1,N). 132C 133C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 134C On entry, the leading M-by-N part of this array must 135C contain the right-hand side matrix C of the first equation 136C in (1) or (2). 137C On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N 138C part of this array contains the solution matrix R of the 139C problem; if JOBD = '1' or '2' and TRANS = 'N', the leading 140C M-by-N part of this array contains the solution matrix R 141C achieved during the computation of the Dif estimate. 142C 143C LDC INTEGER 144C The leading dimension of array C. LDC >= MAX(1,M). 145C 146C D (input/output) DOUBLE PRECISION array, dimension (LDD,M) 147C On entry, the leading M-by-M part of this array must 148C contain the coefficient matrix D of the equation; D must 149C be in upper triangular form if REDUCE = 'B' or 'N'. 150C On exit, the leading M-by-M part of this array contains 151C the upper triangular form of D. 152C 153C LDD INTEGER 154C The leading dimension of array D. LDD >= MAX(1,M). 155C 156C E (input/output) DOUBLE PRECISION array, dimension (LDE,N) 157C On entry, the leading N-by-N part of this array must 158C contain the coefficient matrix E of the equation; E must 159C be in upper triangular form if REDUCE = 'A' or 'N'. 160C On exit, the leading N-by-N part of this array contains 161C the upper triangular form of E. 162C 163C LDE INTEGER 164C The leading dimension of array E. LDE >= MAX(1,N). 165C 166C F (input/output) DOUBLE PRECISION array, dimension (LDF,N) 167C On entry, the leading M-by-N part of this array must 168C contain the right-hand side matrix F of the second 169C equation in (1) or (2). 170C On exit, if JOBD = 'N', 'D' or 'F', the leading M-by-N 171C part of this array contains the solution matrix L of the 172C problem; if JOBD = '1' or '2' and TRANS = 'N', the leading 173C M-by-N part of this array contains the solution matrix L 174C achieved during the computation of the Dif estimate. 175C 176C LDF INTEGER 177C The leading dimension of array F. LDF >= MAX(1,M). 178C 179C SCALE (output) DOUBLE PRECISION 180C The scaling factor in (1) or (2). If 0 < SCALE < 1, C and 181C F hold the solutions R and L, respectively, to a slightly 182C perturbed system (but the input or computed generalized 183C (real) Schur canonical form matrices A, B, D, and E 184C have not been changed). If SCALE = 0, C and F hold the 185C solutions R and L, respectively, to the homogeneous system 186C with C = F = 0. Normally, SCALE = 1. 187C 188C DIF (output) DOUBLE PRECISION 189C If TRANS = 'N' and JOBD <> 'N', then DIF contains the 190C value of the Dif estimator, which is an upper bound of 191C -1 192C Dif[(A,D),(B,E)] = sigma_min(Z) = 1/||Z ||, in either the 193C one-norm, or Frobenius norm, respectively (see METHOD). 194C Otherwise, DIF is not referenced. 195C 196C P (output) DOUBLE PRECISION array, dimension (LDP,*) 197C If REDUCE = 'R' or 'A', then the leading M-by-M part of 198C this array contains the (left) transformation matrix used 199C to reduce (A,D) to generalized Schur form. 200C Otherwise, P is not referenced and can be supplied as a 201C dummy array (i.e. set parameter LDP = 1 and declare this 202C array to be P(1,1) in the calling program). 203C 204C LDP INTEGER 205C The leading dimension of array P. 206C LDP >= MAX(1,M) if REDUCE = 'R' or 'A', 207C LDP >= 1 if REDUCE = 'B' or 'N'. 208C 209C Q (output) DOUBLE PRECISION array, dimension (LDQ,*) 210C If REDUCE = 'R' or 'A', then the leading M-by-M part of 211C this array contains the (right) transformation matrix used 212C to reduce (A,D) to generalized Schur form. 213C Otherwise, Q is not referenced and can be supplied as a 214C dummy array (i.e. set parameter LDQ = 1 and declare this 215C array to be Q(1,1) in the calling program). 216C 217C LDQ INTEGER 218C The leading dimension of array Q. 219C LDQ >= MAX(1,M) if REDUCE = 'R' or 'A', 220C LDQ >= 1 if REDUCE = 'B' or 'N'. 221C 222C U (output) DOUBLE PRECISION array, dimension (LDU,*) 223C If REDUCE = 'R' or 'B', then the leading N-by-N part of 224C this array contains the (left) transformation matrix used 225C to reduce (B,E) to generalized Schur form. 226C Otherwise, U is not referenced and can be supplied as a 227C dummy array (i.e. set parameter LDU = 1 and declare this 228C array to be U(1,1) in the calling program). 229C 230C LDU INTEGER 231C The leading dimension of array U. 232C LDU >= MAX(1,N) if REDUCE = 'R' or 'B', 233C LDU >= 1 if REDUCE = 'A' or 'N'. 234C 235C V (output) DOUBLE PRECISION array, dimension (LDV,*) 236C If REDUCE = 'R' or 'B', then the leading N-by-N part of 237C this array contains the (right) transformation matrix used 238C to reduce (B,E) to generalized Schur form. 239C Otherwise, V is not referenced and can be supplied as a 240C dummy array (i.e. set parameter LDV = 1 and declare this 241C array to be V(1,1) in the calling program). 242C 243C LDV INTEGER 244C The leading dimension of array V. 245C LDV >= MAX(1,N) if REDUCE = 'R' or 'B', 246C LDV >= 1 if REDUCE = 'A' or 'N'. 247C 248C Workspace 249C 250C IWORK INTEGER array, dimension (M+N+6) 251C 252C DWORK DOUBLE PRECISION array, dimension (LDWORK) 253C On exit, if INFO = 0, DWORK(1) returns the optimal value 254C of LDWORK. 255C 256C LDWORK INTEGER 257C The length of the array DWORK. 258C If TRANS = 'N' and JOBD = 'D' or 'F', then 259C LDWORK = MAX(1,7*M,7*N,2*M*N) if REDUCE = 'R'; 260C LDWORK = MAX(1,7*M,2*M*N) if REDUCE = 'A'; 261C LDWORK = MAX(1,7*N,2*M*N) if REDUCE = 'B'; 262C LDWORK = MAX(1,2*M*N) if REDUCE = 'N'. 263C Otherwise, the term 2*M*N above should be omitted. 264C For optimum performance LDWORK should be larger. 265C 266C Error Indicator 267C 268C INFO INTEGER 269C = 0: successful exit; 270C < 0: if INFO = -i, the i-th argument had an illegal 271C value; 272C = 1: if REDUCE <> 'N' and either (A,D) and/or (B,E) 273C cannot be reduced to generalized Schur form; 274C = 2: if REDUCE = 'N' and either A or B is not in 275C upper quasi-triangular form; 276C = 3: if a singular matrix was encountered during the 277C computation of the solution matrices R and L, that 278C is (A,D) and (B,E) have common or close eigenvalues. 279C 280C METHOD 281C 282C For the case TRANS = 'N', and REDUCE = 'R' or 'N', the algorithm 283C used by the routine consists of four steps (see [1] and [2]) as 284C follows: 285C 286C (a) if REDUCE = 'R', then the matrix pairs (A,D) and (B,E) are 287C transformed to generalized Schur form, i.e. orthogonal 288C matrices P, Q, U and V are computed such that P' * A * Q 289C and U' * B * V are in upper quasi-triangular form and 290C P' * D * Q and U' * E * V are in upper triangular form; 291C (b) if REDUCE = 'R', then the matrices C and F are transformed 292C to give P' * C * V and P' * F * V respectively; 293C (c) if REDUCE = 'R', then the transformed system 294C 295C P' * A * Q * R1 - L1 * U' * B * V = scale * P' * C * V 296C P' * D * Q * R1 - L1 * U' * E * V = scale * P' * F * V 297C 298C is solved to give R1 and L1; otherwise, equation (1) is 299C solved to give R and L directly. The Dif estimator 300C is also computed if JOBD <> 'N'. 301C (d) if REDUCE = 'R', then the solution is transformed back 302C to give R = Q * R1 * V' and L = P * L1 * U'. 303C 304C By using Kronecker products, equation (1) can also be written as 305C the system of linear equations Z * x = scale*y (see [1]), where 306C 307C | I*A I*D | 308C Z = | |. 309C |-B'*I -E'*I | 310C 311C -1 312C If JOBD <> 'N', then a lower bound on ||Z ||, in either the one- 313C norm or Frobenius norm, is computed, which in most cases is 314C a reliable estimate of the true value. Notice that since Z is a 315C matrix of order 2 * M * N, the exact value of Dif (i.e., in the 316C Frobenius norm case, the smallest singular value of Z) may be very 317C expensive to compute. 318C 319C The case TRANS = 'N', and REDUCE = 'A' or 'B', is similar, but 320C only one of the matrix pairs should be reduced and the 321C calculations simplify. 322C 323C For the case TRANS = 'T', and REDUCE = 'R' or 'N', the algorithm 324C is similar, but the steps (b), (c), and (d) are as follows: 325C 326C (b) if REDUCE = 'R', then the matrices C and F are transformed 327C to give Q' * C * V and P' * F * U respectively; 328C (c) if REDUCE = 'R', then the transformed system 329C 330C Q' * A' * P * R1 + Q' * D' * P * L1 = scale * Q' * C * V 331C R1 * V' * B' * U + L1 * V' * E' * U = -scale * P' * F * U 332C 333C is solved to give R1 and L1; otherwise, equation (2) is 334C solved to give R and L directly. 335C (d) if REDUCE = 'R', then the solution is transformed back 336C to give R = P * R1 * V' and L = P * L1 * V'. 337C 338C REFERENCES 339C 340C [1] Kagstrom, B. and Westin, L. 341C Generalized Schur Methods with Condition Estimators for 342C Solving the Generalized Sylvester Equation. 343C IEEE Trans. Auto. Contr., 34, pp. 745-751, 1989. 344C [2] Kagstrom, B. and Westin, L. 345C GSYLV - Fortran Routines for the Generalized Schur Method with 346C Dif Estimators for Solving the Generalized Sylvester 347C Equation. 348C Report UMINF-132.86, Institute of Information Processing, 349C Univ. of Umea, Sweden, July 1987. 350C [3] Golub, G.H., Nash, S. and Van Loan, C.F. 351C A Hessenberg-Schur Method for the Problem AX + XB = C. 352C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. 353C [4] Kagstrom, B. and Van Dooren, P. 354C Additive Decomposition of a Transfer Function with respect to 355C a Specified Region. 356C In: "Signal Processing, Scattering and Operator Theory, and 357C Numerical Methods" (Eds. M.A. Kaashoek et al.). 358C Proceedings of MTNS-89, Vol. 3, pp. 469-477, Birkhauser Boston 359C Inc., 1990. 360C [5] Kagstrom, B. and Van Dooren, P. 361C A Generalized State-space Approach for the Additive 362C Decomposition of a Transfer Matrix. 363C Report UMINF-91.12, Institute of Information Processing, Univ. 364C of Umea, Sweden, April 1991. 365C 366C NUMERICAL ASPECTS 367C 368C The algorithm is backward stable. A reliable estimate for the 369C condition number of Z in the Frobenius norm, is (see [1]) 370C 371C K(Z) = SQRT( ||A||**2 + ||B||**2 + ||C||**2 + ||D||**2 )/DIF. 372C 373C If mu is an upper bound on the relative error of the elements of 374C the matrices A, B, C, D, E and F, then the relative error in the 375C actual solution is approximately mu * K(Z). 376C 377C The relative error in the computed solution (due to rounding 378C errors) is approximately EPS * K(Z), where EPS is the machine 379C precision (see LAPACK Library routine DLAMCH). 380C 381C FURTHER COMMENTS 382C 383C For applications of the generalized Sylvester equation in control 384C theory, see [4] and [5]. 385C 386C CONTRIBUTORS 387C 388C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Aug. 1997. 389C Supersedes Release 2.0 routine SB04CD by Bo Kagstrom and Lars 390C Westin. 391C 392C REVISIONS 393C 394C V. Sima, Katholieke Univ. Leuven, Belgium, May 1999, Dec. 1999, 395C May 2009. 396C 397C KEYWORDS 398C 399C Generalized eigenvalue problem, orthogonal transformation, real 400C Schur form, Sylvester equation. 401C 402C ****************************************************************** 403C 404C .. Parameters .. 405 DOUBLE PRECISION ZERO, ONE 406 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 407C .. Scalar Arguments .. 408 CHARACTER JOBD, REDUCE, TRANS 409 INTEGER INFO, LDA, LDB, LDC, LDD, LDE, LDF, LDP, LDQ, 410 $ LDU, LDV, LDWORK, M, N 411 DOUBLE PRECISION DIF, SCALE 412C .. Array Arguments .. 413 INTEGER IWORK(*) 414 DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), 415 $ DWORK(*), E(LDE,*), F(LDF,*), P(LDP,*), 416 $ Q(LDQ,*), U(LDU,*), V(LDV,*) 417C .. Local Scalars .. 418 LOGICAL ILASCL, ILBSCL, ILDSCL, ILESCL, LJOB1, LJOB2, 419 $ LJOBD, LJOBDF, LJOBF, LREDRA, LREDRB, LREDUA, 420 $ LREDUB, LREDUC, LREDUR, LTRANN, SUFWRK 421 INTEGER I, IERR, IJOB, MINWRK, MN, WRKOPT, SDIM 422 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, DNRM, 423 $ DNRMTO, ENRM, ENRMTO, SAFMAX, SAFMIN, SMLNUM 424C .. External Functions .. 425 LOGICAL LSAME 426 DOUBLE PRECISION DLAMCH, DLANGE 427 EXTERNAL DLAMCH, DLANGE, LSAME 428C .. External Subroutines .. 429 EXTERNAL DCOPY, DGGES, DGEMM, DGEMV, DLABAD, DLACPY, 430 $ DLASCL, DTGSYL, XERBLA 431C .. Intrinsic Functions .. 432 INTRINSIC INT, MAX, SQRT 433C .. Executable Statements .. 434C 435 INFO = 0 436 MN = MAX( M, N ) 437 LREDUR = LSAME( REDUCE, 'R' ) 438 LREDUA = LSAME( REDUCE, 'A' ) 439 LREDUB = LSAME( REDUCE, 'B' ) 440 LREDRA = LREDUR.OR.LREDUA 441 LREDRB = LREDUR.OR.LREDUB 442 LREDUC = LREDRA.OR.LREDUB 443 IF ( LREDUR ) THEN 444 MINWRK = MAX( 1, 7*MN ) 445 ELSE IF ( LREDUA ) THEN 446 MINWRK = MAX( 1, 7*M ) 447 ELSE IF ( LREDUB ) THEN 448 MINWRK = MAX( 1, 7*N ) 449 ELSE 450 MINWRK = 1 451 END IF 452 LTRANN = LSAME( TRANS, 'N' ) 453 IF ( LTRANN ) THEN 454 LJOB1 = LSAME( JOBD, '1' ) 455 LJOB2 = LSAME( JOBD, '2' ) 456 LJOBD = LSAME( JOBD, 'D' ) 457 LJOBF = LSAME( JOBD, 'F' ) 458 LJOBDF = LJOB1.OR.LJOB2.OR.LJOBD.OR.LJOBF 459 IF ( LJOBD.OR.LJOBF ) MINWRK = MAX( MINWRK, 2*M*N ) 460 END IF 461C 462C Test the input scalar arguments. 463C 464 IF( .NOT.LREDUC .AND. .NOT.LSAME( REDUCE, 'N' ) ) THEN 465 INFO = -1 466 ELSE IF( .NOT.LTRANN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 467 INFO = -2 468 ELSE IF( LTRANN ) THEN 469 IF( .NOT.LJOBDF .AND. .NOT.LSAME( JOBD, 'N' ) ) 470 $ INFO = -3 471 END IF 472 IF( M.LT.0 ) THEN 473 INFO = -4 474 ELSE IF( N.LT.0 ) THEN 475 INFO = -5 476 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 477 INFO = -7 478 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 479 INFO = -9 480 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 481 INFO = -11 482 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 483 INFO = -13 484 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 485 INFO = -15 486 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 487 INFO = -17 488 ELSE IF( ( .NOT.LREDRA .AND. LDP.LT.1 ) .OR. 489 $ ( LREDRA .AND. LDP.LT.MAX( 1, M ) ) ) THEN 490 INFO = -21 491 ELSE IF( ( .NOT.LREDRA .AND. LDQ.LT.1 ) .OR. 492 $ ( LREDRA .AND. LDQ.LT.MAX( 1, M ) ) ) THEN 493 INFO = -23 494 ELSE IF( ( .NOT.LREDRB .AND. LDU.LT.1 ) .OR. 495 $ ( LREDRB .AND. LDU.LT.MAX( 1, N ) ) ) THEN 496 INFO = -25 497 ELSE IF( ( .NOT.LREDRB .AND. LDV.LT.1 ) .OR. 498 $ ( LREDRB .AND. LDV.LT.MAX( 1, N ) ) ) THEN 499 INFO = -27 500 ELSE IF( LDWORK.LT.MINWRK ) THEN 501 INFO = -30 502 END IF 503C 504 IF ( INFO.NE.0 ) THEN 505C 506C Error return. 507C 508 CALL XERBLA( 'SB04OD', -INFO ) 509 RETURN 510 END IF 511C 512C Quick return if possible. 513C 514 IF ( N.EQ.0 .OR. M.EQ.0 ) THEN 515 SCALE = ONE 516 DWORK(1) = ONE 517 IF ( LTRANN ) THEN 518 IF ( LJOBDF ) DIF = ONE 519 END IF 520 RETURN 521 END IF 522 WRKOPT = 1 523 SUFWRK = LDWORK.GE.M*N 524C 525C STEP 1: Reduce (A,D) and/or (B,E) to generalized Schur form. 526C 527C (Note: Comments in the code beginning "Workspace:" describe the 528C minimal amount of real workspace needed at that point in the 529C code, as well as the preferred amount for good performance. 530C NB refers to the optimal block size for the immediately 531C following subroutine, as returned by ILAENV.) 532C 533 IF ( LREDUC ) THEN 534C 535C Get machine constants. 536C 537 SAFMIN = DLAMCH( 'Safe minimum' ) 538 SAFMAX = ONE / SAFMIN 539 CALL DLABAD( SAFMIN, SAFMAX ) 540 SMLNUM = SQRT( SAFMIN ) / DLAMCH( 'Precision' ) 541 BIGNUM = ONE / SMLNUM 542C 543 IF ( .NOT.LREDUB ) THEN 544C 545C Scale A if max element outside range [SMLNUM,BIGNUM]. 546C 547 ANRM = DLANGE( 'M', M, M, A, LDA, DWORK ) 548 ILASCL = .FALSE. 549 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 550 ANRMTO = SMLNUM 551 ILASCL = .TRUE. 552 ELSE IF( ANRM.GT.BIGNUM ) THEN 553 ANRMTO = BIGNUM 554 ILASCL = .TRUE. 555 END IF 556 IF( ILASCL ) 557 $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, M, M, A, LDA, 558 $ IERR ) 559C 560C Scale D if max element outside range [SMLNUM,BIGNUM] 561C 562 DNRM = DLANGE( 'M', M, M, D, LDD, DWORK ) 563 ILDSCL = .FALSE. 564 IF( DNRM.GT.ZERO .AND. DNRM.LT.SMLNUM ) THEN 565 DNRMTO = SMLNUM 566 ILDSCL = .TRUE. 567 ELSE IF( DNRM.GT.BIGNUM ) THEN 568 DNRMTO = BIGNUM 569 ILDSCL = .TRUE. 570 END IF 571 IF( ILDSCL ) 572 $ CALL DLASCL( 'G', 0, 0, DNRM, DNRMTO, M, M, D, LDD, 573 $ IERR ) 574C 575C Reduce (A,D) to generalized Schur form. 576C Workspace: need 7*M; 577C prefer 5*M + M*(NB+1). 578C 579 CALL DGGES( 'Vectors left', 'Vectors right', 'N', 0, N, 580 $ A, LDA, D, LDD, SDIM, 581 $ DWORK, DWORK(M+1), DWORK(2*M+1), 582 $ P, LDP, Q, LDQ, 583 $ DWORK(3*M+1), LDWORK-3*M, 584 $ 0, INFO ) 585 586C 587C Undo scaling 588C 589 IF( ILASCL ) 590 $ CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, M, M, A, LDA, 591 $ IERR ) 592C 593 IF( ILDSCL ) 594 $ CALL DLASCL( 'U', 0, 0, DNRMTO, DNRM, M, M, D, LDD, 595 $ IERR ) 596C 597 IF ( INFO.NE.0 ) THEN 598 INFO = 1 599 RETURN 600 END IF 601 WRKOPT = MAX( WRKOPT, INT( DWORK(3*M+1) ) + 3*M ) 602 END IF 603 IF ( .NOT.LREDUA ) THEN 604C 605C Scale B if max element outside range [SMLNUM,BIGNUM] 606C 607 BNRM = DLANGE( 'M', N, N, B, LDB, DWORK ) 608 ILBSCL = .FALSE. 609 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 610 BNRMTO = SMLNUM 611 ILBSCL = .TRUE. 612 ELSE IF( BNRM.GT.BIGNUM ) THEN 613 BNRMTO = BIGNUM 614 ILBSCL = .TRUE. 615 END IF 616 IF( ILBSCL ) 617 $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, 618 $ IERR ) 619C 620C Scale E if max element outside range [SMLNUM,BIGNUM] 621C 622 ENRM = DLANGE( 'M', N, N, E, LDE, DWORK ) 623 ILESCL = .FALSE. 624 IF( ENRM.GT.ZERO .AND. ENRM.LT.SMLNUM ) THEN 625 ENRMTO = SMLNUM 626 ILESCL = .TRUE. 627 ELSE IF( ENRM.GT.BIGNUM ) THEN 628 ENRMTO = BIGNUM 629 ILESCL = .TRUE. 630 END IF 631 IF( ILESCL ) 632 $ CALL DLASCL( 'G', 0, 0, ENRM, ENRMTO, N, N, E, LDE, 633 $ IERR ) 634C 635C Reduce (B,E) to generalized Schur form. 636C Workspace: need 7*N; 637C prefer 5*N + N*(NB+1). 638C 639 CALL DGGES( 'Vectors left', 'Vectors right', 'N', 0, N, 640 $ B, LDB, E, LDE, SDIM, 641 $ DWORK, DWORK(N+1), DWORK(2*N+1), 642 $ U, LDU, V, LDV, 643 $ DWORK(3*N+1), LDWORK-3*N, 644 $ 0, INFO ) 645C 646C Undo scaling 647C 648 IF( ILBSCL ) 649 $ CALL DLASCL( 'H', 0, 0, BNRMTO, BNRM, N, N, B, LDB, 650 $ IERR ) 651C 652 IF( ILESCL ) 653 $ CALL DLASCL( 'U', 0, 0, ENRMTO, ENRM, N, N, E, LDE, 654 $ IERR ) 655C 656 IF ( INFO.NE.0 ) THEN 657 INFO = 1 658 RETURN 659 END IF 660 WRKOPT = MAX( WRKOPT, INT( DWORK(3*N+1) ) + 3*N ) 661 END IF 662 END IF 663C 664 IF (.NOT.LREDUR ) THEN 665C 666C Set INFO = 2 if A and/or B are/is not in quasi-triangular form. 667C 668 IF (.NOT.LREDUA ) THEN 669 I = 1 670C 671 20 CONTINUE 672 IF ( I.LE.M-2 ) THEN 673 IF ( A(I+1,I).NE.ZERO ) THEN 674 IF ( A(I+2,I+1).NE.ZERO ) THEN 675 INFO = 2 676 RETURN 677 ELSE 678 I = I + 1 679 END IF 680 END IF 681 I = I + 1 682 GO TO 20 683 END IF 684 END IF 685C 686 IF (.NOT.LREDUB ) THEN 687 I = 1 688C 689 40 CONTINUE 690 IF ( I.LE.N-2 ) THEN 691 IF ( B(I+1,I).NE.ZERO ) THEN 692 IF ( B(I+2,I+1).NE.ZERO ) THEN 693 INFO = 2 694 RETURN 695 ELSE 696 I = I + 1 697 END IF 698 END IF 699 I = I + 1 700 GO TO 40 701 END IF 702 END IF 703 END IF 704C 705C STEP 2: Modify right hand sides (C,F). 706C 707 IF ( LREDUC ) THEN 708 WRKOPT = MAX( WRKOPT, M*N ) 709 IF ( SUFWRK ) THEN 710C 711C Enough workspace for a BLAS 3 calculation. 712C 713 IF ( LTRANN ) THEN 714C 715C Equation (1). 716C 717 IF ( .NOT.LREDUB ) THEN 718 CALL DGEMM( 'Transpose', 'No transpose', M, N, M, ONE, 719 $ P, LDP, C, LDC, ZERO, DWORK, M ) 720 ELSE 721 CALL DLACPY( 'Full', M, N, C, LDC, DWORK, M ) 722 END IF 723 IF ( .NOT.LREDUA ) THEN 724 CALL DGEMM( 'No transpose', 'No transpose', M, N, N, 725 $ ONE, DWORK, M, V, LDV, ZERO, C, LDC ) 726 ELSE 727 CALL DLACPY( 'Full', M, N, DWORK, M, C, LDC ) 728 END IF 729 IF ( .NOT.LREDUB ) THEN 730 CALL DGEMM( 'Transpose', 'No transpose', M, N, M, ONE, 731 $ P, LDP, F, LDF, ZERO, DWORK, M ) 732 ELSE 733 CALL DLACPY( 'Full', M, N, F, LDF, DWORK, M ) 734 END IF 735 IF ( .NOT.LREDUA ) THEN 736 CALL DGEMM( 'No transpose', 'No transpose', M, N, N, 737 $ ONE, DWORK, M, V, LDV, ZERO, F, LDF ) 738 ELSE 739 CALL DLACPY( 'Full', M, N, DWORK, M, F, LDF ) 740 END IF 741 ELSE 742C 743C Equation (2). 744C 745 IF ( .NOT.LREDUB ) THEN 746 CALL DGEMM( 'Transpose', 'No transpose', M, N, M, ONE, 747 $ Q, LDQ, C, LDC, ZERO, DWORK, M ) 748 ELSE 749 CALL DLACPY( 'Full', M, N, C, LDC, DWORK, M ) 750 END IF 751 IF ( .NOT.LREDUA ) THEN 752 CALL DGEMM( 'No transpose', 'No transpose', M, N, N, 753 $ ONE, DWORK, M, V, LDV, ZERO, C, LDC ) 754 ELSE 755 CALL DLACPY( 'Full', M, N, DWORK, M, C, LDC ) 756 END IF 757 IF ( .NOT.LREDUB ) THEN 758 CALL DGEMM( 'Transpose', 'No transpose', M, N, M, ONE, 759 $ P, LDP, F, LDF, ZERO, DWORK, M ) 760 ELSE 761 CALL DLACPY( 'Full', M, N, F, LDF, DWORK, M ) 762 END IF 763 IF ( .NOT.LREDUA ) THEN 764 CALL DGEMM( 'No transpose', 'No transpose', M, N, N, 765 $ ONE, DWORK, M, U, LDU, ZERO, F, LDF ) 766 ELSE 767 CALL DLACPY( 'Full', M, N, DWORK, M, F, LDF ) 768 END IF 769 END IF 770 ELSE 771C 772C Use a BLAS 2 calculation. 773C 774 IF ( LTRANN ) THEN 775C 776C Equation (1). 777C 778 IF ( .NOT.LREDUB ) THEN 779C 780 DO 60 I = 1, N 781 CALL DGEMV( 'Transpose', M, M, ONE, P, LDP, C(1,I), 782 $ 1, ZERO, DWORK, 1 ) 783 CALL DCOPY( M, DWORK, 1, C(1,I), 1 ) 784 60 CONTINUE 785C 786 END IF 787 IF ( .NOT.LREDUA ) THEN 788C 789 DO 80 I = 1, M 790 CALL DGEMV( 'Transpose', N, N, ONE, V, LDV, C(I,1), 791 $ LDC, ZERO, DWORK, 1 ) 792 CALL DCOPY( N, DWORK, 1, C(I,1), LDC ) 793 80 CONTINUE 794C 795 END IF 796 IF ( .NOT.LREDUB ) THEN 797C 798 DO 100 I = 1, N 799 CALL DGEMV( 'Transpose', M, M, ONE, P, LDP, F(1,I), 800 $ 1, ZERO, DWORK, 1 ) 801 CALL DCOPY( M, DWORK, 1, F(1,I), 1 ) 802 100 CONTINUE 803C 804 END IF 805 IF ( .NOT.LREDUA ) THEN 806C 807 DO 120 I = 1, M 808 CALL DGEMV( 'Transpose', N, N, ONE, V, LDV, F(I,1), 809 $ LDF, ZERO, DWORK, 1 ) 810 CALL DCOPY( N, DWORK, 1, F(I,1), LDF ) 811 120 CONTINUE 812C 813 END IF 814 ELSE 815C 816C Equation (2). 817C 818 IF ( .NOT.LREDUB ) THEN 819C 820 DO 140 I = 1, N 821 CALL DGEMV( 'Transpose', M, M, ONE, Q, LDQ, C(1,I), 822 $ 1, ZERO, DWORK, 1 ) 823 CALL DCOPY( M, DWORK, 1, C(1,I), 1 ) 824 140 CONTINUE 825C 826 END IF 827 IF ( .NOT.LREDUA ) THEN 828C 829 DO 160 I = 1, M 830 CALL DGEMV( 'Transpose', N, N, ONE, V, LDV, C(I,1), 831 $ LDC, ZERO, DWORK, 1 ) 832 CALL DCOPY( N, DWORK, 1, C(I,1), LDC ) 833 160 CONTINUE 834C 835 END IF 836 IF ( .NOT.LREDUB ) THEN 837C 838 DO 180 I = 1, N 839 CALL DGEMV( 'Transpose', M, M, ONE, P, LDP, F(1,I), 840 $ 1, ZERO, DWORK, 1 ) 841 CALL DCOPY( M, DWORK, 1, F(1,I), 1 ) 842 180 CONTINUE 843C 844 END IF 845 IF ( .NOT.LREDUA ) THEN 846C 847 DO 200 I = 1, M 848 CALL DGEMV( 'Transpose', N, N, ONE, U, LDU, F(I,1), 849 $ LDF, ZERO, DWORK, 1 ) 850 CALL DCOPY( N, DWORK, 1, F(I,1), LDF ) 851 200 CONTINUE 852C 853 END IF 854 END IF 855 END IF 856 END IF 857C 858C STEP 3: Solve the transformed system and compute the Dif 859C estimator. 860C 861 IF ( LTRANN ) THEN 862 IF ( LJOBD ) THEN 863 IJOB = 1 864 ELSE IF ( LJOBF ) THEN 865 IJOB = 2 866 ELSE IF ( LJOB1 ) THEN 867 IJOB = 3 868 ELSE IF ( LJOB2 ) THEN 869 IJOB = 4 870 ELSE 871 IJOB = 0 872 END IF 873 ELSE 874 IJOB = 0 875 END IF 876C 877C Workspace: need 2*M*N if TRANS = 'N' and JOBD = 'D' or 'F'; 878C 1, otherwise. 879C 880 CALL DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, 881 $ E, LDE, F, LDF, SCALE, DIF, DWORK, LDWORK, IWORK, 882 $ INFO ) 883 IF ( INFO.NE.0 ) THEN 884 INFO = 3 885 RETURN 886 END IF 887 IF ( LTRANN ) THEN 888 IF ( LJOBD.OR.LJOBF ) 889 $ WRKOPT = MAX( WRKOPT, 2*M*N ) 890 END IF 891C 892C STEP 4: Back transformation of the solution. 893C 894 IF ( LREDUC ) THEN 895 IF (SUFWRK ) THEN 896C 897C Enough workspace for a BLAS 3 calculation. 898C 899 IF ( LTRANN ) THEN 900C 901C Equation (1). 902C 903 IF ( .NOT.LREDUB ) THEN 904 CALL DGEMM( 'No transpose', 'No transpose', M, N, M, 905 $ ONE, Q, LDQ, C, LDC, ZERO, DWORK, M ) 906 ELSE 907 CALL DLACPY( 'Full', M, N, C, LDC, DWORK, M ) 908 END IF 909 IF ( .NOT.LREDUA ) THEN 910 CALL DGEMM( 'No transpose', 'Transpose', M, N, N, ONE, 911 $ DWORK, M, V, LDV, ZERO, C, LDC ) 912 ELSE 913 CALL DLACPY( 'Full', M, N, DWORK, M, C, LDC ) 914 END IF 915 IF ( .NOT.LREDUB ) THEN 916 CALL DGEMM( 'No transpose', 'No transpose', M, N, M, 917 $ ONE, P, LDP, F, LDF, ZERO, DWORK, M ) 918 ELSE 919 CALL DLACPY( 'Full', M, N, F, LDF, DWORK, M ) 920 END IF 921 IF ( .NOT.LREDUA ) THEN 922 CALL DGEMM( 'No transpose', 'Transpose', M, N, N, ONE, 923 $ DWORK, M, U, LDU, ZERO, F, LDF ) 924 ELSE 925 CALL DLACPY( 'Full', M, N, DWORK, M, F, LDF ) 926 END IF 927 ELSE 928C 929C Equation (2). 930C 931 IF ( .NOT.LREDUB ) THEN 932 CALL DGEMM( 'No transpose', 'No transpose', M, N, M, 933 $ ONE, P, LDP, C, LDC, ZERO, DWORK, M ) 934 ELSE 935 CALL DLACPY( 'Full', M, N, C, LDC, DWORK, M ) 936 END IF 937 IF ( .NOT.LREDUA ) THEN 938 CALL DGEMM( 'No transpose', 'Transpose', M, N, N, 939 $ ONE, DWORK, M, V, LDV, ZERO, C, LDC ) 940 ELSE 941 CALL DLACPY( 'Full', M, N, DWORK, M, C, LDC ) 942 END IF 943 IF ( .NOT.LREDUB ) THEN 944 CALL DGEMM( 'No transpose', 'No transpose', M, N, M, 945 $ ONE, P, LDP, F, LDF, ZERO, DWORK, M ) 946 ELSE 947 CALL DLACPY( 'Full', M, N, F, LDF, DWORK, M ) 948 END IF 949 IF ( .NOT.LREDUA ) THEN 950 CALL DGEMM( 'No transpose', 'Transpose', M, N, N, 951 $ ONE, DWORK, M, V, LDV, ZERO, F, LDF ) 952 ELSE 953 CALL DLACPY( 'Full', M, N, DWORK, M, F, LDF ) 954 END IF 955 END IF 956 ELSE 957C 958C Use a BLAS 2 calculation. 959C 960 IF ( LTRANN ) THEN 961C 962C Equation (1). 963C 964 IF ( .NOT.LREDUB ) THEN 965C 966 DO 220 I = 1, N 967 CALL DGEMV( 'No transpose', M, M, ONE, Q, LDQ, 968 $ C(1,I), 1, ZERO, DWORK, 1 ) 969 CALL DCOPY( M, DWORK, 1, C(1,I), 1 ) 970 220 CONTINUE 971C 972 END IF 973 IF ( .NOT.LREDUA ) THEN 974C 975 DO 240 I = 1, M 976 CALL DGEMV( 'No transpose', N, N, ONE, V, LDV, 977 $ C(I,1), LDC, ZERO, DWORK, 1 ) 978 CALL DCOPY( N, DWORK, 1, C(I,1), LDC ) 979 240 CONTINUE 980C 981 END IF 982 IF ( .NOT.LREDUB ) THEN 983C 984 DO 260 I = 1, N 985 CALL DGEMV( 'No transpose', M, M, ONE, P, LDP, 986 $ F(1,I), 1, ZERO, DWORK, 1 ) 987 CALL DCOPY( M, DWORK, 1, F(1,I), 1 ) 988 260 CONTINUE 989C 990 END IF 991 IF ( .NOT.LREDUA ) THEN 992C 993 DO 280 I = 1, M 994 CALL DGEMV( 'No transpose', N, N, ONE, U, LDU, 995 $ F(I,1), LDF, ZERO, DWORK, 1 ) 996 CALL DCOPY( N, DWORK, 1, F(I,1), LDF ) 997 280 CONTINUE 998C 999 END IF 1000 ELSE 1001C 1002C Equation (2). 1003C 1004 IF ( .NOT.LREDUB ) THEN 1005C 1006 DO 300 I = 1, N 1007 CALL DGEMV( 'No transpose', M, M, ONE, P, LDP, 1008 $ C(1,I), 1, ZERO, DWORK, 1 ) 1009 CALL DCOPY( M, DWORK, 1, C(1,I), 1 ) 1010 300 CONTINUE 1011C 1012 END IF 1013 IF ( .NOT.LREDUA ) THEN 1014C 1015 DO 320 I = 1, M 1016 CALL DGEMV( 'No transpose', N, N, ONE, V, LDV, 1017 $ C(I,1), LDC, ZERO, DWORK, 1 ) 1018 CALL DCOPY( N, DWORK, 1, C(I,1), LDC ) 1019 320 CONTINUE 1020C 1021 END IF 1022 IF ( .NOT.LREDUB ) THEN 1023C 1024 DO 340 I = 1, N 1025 CALL DGEMV( 'No transpose', M, M, ONE, P, LDP, 1026 $ F(1,I), 1, ZERO, DWORK, 1 ) 1027 CALL DCOPY( M, DWORK, 1, F(1,I), 1 ) 1028 340 CONTINUE 1029C 1030 END IF 1031 IF ( .NOT.LREDUA ) THEN 1032C 1033 DO 360 I = 1, M 1034 CALL DGEMV( 'No transpose', N, N, ONE, V, LDV, 1035 $ F(I,1), LDF, ZERO, DWORK, 1 ) 1036 CALL DCOPY( N, DWORK, 1, F(I,1), LDF ) 1037 360 CONTINUE 1038C 1039 END IF 1040 END IF 1041 END IF 1042 END IF 1043C 1044 DWORK(1) = WRKOPT 1045C 1046 RETURN 1047C *** Last line of SB04OD *** 1048 END 1049