SUBROUTINE LRAUX1( DISCR, UPPER, N, M, S, LDS, A, LDA, C, LDC, \$ SCALE, INFO ) C C PURPOSE C C To compute the solution of the Sylvester equation C C S'*X + X*A = scale*C, if DISCR = .FALSE., C or S'*X*A - X = scale*C, if DISCR = .TRUE., C C where S is an N-by-N block upper or lower triangular matrix with C one-by-one and two-by-two blocks on the diagonal, A is an M-by-M C matrix (M = 1 or M = 2), X and C are each N-by-M matrices, and C scale is an output scale factor, set less than or equal to 1 to C avoid overflow in X. C The solution X is overwritten on C. C C LRAUX1 is a service routine for the Lyapunov solver LRLYAP. C It is a modified version of the SLICOT routine SB03OR. C C ARGUMENTS C C Mode Parameters C C DISCR LOGICAL C Specifies the equation to be solved: C = .FALSE.: S'*X + X*A = scale*C; C = .TRUE. : S'*X*A - X = scale*C. C C UPPER LOGICAL C Specifies the shape of S as follows: C = .TRUE. : S is block upper triangular; C = .FALSE.: S is block lower triangular. C C Input/Output Parameters C C N (input) INTEGER C The order of the matrix S and also the number of rows of C matrices X and C. N >= 0. C C M (input) INTEGER C The order of the matrix A and also the number of columns C of matrices X and C. M = 1 or M = 2. C C S (input) DOUBLE PRECISION array of dimension (LDS,N) C If UPPER = .TRUE., the leading N-by-N upper Hessenberg C part of this array must contain the block upper C triangular matrix S. The elements below the upper C Hessenberg part of the array S are not referenced. C If UPPER = .FALSE., the leading N-by-N lower Hessenberg C part of this array must contain the block lower C triangular matrix S. The elements above the lower C Hessenberg part of the array S are not referenced. C The 2-by-2 blocks of S must only correspond to complex C conjugate pairs of eigenvalues (not to real eigenvalues). C C LDS INTEGER C The leading dimension of array S. LDS >= MAX(1,N). C C A (input) DOUBLE PRECISION array, dimension (LDS,M) C The leading M-by-M part of this array must contain a C given matrix, where M = 1 or M = 2. C C LDA INTEGER C The leading dimension of array A. LDA >= M. C C C (input/output) DOUBLE PRECISION array, dimension (LDC,M) C On entry, C must contain an N-by-M matrix, where M = 1 or C M = 2. C On exit, C contains the N-by-M matrix X, the solution of C the Sylvester equation. C C LDC INTEGER C The leading dimension of array C. LDC >= MAX(1,N). C C SCALE (output) DOUBLE PRECISION C The scale factor, scale, set less than or equal to 1 to C prevent the solution overflowing. C C Error Indicator C C INFO INTEGER C = 0: successful exit; C = 1: if DISCR = .FALSE., and S and -A have common C eigenvalues, or if DISCR = .TRUE., and S and A have C eigenvalues whose product is equal to unity; C a solution has been computed using slightly C perturbed values. C C METHOD C C The LAPACK scheme for solving Sylvester equations is adapted. C C C NUMERICAL ASPECTS C 2 C The algorithm requires 0(N ) operations. C C CONTRIBUTOR C C The SLICOT routine SB03OR, upon which this routine is based, C was originally written by Sven Hammarling, Vasile Sima, and C Andras Varga. C C Daniel Kressner, Umea University, June 2006. C C KEYWORDS C C Sylvester equation, real Schur form. C C ****************************************************************** C C .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) C .. Scalar Arguments .. LOGICAL DISCR, UPPER INTEGER INFO, LDA, LDS, LDC, M, N DOUBLE PRECISION SCALE C .. C .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), S( LDS, * ) C .. Local Scalars .. LOGICAL TBYT INTEGER DL, INFOM, ISGN, J, L, L1, L2, L2P1, LNEXT DOUBLE PRECISION G11, G12, G21, G22, SCALOC, XNORM C .. C .. Local Arrays .. DOUBLE PRECISION AT( 2, 2 ), VEC( 2, 2 ), X( 2, 2 ) C .. C .. External Functions .. DOUBLE PRECISION DDOT EXTERNAL DDOT C .. C .. External Subroutines .. EXTERNAL DLASY2, DSCAL, SB04PX, XERBLA C .. Intrinsic Functions .. INTRINSIC MAX, MIN C .. C .. Executable Statements .. C INFO = 0 C C No test of input scalar arguments for efficiency. Uncomment these C lines to include testing. C C IF( N.LT.0 ) THEN C INFO = -3 C ELSE IF( .NOT.( M.EQ.1 .OR. M.EQ.2 ) ) THEN C INFO = -4 C ELSE IF( LDS.LT.MAX( 1, N ) ) THEN C INFO = -6 C ELSE IF( LDA.LT.M ) THEN C INFO = -8 C ELSE IF( LDC.LT.MAX( 1, N ) ) THEN C INFO = -10 C END IF C C IF ( INFO.NE.0 ) THEN C C Error return. C C CALL XERBLA( 'LRAUX1', -INFO ) C RETURN C END IF C SCALE = ONE C C Quick return if possible. C IF ( N.EQ.0 ) \$ RETURN C ISGN = 1 TBYT = M.EQ.2 INFOM = 0 C C Construct A'. C AT(1,1) = A(1,1) IF ( TBYT ) THEN AT(1,2) = A(2,1) AT(2,1) = A(1,2) AT(2,2) = A(2,2) END IF C IF ( UPPER ) THEN C C Start row loop (index = L). C L1 (L2) : row index of the first (last) row of X(L). C LNEXT = 1 C DO 40 L = 1, N IF( L.LT.LNEXT ) \$ GO TO 40 L1 = L L2 = L IF( L.LT.N ) THEN IF( S( L+1, L ).NE.ZERO ) \$ L2 = L2 + 1 LNEXT = L2 + 1 END IF L2P1 = MIN( L2+1, N ) DL = L2 - L1 + 1 C IF ( DISCR ) THEN C C Solve A'*X'*S - X' = scale*C'. C C The L-th block of X is determined from C C A'*X(L)'*S(L,L) - X(L)' = C(L)' - R(L), C C where C C L-1 C R(L) = A' * SUM [X(J)'*S(J,L)] . C J=1 C G11 = -DDOT( L1-1, C, 1, S( 1, L1 ), 1 ) IF ( TBYT ) THEN G21 = -DDOT( L1-1, C( 1, 2 ), 1, S( 1, L1 ), 1 ) VEC( 1, 1 ) = C( L1, 1 ) + AT(1,1)*G11 + AT(1,2)*G21 VEC( 2, 1 ) = C( L1, 2 ) + AT(2,1)*G11 + AT(2,2)*G21 ELSE VEC (1, 1 ) = C( L1, 1 ) + AT(1,1)*G11 END IF IF ( DL .NE. 1 ) THEN G12 = -DDOT( L1-1, C, 1, S( 1, L2 ), 1 ) IF ( TBYT ) THEN G22 = -DDOT( L1-1, C( 1, 2 ), 1, S( 1, L2 ), 1 ) VEC( 1, 2 ) = C( L2, 1 ) + AT(1,1)*G12 + \$ AT(1,2)*G22 VEC( 2, 2 ) = C( L2, 2 ) + AT(2,1)*G12 + \$ AT(2,2)*G22 ELSE VEC( 1, 2 ) = C( L2, 1 ) + AT(1,1)*G12 END IF END IF CALL SB04PX( .FALSE., .FALSE., -ISGN, M, DL, AT, 2, \$ S( L1, L1 ), LDS, VEC, 2, SCALOC, X, 2, \$ XNORM, INFO ) ELSE C C Solve A'*X' + X'*S = scale*C'. C C The L-th block of X is determined from C C A'*X(L)' + X(L)'*S(L,L) = C(L)' - R(L), C C where C L-1 C R(L) = SUM [X(J)'*S(J,L)]. C J=1 C VEC( 1, 1 ) = C( L1, 1 ) - \$ DDOT( L1-1, C, 1, S( 1, L1 ), 1 ) IF ( TBYT ) \$ VEC( 2, 1 ) = C( L1, 2 ) - \$ DDOT( L1-1, C( 1, 2 ), 1, S( 1, L1 ), 1) C IF ( DL.NE.1 ) THEN VEC( 1, 2 ) = C( L2, 1 ) - \$ DDOT( L1-1, C, 1, S( 1, L2 ), 1 ) IF ( TBYT ) \$ VEC( 2, 2 ) = C( L2, 2 ) - \$ DDOT( L1-1, C( 1, 2 ), 1, S( 1, L2 ), 1) END IF CALL DLASY2( .FALSE., .FALSE., ISGN, M, DL, AT, 2, \$ S( L1, L1 ), LDS, VEC, 2, SCALOC, X, 2, \$ XNORM, INFO ) END IF INFOM = MAX( INFO, INFOM ) IF ( SCALOC.NE.ONE ) THEN C DO 30 J = 1, M CALL DSCAL( N, SCALOC, C( 1, J ), 1 ) 30 CONTINUE C SCALE = SCALE*SCALOC END IF C( L1, 1 ) = X( 1, 1 ) IF ( TBYT ) C( L1, 2 ) = X( 2, 1 ) IF ( DL.NE.1 ) THEN C( L2, 1 ) = X( 1, 2 ) IF ( TBYT ) C( L2, 2 ) = X( 2, 2 ) END IF 40 CONTINUE C ELSE C C Start row loop (index = L). C L1 (L2) : row index of the first (last) row of X(L). C LNEXT = N C DO 20 L = N, 1, -1 IF( L.GT.LNEXT ) \$ GO TO 20 L1 = L L2 = L IF( L.GT.1 ) THEN IF( S( L-1, L ).NE.ZERO ) \$ L1 = L1 - 1 LNEXT = L1 - 1 END IF DL = L2 - L1 + 1 L2P1 = MIN( L2+1, N ) C IF ( DISCR ) THEN C C Solve A'*X'*S - X' = scale*C'. C C The L-th block of X is determined from C C A'*X(L)'*S(L,L) - X(L)' = C(L)' - R(L), C C where C C N C R(L) = A' * SUM [X(J)'*S(J,L)] . C J=L+1 C G11 = -DDOT( N-L2, C(L2P1,1), 1, S(L2P1,L1), 1 ) IF ( TBYT ) THEN G21 = -DDOT( N-L2, C(L2P1,2), 1, S(L2P1,L1), 1 ) VEC( 1, 1 ) = C( L1, 1 ) + AT(1,1)*G11 + AT(1,2)*G21 VEC( 2, 1 ) = C( L1, 2 ) + AT(2,1)*G11 + AT(2,2)*G21 ELSE VEC (1, 1 ) = C( L1, 1 ) + AT(1,1)*G11 END IF IF ( DL .NE. 1 ) THEN G12 = -DDOT( N-L2, C(L2P1,1), 1, S(L2P1,L2), 1 ) IF ( TBYT ) THEN G22 = -DDOT( N-L2, C(L2P1,2), 1, S(L2P1,L2), 1 ) VEC(1,2) = C(L2,1) + AT(1,1)*G12 + AT(1,2)*G22 VEC(2,2) = C(L2,2) + AT(2,1)*G12 + AT(2,2)*G22 ELSE VEC(1,2) = C( L2, 1 ) + AT(1,1)*G12 END IF END IF CALL SB04PX( .FALSE., .FALSE., -ISGN, M, DL, AT, 2, \$ S( L1, L1 ), LDS, VEC, 2, SCALOC, X, 2, \$ XNORM, INFO ) ELSE C C Solve A*X' + X'*S = scale*C'. C C The L-th block of X is determined from C C A'*X(L)' + X(L)'*S(L,L) = C(L)' - R(L), C C where C N C R(L) = SUM [X(J)'*S(J,L)]. C J=L+1 C VEC( 1, 1 ) = C( L1, 1 ) - DDOT( N-L2, C( L2P1,1 ), 1, \$ S( L2P1, L1 ), 1 ) IF ( TBYT ) \$ VEC( 2, 1 ) = C( L1, 2 ) - DDOT( N-L2, C( L2P1, 2 ), \$ 1, S( L2P1, L1 ), 1 ) C IF ( DL.NE.1 ) THEN VEC( 1, 2 ) = C( L2, 1 ) - DDOT( N-L2, C(L2P1,1), 1, \$ S( L2P1, L2 ), 1 ) IF ( TBYT ) \$ VEC( 2, 2 ) = C( L2, 2 ) - \$ DDOT( N-L2, C( L2P1, 2 ), 1, \$ S( L2P1, L2 ), 1) END IF CALL DLASY2( .FALSE., .FALSE., ISGN, M, DL, AT, 2, \$ S( L1, L1 ), LDS, VEC, 2, SCALOC, X, 2, \$ XNORM, INFO ) END IF INFOM = MAX( INFO, INFOM ) IF ( SCALOC.NE.ONE ) THEN C DO 130 J = 1, M CALL DSCAL( N, SCALOC, C( 1, J ), 1 ) 130 CONTINUE C SCALE = SCALE*SCALOC END IF C( L1, 1 ) = X( 1, 1 ) C IF ( TBYT ) C( L1, 2 ) = X( 2, 1 ) IF ( DL.NE.1 ) THEN C( L2, 1 ) = X( 1, 2 ) IF ( TBYT ) C( L2, 2 ) = X( 2, 2 ) END IF C 20 CONTINUE END IF C INFO = INFOM RETURN C *** Last line of LRAUX1 *** END