SB03OY

Solving (for Cholesky factor) stable 2-by-2 continuous- or discrete-time Lyapunov equations, with matrix A having complex conjugate eigenvalues

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