SUBROUTINE SB03OR( DISCR, LTRANS, N, M, S, LDS, A, LDA, C, LDC,
$ SCALE, INFO )
C
C PURPOSE
C
C To compute the solution of the Sylvester equations
C
C op(S)'*X + X*op(A) = scale*C, if DISCR = .FALSE. or
C
C op(S)'*X*op(A) - X = scale*C, if DISCR = .TRUE.
C
C where op(K) = K or K' (i.e., the transpose of the matrix K), S is
C an N-by-N block upper triangular matrix with one-by-one and
C two-by-two blocks on the diagonal, A is an M-by-M matrix (M = 1 or
C M = 2), X and C are each N-by-M matrices, and scale is an output
C scale factor, set less than or equal to 1 to avoid overflow in X.
C The solution X is overwritten on C.
C
C SB03OR is a service routine for the Lyapunov solver SB03OT.
C
C ARGUMENTS
C
C Mode Parameters
C
C DISCR LOGICAL
C Specifies the equation to be solved:
C = .FALSE.: op(S)'*X + X*op(A) = scale*C;
C = .TRUE. : op(S)'*X*op(A) - X = scale*C.
C
C LTRANS LOGICAL
C Specifies the form of op(K) to be used, as follows:
C = .FALSE.: op(K) = K (No transpose);
C = .TRUE. : op(K) = K**T (Transpose).
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, dimension (LDS,N)
C The leading N-by-N upper Hessenberg part of the array S
C must contain the block upper triangular matrix. The
C elements below the upper Hessenberg part of the array S
C are not referenced. The array S must not contain
C diagonal blocks larger than two-by-two and the two-by-two
C blocks must only correspond to complex conjugate pairs of
C 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 REFERENCES
C
C [1] Hammarling, S.J.
C Numerical solution of the stable, non-negative definite
C Lyapunov equation.
C IMA J. Num. Anal., 2, pp. 303-325, 1982.
C
C NUMERICAL ASPECTS
C 2
C The algorithm requires 0(N M) operations and is backward stable.
C
C CONTRIBUTOR
C
C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, May 1997.
C Supersedes Release 2.0 routines SB03CW and SB03CX by
C Sven Hammarling, NAG Ltd, United Kingdom, Oct. 1986.
C Partly based on routine PLYAP4 by A. Varga, University of Bochum,
C May 1992.
C
C REVISIONS
C
C December 1997, April 1998, May 1999, April 2000.
C
C KEYWORDS
C
C Lyapunov equation, orthogonal transformation, 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, LTRANS
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 Test the input scalar arguments.
C
IF( N.LT.0 ) THEN
INFO = -3
ELSE IF( .NOT.( M.EQ.1 .OR. M.EQ.2 ) ) THEN
INFO = -4
ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
INFO = -6
ELSE IF( LDA.LT.M ) THEN
INFO = -8
ELSE IF( LDC.LT.MAX( 1, N ) ) THEN
INFO = -10
END IF
C
IF ( INFO.NE.0 ) THEN
C
C Error return.
C
CALL XERBLA( 'SB03OR', -INFO )
RETURN
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 ( LTRANS ) THEN
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, L-1 ).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 S*X*A' - X = scale*C.
C
C The L-th block of X is determined from
C
C S(L,L)*X(L)*A' - X(L) = C(L) - R(L),
C
C where
C
C N
C R(L) = SUM [S(L,J)*X(J)] * A' .
C J=L+1
C
G11 = -DDOT( N-L2, S( L1, L2P1 ), LDS, C( L2P1, 1 ), 1 )
IF ( TBYT ) THEN
G12 = -DDOT( N-L2, S( L1, L2P1 ), LDS, C( L2P1, 2 ),
$ 1 )
VEC( 1, 1 ) = C( L1, 1 ) + G11*AT(1,1) + G12*AT(2,1)
VEC( 1, 2 ) = C( L1, 2 ) + G11*AT(1,2) + G12*AT(2,2)
ELSE
VEC (1, 1 ) = C( L1, 1 ) + G11*AT(1,1)
END IF
IF ( DL.NE.1 ) THEN
G21 = -DDOT( N-L2, S( L2, L2P1 ), LDS, C( L2P1, 1 ),
$ 1 )
IF ( TBYT ) THEN
G22 = -DDOT( N-L2, S( L2, L2P1 ), LDS,
$ C( L2P1, 2 ), 1 )
VEC( 2, 1 ) = C( L2, 1 ) + G21*AT(1,1) +
$ G22*AT(2,1)
VEC( 2, 2 ) = C( L2, 2 ) + G21*AT(1,2) +
$ G22*AT(2,2)
ELSE
VEC( 2, 1 ) = C( L2, 1 ) + G21*AT(1,1)
END IF
END IF
CALL SB04PX( .FALSE., .FALSE., -ISGN, DL, M, S( L1, L1 ),
$ LDS, AT, 2, VEC, 2, SCALOC, X, 2, XNORM,
$ INFO )
ELSE
C
C Solve S*X + X*A' = scale*C.
C
C The L-th block of X is determined from
C
C S(L,L)*X(L) + X(L)*A' = C(L) - R(L),
C
C where
C N
C R(L) = SUM S(L,J)*X(J) .
C J=L+1
C
VEC( 1, 1 ) = C( L1, 1 ) -
$ DDOT( N-L2, S( L1, L2P1 ), LDS,
$ C( L2P1, 1 ), 1 )
IF ( TBYT )
$ VEC( 1, 2 ) = C( L1, 2 ) -
$ DDOT( N-L2, S( L1, L2P1 ), LDS,
$ C( L2P1, 2 ), 1 )
C
IF ( DL.NE.1 ) THEN
VEC( 2, 1 ) = C( L2, 1 ) -
$ DDOT( N-L2, S( L2, L2P1 ), LDS,
$ C( L2P1, 1 ), 1 )
IF ( TBYT )
$ VEC( 2, 2 ) = C( L2, 2 ) -
$ DDOT( N-L2, S( L2, L2P1 ), LDS,
$ C( L2P1, 2 ), 1 )
END IF
CALL DLASY2( .FALSE., .FALSE., ISGN, DL, M, S( L1, L1 ),
$ LDS, AT, 2, VEC, 2, SCALOC, X, 2, XNORM,
$ INFO )
END IF
INFOM = MAX( INFO, INFOM )
IF ( SCALOC.NE.ONE ) THEN
C
DO 10 J = 1, M
CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
10 CONTINUE
C
SCALE = SCALE*SCALOC
END IF
C( L1, 1 ) = X( 1, 1 )
IF ( TBYT ) C( L1, 2 ) = X( 1, 2 )
IF ( DL.NE.1 ) THEN
C( L2, 1 ) = X( 2, 1 )
IF ( TBYT ) C( L2, 2 ) = X( 2, 2 )
END IF
20 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 = 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
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
END IF
C
INFO = INFOM
RETURN
C *** Last line of SB03OR ***
END