**Purpose**

To solve for the Cholesky factor U of X, op(U)'*op(U) = X, where U is a two-by-two upper triangular matrix, either the continuous-time two-by-two Lyapunov equation 2 op(S)'*X + X*op(S) = -ISGN*scale *op(R)'*op(R), when DISCR = .FALSE., or the discrete-time two-by-two Lyapunov equation 2 op(S)'*X*op(S) - X = -ISGN*scale *op(R)'*op(R), when DISCR = .TRUE., where op(K) = K or K' (i.e., the transpose of the matrix K), S is a two-by-two matrix with complex conjugate eigenvalues, R is a two-by-two upper triangular matrix, ISGN = -1 or 1, and scale is an output scale factor, set less than or equal to 1 to avoid overflow in X. The routine also computes two matrices, B and A, so that 2 B*U = U*S and A*U = scale *R, if LTRANS = .FALSE., or 2 U*B = S*U and U*A = scale *R, if LTRANS = .TRUE., which are used by the general Lyapunov solver. In the continuous-time case ISGN*S must be stable, so that its eigenvalues must have strictly negative real parts. In the discrete-time case S must be convergent if ISGN = 1, that is, its eigenvalues must have moduli less than unity, or S must be completely divergent if ISGN = -1, that is, its eigenvalues must have moduli greater than unity.

SUBROUTINE SB03OY( DISCR, LTRANS, ISGN, S, LDS, R, LDR, A, LDA, $ SCALE, INFO ) C .. Scalar Arguments .. LOGICAL DISCR, LTRANS INTEGER INFO, ISGN, LDA, LDR, LDS DOUBLE PRECISION SCALE C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), R(LDR,*), S(LDS,*)

**Mode Parameters**

DISCR LOGICAL Specifies the equation to be solved: 2 = .FALSE.: op(S)'*X + X*op(S) = -ISGN*scale *op(R)'*op(R); 2 = .TRUE. : op(S)'*X*op(S) - X = -ISGN*scale *op(R)'*op(R). LTRANS LOGICAL Specifies the form of op(K) to be used, as follows: = .FALSE.: op(K) = K (No transpose); = .TRUE. : op(K) = K**T (Transpose). ISGN INTEGER Specifies the sign of the equation as described before. ISGN may only be 1 or -1.

S (input/output) DOUBLE PRECISION array, dimension (LDS,2) On entry, S must contain a 2-by-2 matrix. On exit, S contains a 2-by-2 matrix B such that B*U = U*S, if LTRANS = .FALSE., or U*B = S*U, if LTRANS = .TRUE.. Notice that if U is nonsingular then B = U*S*inv( U ), if LTRANS = .FALSE. B = inv( U )*S*U, if LTRANS = .TRUE.. LDS INTEGER The leading dimension of array S. LDS >= 2. R (input/output) DOUBLE PRECISION array, dimension (LDR,2) On entry, R must contain a 2-by-2 upper triangular matrix. The element R( 2, 1 ) is not referenced. On exit, R contains U, the 2-by-2 upper triangular Cholesky factor of the solution X, X = op(U)'*op(U). LDR INTEGER The leading dimension of array R. LDR >= 2. A (output) DOUBLE PRECISION array, dimension (LDA,2) A contains a 2-by-2 upper triangular matrix A satisfying A*U/scale = scale*R, if LTRANS = .FALSE., or U*A/scale = scale*R, if LTRANS = .TRUE.. Notice that if U is nonsingular then A = scale*scale*R*inv( U ), if LTRANS = .FALSE. A = scale*scale*inv( U )*R, if LTRANS = .TRUE.. LDA INTEGER The leading dimension of array A. LDA >= 2. SCALE (output) DOUBLE PRECISION The scale factor, scale, set less than or equal to 1 to prevent the solution overflowing.

INFO INTEGER = 0: successful exit; = 1: if the Lyapunov equation is (nearly) singular (warning indicator); if DISCR = .FALSE., this means that while the matrix S has computed eigenvalues with negative real parts, it is only just stable in the sense that small perturbations in S can make one or more of the eigenvalues have a non-negative real part; if DISCR = .TRUE., this means that while the matrix S has computed eigenvalues inside the unit circle, it is nevertheless only just convergent, in the sense that small perturbations in S can make one or more of the eigenvalues lie outside the unit circle; perturbed values were used to solve the equation (but the matrix S is unchanged); = 2: if DISCR = .FALSE., and ISGN*S is not stable or if DISCR = .TRUE., ISGN = 1 and S is not convergent or if DISCR = .TRUE., ISGN = -1 and S is not completely divergent; = 4: if S has real eigenvalues. NOTE: In the interests of speed, this routine does not check all inputs for errors.

The LAPACK scheme for solving 2-by-2 Sylvester equations is adapted for 2-by-2 Lyapunov equations, but directly computing the Cholesky factor of the solution.

[1] Hammarling S. J. Numerical solution of the stable, non-negative definite Lyapunov equation. IMA J. Num. Anal., 2, pp. 303-325, 1982.

The algorithm is backward stable.

None

**Program Text**

None

None

None