1 SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q, 2 $ LDQ, Z, LDZ, B, LDB, SCALE, ALPHAR, ALPHAI, 3 $ BETA, DWORK, LDWORK, INFO ) 4C 5C WARNING 6C 7C This file alters the SLICOT routine SG03BD. 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 DGEGS 15C is replaced by DGGES. 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 compute the Cholesky factor U of the matrix X, 39C 40C T 41C X = op(U) * op(U), 42C 43C which is the solution of either the generalized 44C c-stable continuous-time Lyapunov equation 45C 46C T T 47C op(A) * X * op(E) + op(E) * X * op(A) 48C 49C 2 T 50C = - SCALE * op(B) * op(B), (1) 51C 52C or the generalized d-stable discrete-time Lyapunov equation 53C 54C T T 55C op(A) * X * op(A) - op(E) * X * op(E) 56C 57C 2 T 58C = - SCALE * op(B) * op(B), (2) 59C 60C without first finding X and without the need to form the matrix 61C op(B)**T * op(B). 62C 63C op(K) is either K or K**T for K = A, B, E, U. A and E are N-by-N 64C matrices, op(B) is an M-by-N matrix. The resulting matrix U is an 65C N-by-N upper triangular matrix with non-negative entries on its 66C main diagonal. SCALE is an output scale factor set to avoid 67C overflow in U. 68C 69C In the continuous-time case (1) the pencil A - lambda * E must be 70C c-stable (that is, all eigenvalues must have negative real parts). 71C In the discrete-time case (2) the pencil A - lambda * E must be 72C d-stable (that is, the moduli of all eigenvalues must be smaller 73C than one). 74C 75C ARGUMENTS 76C 77C Mode Parameters 78C 79C DICO CHARACTER*1 80C Specifies which type of the equation is considered: 81C = 'C': Continuous-time equation (1); 82C = 'D': Discrete-time equation (2). 83C 84C FACT CHARACTER*1 85C Specifies whether the generalized real Schur 86C factorization of the pencil A - lambda * E is supplied 87C on entry or not: 88C = 'N': Factorization is not supplied; 89C = 'F': Factorization is supplied. 90C 91C TRANS CHARACTER*1 92C Specifies whether the transposed equation is to be solved 93C or not: 94C = 'N': op(A) = A, op(E) = E; 95C = 'T': op(A) = A**T, op(E) = E**T. 96C 97C Input/Output Parameters 98C 99C N (input) INTEGER 100C The order of the matrix A. N >= 0. 101C 102C M (input) INTEGER 103C The number of rows in the matrix op(B). M >= 0. 104C 105C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 106C On entry, if FACT = 'F', then the leading N-by-N upper 107C Hessenberg part of this array must contain the 108C generalized Schur factor A_s of the matrix A (see 109C definition (3) in section METHOD). A_s must be an upper 110C quasitriangular matrix. The elements below the upper 111C Hessenberg part of the array A are not referenced. 112C If FACT = 'N', then the leading N-by-N part of this 113C array must contain the matrix A. 114C On exit, the leading N-by-N part of this array contains 115C the generalized Schur factor A_s of the matrix A. (A_s is 116C an upper quasitriangular matrix.) 117C 118C LDA INTEGER 119C The leading dimension of the array A. LDA >= MAX(1,N). 120C 121C E (input/output) DOUBLE PRECISION array, dimension (LDE,N) 122C On entry, if FACT = 'F', then the leading N-by-N upper 123C triangular part of this array must contain the 124C generalized Schur factor E_s of the matrix E (see 125C definition (4) in section METHOD). The elements below the 126C upper triangular part of the array E are not referenced. 127C If FACT = 'N', then the leading N-by-N part of this 128C array must contain the coefficient matrix E of the 129C equation. 130C On exit, the leading N-by-N part of this array contains 131C the generalized Schur factor E_s of the matrix E. (E_s is 132C an upper triangular matrix.) 133C 134C LDE INTEGER 135C The leading dimension of the array E. LDE >= MAX(1,N). 136C 137C Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 138C On entry, if FACT = 'F', then the leading N-by-N part of 139C this array must contain the orthogonal matrix Q from 140C the generalized Schur factorization (see definitions (3) 141C and (4) in section METHOD). 142C If FACT = 'N', Q need not be set on entry. 143C On exit, the leading N-by-N part of this array contains 144C the orthogonal matrix Q from the generalized Schur 145C factorization. 146C 147C LDQ INTEGER 148C The leading dimension of the array Q. LDQ >= MAX(1,N). 149C 150C Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) 151C On entry, if FACT = 'F', then the leading N-by-N part of 152C this array must contain the orthogonal matrix Z from 153C the generalized Schur factorization (see definitions (3) 154C and (4) in section METHOD). 155C If FACT = 'N', Z need not be set on entry. 156C On exit, the leading N-by-N part of this array contains 157C the orthogonal matrix Z from the generalized Schur 158C factorization. 159C 160C LDZ INTEGER 161C The leading dimension of the array Z. LDZ >= MAX(1,N). 162C 163C B (input/output) DOUBLE PRECISION array, dimension (LDB,N1) 164C On entry, if TRANS = 'T', the leading N-by-M part of this 165C array must contain the matrix B and N1 >= MAX(M,N). 166C If TRANS = 'N', the leading M-by-N part of this array 167C must contain the matrix B and N1 >= N. 168C On exit, the leading N-by-N part of this array contains 169C the Cholesky factor U of the solution matrix X of the 170C problem, X = op(U)**T * op(U). 171C If M = 0 and N > 0, then U is set to zero. 172C 173C LDB INTEGER 174C The leading dimension of the array B. 175C If TRANS = 'T', LDB >= MAX(1,N). 176C If TRANS = 'N', LDB >= MAX(1,M,N). 177C 178C SCALE (output) DOUBLE PRECISION 179C The scale factor set to avoid overflow in U. 180C 0 < SCALE <= 1. 181C 182C ALPHAR (output) DOUBLE PRECISION array, dimension (N) 183C ALPHAI (output) DOUBLE PRECISION array, dimension (N) 184C BETA (output) DOUBLE PRECISION array, dimension (N) 185C If INFO = 0, 3, 5, 6, or 7, then 186C (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, are the 187C eigenvalues of the matrix pencil A - lambda * E. 188C 189C Workspace 190C 191C DWORK DOUBLE PRECISION array, dimension (LDWORK) 192C On exit, if INFO = 0, DWORK(1) returns the optimal value 193C of LDWORK. 194C 195C LDWORK INTEGER 196C The dimension of the array DWORK. 197C LDWORK >= MAX(1,4*N,6*N-6), if FACT = 'N'; 198C LDWORK >= MAX(1,2*N,6*N-6), if FACT = 'F'. 199C For good performance, LDWORK should be larger. 200C 201C Error indicator 202C 203C INFO INTEGER 204C = 0: successful exit; 205C < 0: if INFO = -i, the i-th argument had an illegal 206C value; 207C = 1: the pencil A - lambda * E is (nearly) singular; 208C perturbed values were used to solve the equation 209C (but the reduced (quasi)triangular matrices A and E 210C are unchanged); 211C = 2: FACT = 'F' and the matrix contained in the upper 212C Hessenberg part of the array A is not in upper 213C quasitriangular form; 214C = 3: FACT = 'F' and there is a 2-by-2 block on the main 215C diagonal of the pencil A_s - lambda * E_s whose 216C eigenvalues are not conjugate complex; 217C = 4: FACT = 'N' and the pencil A - lambda * E cannot be 218C reduced to generalized Schur form: LAPACK routine 219C DGEGS has failed to converge; 220C = 5: DICO = 'C' and the pencil A - lambda * E is not 221C c-stable; 222C = 6: DICO = 'D' and the pencil A - lambda * E is not 223C d-stable; 224C = 7: the LAPACK routine DSYEVX utilized to factorize M3 225C failed to converge in the discrete-time case (see 226C section METHOD for SLICOT Library routine SG03BU). 227C This error is unlikely to occur. 228C 229C METHOD 230C 231C An extension [2] of Hammarling's method [1] to generalized 232C Lyapunov equations is utilized to solve (1) or (2). 233C 234C First the pencil A - lambda * E is reduced to real generalized 235C Schur form A_s - lambda * E_s by means of orthogonal 236C transformations (QZ-algorithm): 237C 238C A_s = Q**T * A * Z (upper quasitriangular) (3) 239C 240C E_s = Q**T * E * Z (upper triangular). (4) 241C 242C If the pencil A - lambda * E has already been factorized prior to 243C calling the routine however, then the factors A_s, E_s, Q and Z 244C may be supplied and the initial factorization omitted. 245C 246C Depending on the parameters TRANS and M the N-by-N upper 247C triangular matrix B_s is defined as follows. In any case Q_B is 248C an M-by-M orthogonal matrix, which need not be accumulated. 249C 250C 1. If TRANS = 'N' and M < N, B_s is the upper triangular matrix 251C from the QR-factorization 252C 253C ( Q_B O ) ( B * Z ) 254C ( ) * B_s = ( ), 255C ( O I ) ( O ) 256C 257C where the O's are zero matrices of proper size and I is the 258C identity matrix of order N-M. 259C 260C 2. If TRANS = 'N' and M >= N, B_s is the upper triangular matrix 261C from the (rectangular) QR-factorization 262C 263C ( B_s ) 264C Q_B * ( ) = B * Z, 265C ( O ) 266C 267C where O is the (M-N)-by-N zero matrix. 268C 269C 3. If TRANS = 'T' and M < N, B_s is the upper triangular matrix 270C from the RQ-factorization 271C 272C ( Q_B O ) 273C (B_s O ) * ( ) = ( Q**T * B O ). 274C ( O I ) 275C 276C 4. If TRANS = 'T' and M >= N, B_s is the upper triangular matrix 277C from the (rectangular) RQ-factorization 278C 279C ( B_s O ) * Q_B = Q**T * B, 280C 281C where O is the N-by-(M-N) zero matrix. 282C 283C Assuming SCALE = 1, the transformation of A, E and B described 284C above leads to the reduced continuous-time equation 285C 286C T T 287C op(A_s) op(U_s) op(U_s) op(E_s) 288C 289C T T 290C + op(E_s) op(U_s) op(U_s) op(A_s) 291C 292C T 293C = - op(B_s) op(B_s) (5) 294C 295C or to the reduced discrete-time equation 296C 297C T T 298C op(A_s) op(U_s) op(U_s) op(A_s) 299C 300C T T 301C - op(E_s) op(U_s) op(U_s) op(E_s) 302C 303C T 304C = - op(B_s) op(B_s). (6) 305C 306C For brevity we restrict ourself to equation (5) and the case 307C TRANS = 'N'. The other three cases can be treated in a similar 308C fashion. 309C 310C We use the following partitioning for the matrices A_s, E_s, B_s 311C and U_s 312C 313C ( A11 A12 ) ( E11 E12 ) 314C A_s = ( ), E_s = ( ), 315C ( 0 A22 ) ( 0 E22 ) 316C 317C ( B11 B12 ) ( U11 U12 ) 318C B_s = ( ), U_s = ( ). (7) 319C ( 0 B22 ) ( 0 U22 ) 320C 321C The size of the (1,1)-blocks is 1-by-1 (iff A_s(2,1) = 0.0) or 322C 2-by-2. 323C 324C We compute U11 and U12**T in three steps. 325C 326C Step I: 327C 328C From (5) and (7) we get the 1-by-1 or 2-by-2 equation 329C 330C T T T T 331C A11 * U11 * U11 * E11 + E11 * U11 * U11 * A11 332C 333C T 334C = - B11 * B11. 335C 336C For brevity, details are omitted here. See [2]. The technique 337C for computing U11 is similar to those applied to standard 338C Lyapunov equations in Hammarling's algorithm ([1], section 6). 339C 340C Furthermore, the auxiliary matrices M1 and M2 defined as 341C follows 342C 343C -1 -1 344C M1 = U11 * A11 * E11 * U11 345C 346C -1 -1 347C M2 = B11 * E11 * U11 348C 349C are computed in a numerically reliable way. 350C 351C Step II: 352C 353C The generalized Sylvester equation 354C 355C T T T T 356C A22 * U12 + E22 * U12 * M1 = 357C 358C T T T T T 359C - B12 * M2 - A12 * U11 - E12 * U11 * M1 360C 361C is solved for U12**T. 362C 363C Step III: 364C 365C It can be shown that 366C 367C T T T T 368C A22 * U22 * U22 * E22 + E22 * U22 * U22 * A22 = 369C 370C T T 371C - B22 * B22 - y * y (8) 372C 373C holds, where y is defined as 374C 375C T T T T T T 376C y = B12 - ( E12 * U11 + E22 * U12 ) * M2 . 377C 378C If B22_tilde is the square triangular matrix arising from the 379C (rectangular) QR-factorization 380C 381C ( B22_tilde ) ( B22 ) 382C Q_B_tilde * ( ) = ( ), 383C ( O ) ( y**T ) 384C 385C where Q_B_tilde is an orthogonal matrix of order N, then 386C 387C T T T 388C - B22 * B22 - y * y = - B22_tilde * B22_tilde. 389C 390C Replacing the right hand side in (8) by the term 391C - B22_tilde**T * B22_tilde leads to a reduced generalized 392C Lyapunov equation of lower dimension compared to (5). 393C 394C The recursive application of the steps I to III yields the 395C solution U_s of the equation (5). 396C 397C It remains to compute the solution matrix U of the original 398C problem (1) or (2) from the matrix U_s. To this end we transform 399C the solution back (with respect to the transformation that led 400C from (1) to (5) (from (2) to (6)) and apply the QR-factorization 401C (RQ-factorization). The upper triangular solution matrix U is 402C obtained by 403C 404C Q_U * U = U_s * Q**T (if TRANS = 'N') 405C 406C or 407C 408C U * Q_U = Z * U_s (if TRANS = 'T') 409C 410C where Q_U is an N-by-N orthogonal matrix. Again, the orthogonal 411C matrix Q_U need not be accumulated. 412C 413C REFERENCES 414C 415C [1] Hammarling, S.J. 416C Numerical solution of the stable, non-negative definite 417C Lyapunov equation. 418C IMA J. Num. Anal., 2, pp. 303-323, 1982. 419C 420C [2] Penzl, T. 421C Numerical solution of generalized Lyapunov equations. 422C Advances in Comp. Math., vol. 8, pp. 33-48, 1998. 423C 424C NUMERICAL ASPECTS 425C 426C The number of flops required by the routine is given by the 427C following table. Note that we count a single floating point 428C arithmetic operation as one flop. 429C 430C | FACT = 'F' FACT = 'N' 431C ---------+-------------------------------------------------- 432C M <= N | (13*N**3+6*M*N**2 (211*N**3+6*M*N**2 433C | +6*M**2*N-2*M**3)/3 +6*M**2*N-2*M**3)/3 434C | 435C M > N | (11*N**3+12*M*N**2)/3 (209*N**3+12*M*N**2)/3 436C 437C FURTHER COMMENTS 438C 439C The Lyapunov equation may be very ill-conditioned. In particular, 440C if DICO = 'D' and the pencil A - lambda * E has a pair of almost 441C reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost 442C degenerate pair of eigenvalues, then the Lyapunov equation will be 443C ill-conditioned. Perturbed values were used to solve the equation. 444C A condition estimate can be obtained from the routine SG03AD. 445C When setting the error indicator INFO, the routine does not test 446C for near instability in the equation but only for exact 447C instability. 448C 449C CONTRIBUTOR 450C 451C T. Penzl, Technical University Chemnitz, Germany, Aug. 1998. 452C 453C REVISIONS 454C 455C Sep. 1998 (V. Sima). 456C May 1999 (V. Sima). 457C March 2002 (A. Varga). 458C Feb. 2004 (V. Sima). 459C 460C KEYWORDS 461C 462C Lyapunov equation 463C 464C ****************************************************************** 465C 466C .. Parameters .. 467 DOUBLE PRECISION MONE, ONE, TWO, ZERO 468 PARAMETER ( MONE = -1.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, 469 $ ZERO = 0.0D+0 ) 470C .. Scalar Arguments .. 471 DOUBLE PRECISION SCALE 472 INTEGER INFO, LDA, LDB, LDE, LDQ, LDWORK, LDZ, M, N 473 CHARACTER DICO, FACT, TRANS 474C .. Array Arguments .. 475 DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*), 476 $ BETA(*), DWORK(*), E(LDE,*), Q(LDQ,*), Z(LDZ,*) 477C .. Local Scalars .. 478 DOUBLE PRECISION S1, S2, SAFMIN, WI, WR1, WR2 479 INTEGER I, INFO1, MINMN, MINWRK, OPTWRK, SDIM 480 LOGICAL ISDISC, ISFACT, ISTRAN 481C .. Local Arrays .. 482 DOUBLE PRECISION E1(2,2) 483C .. External Functions .. 484 DOUBLE PRECISION DLAMCH, DLAPY2 485 LOGICAL LSAME 486 EXTERNAL DLAMCH, DLAPY2, LSAME 487C .. External Subroutines .. 488 EXTERNAL DCOPY, DGGES, DGEMM, DGEMV, DGEQRF, DGERQF, 489 $ DLACPY, DLAG2, DLASET, DSCAL, DTRMM, SG03BU, 490 $ SG03BV, XERBLA 491C .. Intrinsic Functions .. 492 INTRINSIC ABS, DBLE, INT, MAX, MIN, SIGN 493C .. Executable Statements .. 494C 495C Decode input parameters. 496C 497 ISDISC = LSAME( DICO, 'D' ) 498 ISFACT = LSAME( FACT, 'F' ) 499 ISTRAN = LSAME( TRANS, 'T' ) 500C 501C Compute minimal workspace. 502C 503 IF (ISFACT ) THEN 504 MINWRK = MAX( 1, 2*N, 6*N-6 ) 505 ELSE 506 MINWRK = MAX( 1, 8*N, 6*N+16 ) 507 END IF 508C 509C Check the scalar input parameters. 510C 511 IF ( .NOT.( ISDISC .OR. LSAME( DICO, 'C' ) ) ) THEN 512 INFO = -1 513 ELSEIF ( .NOT.( ISFACT .OR. LSAME( FACT, 'N' ) ) ) THEN 514 INFO = -2 515 ELSEIF ( .NOT.( ISTRAN .OR. LSAME( TRANS, 'N' ) ) ) THEN 516 INFO = -3 517 ELSEIF ( N .LT. 0 ) THEN 518 INFO = -4 519 ELSEIF ( M .LT. 0 ) THEN 520 INFO = -5 521 ELSEIF ( LDA .LT. MAX( 1, N ) ) THEN 522 INFO = -7 523 ELSEIF ( LDE .LT. MAX( 1, N ) ) THEN 524 INFO = -9 525 ELSEIF ( LDQ .LT. MAX( 1, N ) ) THEN 526 INFO = -11 527 ELSEIF ( LDZ .LT. MAX( 1, N ) ) THEN 528 INFO = -13 529 ELSEIF ( ( ISTRAN .AND. ( LDB .LT. MAX( 1, N ) ) ) .OR. 530 $ ( .NOT.ISTRAN .AND. ( LDB .LT. MAX( 1, M, N ) ) ) ) THEN 531 INFO = -15 532 ELSEIF ( LDWORK .LT. MINWRK ) THEN 533 INFO = -21 534 ELSE 535 INFO = 0 536 END IF 537 IF ( INFO .NE. 0 ) THEN 538 CALL XERBLA( 'SG03BD', -INFO ) 539 RETURN 540 END IF 541C 542 SCALE = ONE 543C 544C Quick return if possible. 545C 546 MINMN = MIN( M, N ) 547 IF ( MINMN .EQ. 0 ) THEN 548 IF ( N.GT.0 ) 549 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDB ) 550 DWORK(1) = ONE 551 RETURN 552 ENDIF 553C 554 IF ( ISFACT ) THEN 555C 556C Make sure the upper Hessenberg part of A is quasitriangular. 557C 558 DO 20 I = 1, N-2 559 IF ( A(I+1,I).NE.ZERO .AND. A(I+2,I+1).NE.ZERO ) THEN 560 INFO = 2 561 RETURN 562 END IF 563 20 CONTINUE 564 END IF 565C 566 IF ( .NOT.ISFACT ) THEN 567C 568C Reduce the pencil A - lambda * E to generalized Schur form. 569C 570C A := Q**T * A * Z (upper quasitriangular) 571C E := Q**T * E * Z (upper triangular) 572C 573C ( Workspace: >= MAX(1,4*N) ) 574C 575 CALL DGGES( 'Vectors', 'Vectors', 'N', 0, N, 576 $ A, LDA, E, LDE, SDIM, ALPHAR, ALPHAI, BETA, 577 $ Q, LDQ, Z, LDZ, DWORK, LDWORK, 0, INFO1) 578 IF ( INFO1 .NE. 0 ) THEN 579 INFO = 4 580 RETURN 581 END IF 582 OPTWRK = INT( DWORK(1) ) 583 ELSE 584 OPTWRK = MINWRK 585 END IF 586C 587 IF ( ISFACT ) THEN 588C 589C If the matrix pencil A - lambda * E has been in generalized 590C Schur form on entry, compute its eigenvalues. 591C 592 SAFMIN = DLAMCH( 'Safe minimum' ) 593 E1(2,1) = ZERO 594 I = 1 595C WHILE ( I .LE. N ) DO 596 30 IF ( I .LE. N ) THEN 597 IF ( ( I.EQ.N ) .OR. ( A(MIN( I+1, N ),I).EQ.ZERO ) ) THEN 598 ALPHAR(I) = A(I,I) 599 ALPHAI(I) = ZERO 600 BETA(I) = E(I,I) 601 I = I+1 602 ELSE 603 E1(1,1) = E(I,I) 604 E1(1,2) = E(I,I+1) 605 E1(2,2) = E(I+1,I+1) 606 CALL DLAG2( A(I,I), LDA, E1, 2, SAFMIN, S1, S2, WR1, WR2, 607 $ WI ) 608 IF ( WI .EQ. ZERO ) INFO = 3 609 ALPHAR(I) = WR1 610 ALPHAI(I) = WI 611 BETA(I) = S1 612 ALPHAR(I+1) = WR2 613 ALPHAI(I+1) = -WI 614 BETA(I+1) = S2 615 I = I+2 616 END IF 617 GOTO 30 618 END IF 619C END WHILE 30 620 IF ( INFO.NE.0 ) RETURN 621 END IF 622C 623C Check on the stability of the matrix pencil A - lambda * E. 624C 625 DO 40 I = 1, N 626 IF ( ISDISC ) THEN 627 IF ( DLAPY2( ALPHAR(I), ALPHAI(I) ) .GE. ABS( BETA(I) ) ) 628 $ THEN 629 INFO = 6 630 RETURN 631 END IF 632 ELSE 633 IF ( ( ALPHAR(I).EQ.ZERO ) .OR. ( BETA(I).EQ.ZERO ) .OR. 634 $ ( SIGN( ONE,ALPHAR(I) )*SIGN( ONE, BETA(I) ) .GE. ZERO) ) 635 $ THEN 636 INFO = 5 637 RETURN 638 END IF 639 END IF 640 40 CONTINUE 641C 642C Transformation of the right hand side. 643C 644C B := B * Z or B := Q**T * B 645C 646C Use BLAS 3 if there is enough workspace. Otherwise, use BLAS 2. 647C 648C ( Workspace: max(1,N) ) 649C 650 IF ( .NOT.ISTRAN ) THEN 651 IF ( LDWORK .GE. N*M ) THEN 652 CALL DGEMM( 'NoTranspose', 'NoTranspose', M, N, N, ONE, B, 653 $ LDB, Z, LDZ, ZERO, DWORK, M ) 654 CALL DLACPY( 'All', M, N, DWORK, M, B, LDB ) 655 ELSE 656 DO 60 I = 1, M 657 CALL DCOPY( N, B(I,1), LDB, DWORK, 1 ) 658 CALL DGEMV( 'Transpose', N, N, ONE, Z, LDZ, DWORK, 1, 659 $ ZERO, B(I,1), LDB ) 660 60 CONTINUE 661 END IF 662 IF ( M .LT. N ) 663 $ CALL DLASET( 'All', N-M, N, ZERO, ZERO, B(M+1,1), LDB ) 664 ELSE 665 IF ( LDWORK .GE. N*M ) THEN 666 CALL DLACPY( 'All', N, M, B, LDB, DWORK, N ) 667 CALL DGEMM( 'Transpose', 'NoTranspose', N, M, N, ONE, Q, 668 $ LDQ, DWORK, N, ZERO, B, LDB ) 669 ELSE 670 DO 80 I = 1, M 671 CALL DCOPY( N, B(1,I), 1, DWORK, 1 ) 672 CALL DGEMV( 'Transpose', N, N, ONE, Q, LDQ, DWORK, 1, 673 $ ZERO, B(1,I), 1 ) 674 80 CONTINUE 675 END IF 676 IF ( M .LT. N ) 677 $ CALL DLASET( 'All', N, N-M, ZERO, ZERO, B(1,M+1), LDB ) 678 END IF 679 OPTWRK = MAX( OPTWRK, N*M ) 680C 681C Overwrite B with the triangular matrix of its QR-factorization 682C or its RQ-factorization. 683C (The entries on the main diagonal are non-negative.) 684C 685C ( Workspace: >= max(1,2*N) ) 686C 687 IF ( .NOT.ISTRAN ) THEN 688 IF ( M .GE. 2 ) THEN 689 CALL DGEQRF( M, N, B, LDB, DWORK, DWORK(N+1), LDWORK-N, 690 $ INFO1 ) 691 CALL DLASET( 'Lower', MAX( M, N )-1, MIN( M, N ), ZERO, 692 $ ZERO, B(2,1), LDB ) 693 END IF 694 DO 100 I = 1, MINMN 695 IF ( B(I,I) .LT. ZERO ) 696 $ CALL DSCAL( N+1-I, MONE, B(I,I), LDB ) 697 100 CONTINUE 698 ELSE 699 IF ( M .GE. 2 ) THEN 700 CALL DGERQF( N, M, B, LDB, DWORK, DWORK(N+1), LDWORK-N, 701 $ INFO1 ) 702 IF ( N .GE. M ) THEN 703 CALL DLASET( 'Lower', M-1, M-1, ZERO, ZERO, B(N-M+2,1), 704 $ LDB ) 705 IF ( N .GT. M ) THEN 706 DO 120 I = M, 1, -1 707 CALL DCOPY( N, B(1,I), 1, B(1,I+N-M), 1 ) 708 120 CONTINUE 709 CALL DLASET( 'All', N, N-M, ZERO, ZERO, B(1,1), LDB ) 710 END IF 711 ELSE 712 IF ( N .GT. 1 ) 713 $ CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, 714 $ B(2,M-N+1), LDB ) 715 DO 140 I = 1, N 716 CALL DCOPY( N, B(1,M-N+I), 1, B(1,I), 1 ) 717 140 CONTINUE 718 CALL DLASET( 'All', N, M-N, ZERO, ZERO, B(1,N+1), LDB ) 719 END IF 720 ELSE 721 IF ( N .NE. 1 ) THEN 722 CALL DCOPY( N, B(1,1), 1, B(1,N), 1 ) 723 CALL DLASET( 'All', N, 1, ZERO, ZERO, B(1,1), LDB ) 724 END IF 725 END IF 726 DO 160 I = N - MINMN + 1, N 727 IF ( B(I,I) .LT. ZERO ) 728 $ CALL DSCAL( I, MONE, B(1,I), 1 ) 729 160 CONTINUE 730 END IF 731 OPTWRK = MAX( OPTWRK, INT( DWORK(N+1) ) + N ) 732C 733C Solve the reduced generalized Lyapunov equation. 734C 735C ( Workspace: 6*N-6 ) 736C 737 IF ( ISDISC ) THEN 738 CALL SG03BU( TRANS, N, A, LDA, E, LDE, B, LDB, SCALE, DWORK, 739 $ INFO1 ) 740 IF ( INFO1 .NE. 0 ) THEN 741 IF ( INFO1 .EQ. 1 ) INFO = 1 742 IF ( INFO1 .EQ. 2 ) INFO = 3 743 IF ( INFO1 .EQ. 3 ) INFO = 6 744 IF ( INFO1 .EQ. 4 ) INFO = 7 745 IF ( INFO .NE. 1 ) 746 $ RETURN 747 END IF 748 ELSE 749 CALL SG03BV( TRANS, N, A, LDA, E, LDE, B, LDB, SCALE, DWORK, 750 $ INFO1 ) 751 IF ( INFO1 .NE. 0 ) THEN 752 IF ( INFO1 .EQ. 1 ) INFO = 1 753 IF ( INFO1 .GE. 2 ) INFO = 3 754 IF ( INFO1 .EQ. 3 ) INFO = 5 755 IF ( INFO .NE. 1 ) 756 $ RETURN 757 END IF 758 END IF 759C 760C Transform the solution matrix back. 761C 762C U := U * Q**T or U := Z * U 763C 764C Use BLAS 3 if there is enough workspace. Otherwise, use BLAS 2. 765C 766C ( Workspace: max(1,N) ) 767C 768 IF ( .NOT.ISTRAN ) THEN 769 IF ( LDWORK .GE. N*N ) THEN 770 CALL DLACPY( 'All', N, N, Q, LDQ, DWORK, N ) 771 CALL DTRMM( 'Right', 'Upper', 'Transpose', 'NonUnit', N, N, 772 $ ONE, B, LDB, DWORK, N) 773 DO 170 I = 1, N 774 CALL DCOPY( N, DWORK(N*(I-1)+1), 1, B(I,1), LDB ) 775 170 CONTINUE 776 ELSE 777 DO 180 I = 1, N 778 CALL DCOPY( N-I+1, B(I,I), LDB, DWORK, 1 ) 779 CALL DGEMV( 'NoTranspose', N, N-I+1, ONE, Q(1,I), LDQ, 780 $ DWORK, 1, ZERO, B(I,1), LDB ) 781 180 CONTINUE 782 END IF 783 ELSE 784 IF ( LDWORK .GE. N*N ) THEN 785 CALL DLACPY( 'All', N, N, Z, LDZ, DWORK, N ) 786 CALL DTRMM( 'Right', 'Upper', 'NoTranspose', 'NonUnit', N, 787 $ N, ONE, B, LDB, DWORK, N ) 788 CALL DLACPY( 'All', N, N, DWORK, N, B, LDB ) 789 ELSE 790 DO 200 I = 1, N 791 CALL DCOPY( I, B(1,I), 1, DWORK, 1 ) 792 CALL DGEMV( 'NoTranspose', N, I, ONE, Z, LDZ, DWORK, 1, 793 $ ZERO, B(1,I), 1 ) 794 200 CONTINUE 795 END IF 796 END IF 797 OPTWRK = MAX( OPTWRK, N*N ) 798C 799C Overwrite U with the triangular matrix of its QR-factorization 800C or its RQ-factorization. 801C (The entries on the main diagonal are non-negative.) 802C 803C ( Workspace: >= max(1,2*N) ) 804C 805 IF ( .NOT.ISTRAN ) THEN 806 CALL DGEQRF( N, N, B, LDB, DWORK, DWORK(N+1), LDWORK-N, INFO1 ) 807 IF ( N .GT. 1 ) 808 $ CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2,1), LDB ) 809 DO 220 I = 1, N 810 IF ( B(I,I) .LT. ZERO ) 811 $ CALL DSCAL( N+1-I, MONE, B(I,I), LDB ) 812 220 CONTINUE 813 ELSE 814 CALL DGERQF( N, N, B, LDB, DWORK, DWORK(N+1), LDWORK-N, INFO1 ) 815 IF ( N .GT. 1 ) 816 $ CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2,1), LDB ) 817 DO 240 I = 1, N 818 IF ( B(I,I) .LT. ZERO ) 819 $ CALL DSCAL( I, MONE, B(1,I), 1 ) 820 240 CONTINUE 821 END IF 822 OPTWRK = MAX( OPTWRK, INT( DWORK(N+1) ) + N ) 823C 824 DWORK(1) = DBLE( MAX( OPTWRK, MINWRK ) ) 825 RETURN 826C *** Last line of SG03BD *** 827 END 828