path: root/modules/cacsd/src/slicot/sb03ot.f

      SUBROUTINE SB03OT( DISCR, LTRANS, N, S, LDS, R, LDR, SCALE, DWORK,
     $                   INFO )
C
C     RELEASE 4.0, WGS COPYRIGHT 1999.
C
C     PURPOSE
C
C     To solve for X = op(U)'*op(U) either the stable non-negative
C     definite continuous-time Lyapunov equation
C                                   2
C        op(S)'*X + X*op(S) = -scale *op(R)'*op(R)                   (1)
C
C     or the convergent non-negative definite discrete-time Lyapunov
C     equation
C                                   2
C        op(S)'*X*op(S) - X = -scale *op(R)'*op(R)                   (2)
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 or
C     two-by-two blocks on the diagonal, R is an N-by-N upper triangular
C     matrix, and scale is an output scale factor, set less than or
C     equal to 1 to avoid overflow in X.
C
C     In the case of equation (1) the matrix S must be stable (that
C     is, all the eigenvalues of S must have negative real parts),
C     and for equation (2) the matrix S must be convergent (that is,
C     all the eigenvalues of S must lie inside the unit circle).
C
C     ARGUMENTS 
C
C     Mode Parameters
C
C     DISCR   LOGICAL
C             Specifies the type of Lyapunov equation to be solved as
C             follows:
C             = .TRUE. :  Equation (2), discrete-time case;
C             = .FALSE.:  Equation (1), continuous-time case.
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 matrices S and R.  N >= 0.
C
C     S       (input) DOUBLE PRECISION array of dimension (LDS,N)
C             The leading N-by-N upper Hessenberg part of this array
C             must contain the block upper triangular matrix.
C             The elements below the upper Hessenberg part of the array
C             S are not referenced. The 2-by-2 blocks must only
C             correspond to complex conjugate pairs of eigenvalues (not
C             to real eigenvalues).
C
C     LDS     INTEGER
C             The leading dimension of array S.  LDS >= MAX(1,N).
C
C     R       (input/output) DOUBLE PRECISION array of dimension (LDR,N)
C             On entry, the leading N-by-N upper triangular part of this
C             array must contain the upper triangular matrix R.
C             On exit, the leading N-by-N upper triangular part of this
C             array contains the upper triangular matrix U.
C             The strict lower triangle of R is not referenced.
C
C     LDR     INTEGER
C             The leading dimension of array R.  LDR >= 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     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (4*N)
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal 
C                   value;
C             = 1:  if the Lyapunov equation is (nearly) singular
C                   (warning indicator);
C                   if DISCR = .FALSE., this means that while the
C                   matrix S has computed eigenvalues with negative real
C                   parts, it is only just stable in the sense that
C                   small perturbations in S can make one or more of the
C                   eigenvalues have a non-negative real part;
C                   if DISCR = .TRUE., this means that while the
C                   matrix S has computed eigenvalues inside the unit
C                   circle, it is nevertheless only just convergent, in
C                   the sense that small perturbations in S can make one
C                   or more of the eigenvalues lie outside the unit
C                   circle;
C                   perturbed values were used to solve the equation 
C                   (but the matrix S is unchanged);
C             = 2:  if the matrix S is not stable (that is, one or more
C                   of the eigenvalues of S has a non-negative real
C                   part), if DISCR = .FALSE., or not convergent (that
C                   is, one or more of the eigenvalues of S lies outside
C                   the unit circle), if DISCR = .TRUE.;
C             = 3:  if the matrix S has two or more consecutive non-zero
C                   elements on the first sub-diagonal, so that there is
C                   a block larger than 2-by-2 on the diagonal;
C             = 4:  if the matrix S has a 2-by-2 diagonal block with
C                   real eigenvalues instead of a complex conjugate
C                   pair.
C
C     METHOD
C
C     The method used by the routine is based on a variant of the
C     Bartels and Stewart backward substitution method [1], that finds
C     the Cholesky factor op(U) directly without first finding X and
C     without the need to form the normal matrix op(R)'*op(R) [2].
C
C     The continuous-time Lyapunov equation in the canonical form
C                                                        2
C       op(S)'*op(U)'*op(U) + op(U)'*op(U)*op(S) = -scale *op(R)'*op(R),
C
C     or the discrete-time Lyapunov equation in the canonical form
C                                                        2
C       op(S)'*op(U)'*op(U)*op(S) - op(U)'*op(U) = -scale *op(R)'*op(R),
C
C     where U and R are upper triangular, is solved for U.
C
C     REFERENCES
C
C     [1] Bartels, R.H. and Stewart, G.W.
C         Solution of the matrix equation  A'X + XB = C.
C         Comm. A.C.M., 15, pp. 820-826, 1972.
C
C     [2] 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                               3
C     The algorithm requires 0(N ) operations and is backward stable.
C
C     FURTHER COMMENTS
C
C     The Lyapunov equation may be very ill-conditioned. In particular
C     if S is only just stable (or convergent) then the Lyapunov
C     equation will be ill-conditioned. "Large" elements in U relative
C     to those of S and R, or a "small" value for scale, is a symptom
C     of ill-conditioning. A condition estimate can be computed using
C     SLICOT Library routine SB03MD.
C
C     CONTRIBUTOR
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, May 1997.
C     Supersedes Release 2.0 routine SB03CZ by Sven Hammarling,
C     NAG Ltd, United Kingdom, Oct. 1986.
C     Partly based on SB03CZ and PLYAP1 by A. Varga, University of
C     Bochum, May 1992.
C
C     REVISIONS
C
C     Dec. 1997, April 1998, May 1999.
C
C     KEYWORDS
C
C     Lyapunov equation, orthogonal transformation, real Schur form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE, TWO
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
C     .. Scalar Arguments ..
      LOGICAL           DISCR, LTRANS
      INTEGER           INFO, LDR, LDS, N
      DOUBLE PRECISION  SCALE
C     .. Array Arguments ..
      DOUBLE PRECISION  DWORK(*), R(LDR,*), S(LDS,*)
C     .. Local Scalars ..
      LOGICAL           CONT, TBYT
      INTEGER           INFOM, ISGN, J, J1, J2, J3, K, K1, K2, K3, 
     $                  KOUNT, KSIZE
      DOUBLE PRECISION  ABSSKK, ALPHA, BIGNUM, D1, D2, DR, EPS, SCALOC,
     $                  SMIN, SMLNUM, SUM, T1, T2, T3, T4, TAU1, TAU2,
     $                  TEMP, V1, V2, V3, V4
C     .. Local Arrays ..
      DOUBLE PRECISION  A(2,2), B(2,2), U(2,2)
C     .. External Functions ..
      DOUBLE PRECISION  DLAMCH, DLANHS
      EXTERNAL          DLAMCH, DLANHS
C     .. External Subroutines ..
      EXTERNAL          DAXPY, DCOPY, DLABAD, DLARFG, DSCAL, DSWAP,
     $                  DTRMM, DTRMV, MB04ND, MB04OD, SB03OR, SB03OY,
     $                  XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, SIGN, SQRT
C     .. Executable Statements ..
C
      INFO = 0
C
C     Test the input scalar arguments.
C
      IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDR.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
C
      IF ( INFO.NE.0 ) THEN      
C
C        Error return.
C
         CALL XERBLA( 'SB03OT', -INFO )
         RETURN
      END IF
C
      SCALE = ONE
C
C     Quick return if possible.
C
      IF (N.EQ.0)
     $   RETURN
C
C     Set constants to control overflow.
C
      EPS = DLAMCH( 'P' )
      SMLNUM = DLAMCH( 'S' )
      BIGNUM = ONE / SMLNUM
      CALL DLABAD( SMLNUM, BIGNUM )
      SMLNUM = SMLNUM*DBLE( N*N ) / EPS
      BIGNUM = ONE / SMLNUM
C
      SMIN = MAX( SMLNUM, EPS*DLANHS( 'Max', N, S, LDS, DWORK ) )
      INFOM = 0
C
C     Start the solution. Most of the comments refer to notation and
C     equations in sections 5 and 10 of the second reference above.
C
C     Determine whether or not the current block is two-by-two.
C     K gives the position of the start of the current block and
C     TBYT is true if the block is two-by-two.
C
      CONT = .NOT.DISCR
      ISGN = 1
      IF ( .NOT.LTRANS ) THEN
C
C        Case op(M) = M.
C
         KOUNT = 1
C
   10    CONTINUE
C        WHILE( KOUNT.LE.N )LOOP
         IF ( KOUNT.LE.N ) THEN
            K = KOUNT
            IF ( KOUNT.GE.N ) THEN
               TBYT  = .FALSE.
               KOUNT = KOUNT + 1
            ELSE IF ( S(K+1,K).EQ.ZERO ) THEN
               TBYT  = .FALSE.
               KOUNT = KOUNT + 1
            ELSE
               TBYT = .TRUE.
               IF ( (K+1).LT.N ) THEN
                  IF ( S(K+2,K+1).NE.ZERO ) THEN
                     INFO = 3
                     RETURN
                  END IF
               END IF
               KOUNT = KOUNT + 2
            END IF
            IF ( TBYT ) THEN
C
C              Solve the two-by-two Lyapunov equation (6.1) or (10.19),
C              using the routine SB03OY.
C
               B(1,1) = S(K,K)
               B(2,1) = S(K+1,K)
               B(1,2) = S(K,K+1)
               B(2,2) = S(K+1,K+1)
               U(1,1) = R(K,K)
               U(1,2) = R(K,K+1)
               U(2,2) = R(K+1,K+1)
C
               CALL SB03OY( DISCR, LTRANS, ISGN, B, 2, U, 2, A, 2,
     $                      SCALOC, INFO )
               IF ( INFO.GT.1 )
     $            RETURN
               INFOM = MAX( INFO, INFOM )
               IF( SCALOC.NE.ONE ) THEN
C
                  DO 20 J = 1, N
                     CALL DSCAL( J, SCALOC, R(1,J), 1 )
   20             CONTINUE
C
                  SCALE = SCALE*SCALOC
               END IF
               R(K,K)     = U(1,1)
               R(K,K+1)   = U(1,2)
               R(K+1,K+1) = U(2,2)
C
C              If we are not at the end of S then set up and solve
C              equation (6.2) or (10.20).
C
C              Note that  SB03OY  returns  ( u11*s11*inv( u11 ) ) in  B
C              and returns scaled alpha in  A.  ksize is the order of 
C              the remainder of  S.  k1, k2 and k3  point to the start
C              of vectors in  DWORK.
C
               IF ( KOUNT.LE.N ) THEN
                  KSIZE = N - K - 1
                  K1 = KSIZE + 1
                  K2 = KSIZE + K1
                  K3 = KSIZE + K2
C
C                 Form the right-hand side of (6.2) or (10.20), the
C                 first column in DWORK( 1 ) ,..., DWORK( n - k - 1 )
C                 the second in DWORK( n - k ) ,...,
C                 DWORK( 2*( n - k - 1 ) ).
C
                  CALL DCOPY( KSIZE, R(K,K+2), LDR, DWORK, 1 )
                  CALL DCOPY( KSIZE, R(K+1,K+2), LDR, DWORK(K1), 1 )
                  CALL DTRMM( 'Right', 'Upper', 'No transpose',
     $                        'Non-unit', KSIZE, 2, -ONE, A, 2, DWORK,
     $                        KSIZE )
                  IF ( CONT ) THEN
                     CALL DAXPY( KSIZE, -R(K,K), S(K,K+2), LDS, DWORK,
     $                           1 )
                     CALL DAXPY( KSIZE, -R(K,K+1), S(K+1,K+2), LDS,
     $                           DWORK, 1)
                     CALL DAXPY( KSIZE, -R(K+1,K+1), S(K+1,K+2), LDS,
     $                           DWORK(K1), 1 )
                  ELSE
                     CALL DAXPY( KSIZE, -R(K,K)*B(1,1), S(K,K+2), LDS,
     $                           DWORK, 1 )
                     CALL DAXPY( KSIZE, -( R(K,K+1)*B(1,1) + R(K+1,K+1)
     $                           *B(2,1) ), S(K+1,K+2), LDS, DWORK, 1 )
                     CALL DAXPY( KSIZE, -R(K,K)*B(1,2), S(K,K+2), LDS,
     $                           DWORK(K1), 1 )
                     CALL DAXPY( KSIZE, -( R(K,K+1)*B(1,2) + R(K+1,K+1)
     $                           *B(2,2) ), S(K+1,K+2), LDS, DWORK(K1),
     $                           1 )
                  END IF
C
C                 SB03OR  solves the Sylvester equations. The solution
C                 is overwritten on DWORK.
C
                  CALL SB03OR( DISCR, LTRANS, KSIZE, 2, S(K+2,K+2), LDS,
     $                         B, 2, DWORK, KSIZE, SCALOC, INFO )
                  INFOM = MAX( INFO, INFOM )
                  IF( SCALOC.NE.ONE ) THEN
C
                     DO 30 J = 1, N
                        CALL DSCAL( J, SCALOC, R(1,J), 1 )
   30                CONTINUE
C
                     SCALE = SCALE*SCALOC
                  END IF
C
C                 Copy the solution into the next  2*( n - k - 1 )
C                 elements of  DWORK.
C
                  CALL DCOPY( 2*KSIZE, DWORK, 1, DWORK(K2), 1 )
C
C                 Now form the matrix  Rhat  of equation (6.4) or
C                 (10.22). Note that (10.22) is incorrect, so here we
C                 implement a corrected version of (10.22).
C
                  IF ( CONT ) THEN
C
C                    Swap the two rows of R with DWORK.
C
                     CALL DSWAP( KSIZE, DWORK, 1, R(K,K+2), LDR )
                     CALL DSWAP( KSIZE, DWORK(K1), 1, R(K+1,K+2), LDR )
C
C                    1st column:
C
                     CALL DAXPY( KSIZE, -A(1,1), DWORK(K2), 1, DWORK,
     $                           1 )
                     CALL DAXPY( KSIZE, -A(1,2), DWORK(K3), 1, DWORK,
     $                           1 )
C
C                    2nd column:
C
                     CALL DAXPY( KSIZE, -A(2,2), DWORK(K3), 1,
     $                           DWORK(K1), 1 )
                  ELSE
C
C                    Form  v = S1'*u + s*u11', overwriting  v  on DWORK.
C
C                    Compute  S1'*u,  first multiplying by the
C                    triangular part of  S1.
C
                     CALL DTRMM( 'Left', 'Upper', 'Transpose',
     $                           'Non-unit', KSIZE, 2, ONE, S(K+2,K+2),
     $                           LDS, DWORK, KSIZE )
C
C                    Then multiply by the subdiagonal of  S1  and add in
C                    to the above result.
C
                     J1 = K1
                     J2 = K + 2
C
                     DO 40 J = 1, KSIZE-1
                        IF ( S(J2+1,J2).NE.ZERO ) THEN
                           DWORK(J)  = S(J2+1,J2)*DWORK(K2+J) + DWORK(J)
                           DWORK(J1) = S(J2+1,J2)*DWORK(K3+J) +
     $                                 DWORK(J1)
                        END IF
                        J1 = J1 + 1
                        J2 = J2 + 1
   40                CONTINUE
C
C                    Add in s*u11'.
C
                     CALL DAXPY( KSIZE, R(K,K), S(K,K+2), LDS, DWORK,
     $                           1 )
                     CALL DAXPY( KSIZE, R(K,K+1), S(K+1,K+2), LDS,
     $                           DWORK, 1 )
                     CALL DAXPY( KSIZE, R(K+1,K+1), S(K+1,K+2), LDS,
     $                           DWORK(K1), 1 )
C
C                    Next recover r from R, swapping r with u.
C
                     CALL DSWAP( KSIZE, DWORK(K2), 1, R(K,K+2), LDR )
                     CALL DSWAP( KSIZE, DWORK(K3), 1, R(K+1,K+2), LDR )
C
C                    Now we perform the QR factorization.
C
C                    ( a ) = Q*( t ),
C                    ( b )
C
C                    and form
C
C                    ( p' ) = Q'*( r' ).
C                    ( y' )      ( v' )
C
C                    y  is then the correct vector to use in (10.22).
C                    Note that  a  is upper triangular and that  t  and
C                    p  are not required.
C
                     CALL DLARFG( 3, A(1,1), B(1,1), 1, TAU1 )
                     V1  = B(1,1)
                     T1  = TAU1*V1
                     V2  = B(2,1)
                     T2  = TAU1*V2
                     SUM = A(1,2) + V1*B(1,2) + V2*B(2,2)
                     B(1,2) = B(1,2) - SUM*T1
                     B(2,2) = B(2,2) - SUM*T2
                     CALL DLARFG( 3, A(2,2), B(1,2), 1, TAU2 )
                     V3 = B(1,2)
                     T3 = TAU2*V3
                     V4 = B(2,2)
                     T4 = TAU2*V4
                     J1 = K1
                     J2 = K2
                     J3 = K3
C
                     DO 50 J = 1, KSIZE
                        SUM = DWORK(J2) + V1*DWORK(J) + V2*DWORK(J1)
                        D1  = DWORK(J)  - SUM*T1
                        D2  = DWORK(J1) - SUM*T2
                        SUM = DWORK(J3) + V3*D1 + V4*D2
                        DWORK(J)  =  D1 - SUM*T3
                        DWORK(J1) =  D2 - SUM*T4
                        J1 = J1 + 1
                        J2 = J2 + 1
                        J3 = J3 + 1
   50                CONTINUE
C
                  END IF
C
C                 Now update  R1  to give  Rhat.
C
                  CALL DCOPY( KSIZE, DWORK, 1, DWORK(K2), 1 )
                  CALL DCOPY( KSIZE, DWORK(K1), 1, DWORK(2), 2 )
                  CALL DCOPY( KSIZE, DWORK(K2), 1, DWORK(1), 2 )
                  CALL MB04OD( 'Full', KSIZE, 0, 2, R(K+2,K+2), LDR,
     $                         DWORK, 2, DWORK, 1, DWORK, 1, DWORK(K2),
     $                         DWORK(K3) )
               END IF
            ELSE
C
C              1-by-1 block.
C
C              Make sure S is stable or convergent and find u11 in
C              equation (5.13) or (10.15).
C
               IF ( DISCR ) THEN
                  ABSSKK = ABS( S(K,K) )
                  IF ( ( ABSSKK - ONE ).GE.ZERO ) THEN
                     INFO = 2
                     RETURN
                  END IF
                  TEMP = SQRT( ( ONE - ABSSKK )*( ONE + ABSSKK ) )
               ELSE
                  IF ( S(K,K).GE.ZERO ) THEN
                     INFO = 2
                     RETURN
                  END IF
                  TEMP = SQRT( ABS( TWO*S(K,K) ) )
               END IF
C
               SCALOC = ONE
               IF( TEMP.LT.SMIN ) THEN
                  TEMP  = SMIN
                  INFOM = 1
               END IF
               DR = ABS( R(K,K) )
               IF( TEMP.LT.ONE .AND. DR.GT.ONE ) THEN
                  IF( DR.GT.BIGNUM*TEMP )
     $               SCALOC = ONE / DR
               END IF
               ALPHA = SIGN( TEMP, R(K,K) )
               R(K,K) = R(K,K)/ALPHA
               IF( SCALOC.NE.ONE ) THEN
C
                  DO 60 J = 1, N
                     CALL DSCAL( J, SCALOC, R(1,J), 1 )
   60             CONTINUE
C
                  SCALE = SCALE*SCALOC
               END IF
C
C              If we are not at the end of  S  then set up and solve
C              equation (5.14) or (10.16).  ksize is the order of the
C              remainder of  S.  k1 and k2 point to the start of vectors
C              in  DWORK.
C
               IF ( KOUNT.LE.N ) THEN
                  KSIZE = N - K
                  K1 = KSIZE + 1
                  K2 = KSIZE + K1
C
C                 Form the right-hand side in DWORK( 1 ),...,
C                 DWORK( n - k ).
C
                  CALL DCOPY( KSIZE, R(K,K+1), LDR, DWORK, 1 )
                  CALL DSCAL( KSIZE, -ALPHA, DWORK, 1 )
                  IF ( CONT ) THEN
                     CALL DAXPY( KSIZE, -R(K,K), S(K,K+1), LDS, DWORK,
     $                          1 )
                  ELSE
                     CALL DAXPY( KSIZE, -S(K,K)*R(K,K), S(K,K+1), LDS,
     $                          DWORK, 1 )
                  END IF
C
C                 SB03OR solves the Sylvester equations. The solution is
C                 overwritten on  DWORK.
C
                  CALL SB03OR( DISCR, LTRANS, KSIZE, 1, S(K+1,K+1), LDS,
     $                         S(K,K), 1, DWORK, KSIZE, SCALOC, INFO )
                  INFOM = MAX( INFO, INFOM )
                  IF( SCALOC.NE.ONE ) THEN
C
                     DO 70 J = 1, N
                        CALL DSCAL( J, SCALOC, R(1,J), 1 )
   70                CONTINUE
C
                     SCALE = SCALE*SCALOC
                  END IF
C
C                 Copy the solution into the next  ( n - k ) elements
C                 of  DWORK,  copy the solution back into  R  and copy
C                 the row of  R  back into  DWORK.
C
                  CALL DCOPY( KSIZE, DWORK, 1, DWORK(K1), 1 )
                  CALL DSWAP( KSIZE, DWORK, 1, R(K,K+1), LDR )
C
C                 Now form the matrix  Rhat  of equation (5.15) or
C                 (10.17), first computing  y  in  DWORK,  and then
C                 updating  R1.
C
                  IF ( CONT ) THEN
                     CALL DAXPY( KSIZE, -ALPHA, DWORK(K1), 1, DWORK, 1 )
                  ELSE
C
C                    First form  lambda( 1 )*r  and then add in
C                    alpha*u11*s.
C
                     CALL DSCAL( KSIZE, -S(K,K), DWORK, 1 )
                     CALL DAXPY( KSIZE, ALPHA*R(K,K), S(K,K+1), LDS,
     $                           DWORK, 1 )
C
C                    Now form  alpha*S1'*u,  first multiplying by the
C                    sub-diagonal of  S1  and then the triangular part
C                    of  S1,  and add the result in DWORK.
C
                     J1 = K + 1
C
                     DO 80 J = 1, KSIZE-1
                        IF ( S(J1+1,J1).NE.ZERO ) DWORK(J)
     $                         = ALPHA*S(J1+1,J1)*DWORK(K1+J) + DWORK(J)
                        J1 = J1 + 1
   80                CONTINUE
C
                     CALL DTRMV( 'Upper', 'Transpose', 'Non-unit',
     $                           KSIZE, S(K+1,K+1), LDS, DWORK(K1), 1 )
                     CALL DAXPY( KSIZE, ALPHA, DWORK(K1), 1, DWORK, 1 )
                  END IF
                  CALL MB04OD( 'Full', KSIZE, 0, 1, R(K+1,K+1), LDR,
     $                         DWORK, 1, DWORK, 1, DWORK, 1, DWORK(K2),
     $                         DWORK(K1) )
               END IF
            END IF
            GO TO 10
         END IF
C        END WHILE 10
C
      ELSE
C
C        Case op(M) = M'.
C
         KOUNT = N
C
   90    CONTINUE
C        WHILE( KOUNT.GE.1 )LOOP
         IF ( KOUNT.GE.1 ) THEN
            K = KOUNT
            IF ( KOUNT.EQ.1 ) THEN
               TBYT  = .FALSE.
               KOUNT = KOUNT - 1
            ELSE IF ( S(K,K-1).EQ.ZERO ) THEN
               TBYT  = .FALSE.
               KOUNT = KOUNT - 1
            ELSE
               TBYT = .TRUE.
               K = K - 1
               IF ( K.GT.1 ) THEN
                  IF ( S(K,K-1).NE.ZERO ) THEN
                     INFO = 3
                     RETURN
                  END IF
               END IF
               KOUNT = KOUNT - 2
            END IF
            IF ( TBYT ) THEN
C
C              Solve the two-by-two Lyapunov equation corresponding to
C              (6.1) or (10.19), using the routine SB03OY.
C
               B(1,1) = S(K,K)
               B(2,1) = S(K+1,K)
               B(1,2) = S(K,K+1)
               B(2,2) = S(K+1,K+1)
               U(1,1) = R(K,K)
               U(1,2) = R(K,K+1)
               U(2,2) = R(K+1,K+1)
C
               CALL SB03OY( DISCR, LTRANS, ISGN, B, 2, U, 2, A, 2,
     $                      SCALOC, INFO )
               IF ( INFO.GT.1 )
     $            RETURN
               INFOM = MAX( INFO, INFOM )
               IF( SCALOC.NE.ONE ) THEN
C
                  DO 100 J = 1, N
                     CALL DSCAL( J, SCALOC, R(1,J), 1 )
  100             CONTINUE
C
                  SCALE = SCALE*SCALOC
               END IF
               R(K,K)     = U(1,1)
               R(K,K+1)   = U(1,2)
               R(K+1,K+1) = U(2,2)
C
C              If we are not at the front of S then set up and solve
C              equation corresponding to (6.2) or (10.20).
C
C              Note that  SB03OY  returns  ( inv( u11 )*s11*u11 ) in  B
C              and returns scaled alpha, alpha = inv( u11 )*r11, in  A.
C              ksize is the order of the remainder leading part of  S.
C              k1, k2 and k3 point to the start of vectors in  DWORK.
C
               IF ( KOUNT.GE.1 ) THEN
                  KSIZE = K - 1
                  K1 = KSIZE + 1
                  K2 = KSIZE + K1
                  K3 = KSIZE + K2
C
C                 Form the right-hand side of equations corresponding to
C                 (6.2) or (10.20), the first column in DWORK( 1 ) ,...,
C                 DWORK( k - 1 ) the second in DWORK( k ) ,...,
C                 DWORK( 2*( k - 1 ) ).
C
                  CALL DCOPY( KSIZE, R(1,K), 1, DWORK, 1 )
                  CALL DCOPY( KSIZE, R(1,K+1), 1, DWORK(K1), 1 )
                  CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit',
     $                        KSIZE, 2, -ONE, A, 2, DWORK, KSIZE )
                  IF ( CONT ) THEN
                     CALL DAXPY( KSIZE, -R(K,K), S(1,K), 1, DWORK, 1 )
                     CALL DAXPY( KSIZE, -R(K,K+1), S(1,K), 1, DWORK(K1),
     $                           1)
                     CALL DAXPY( KSIZE, -R(K+1,K+1), S(1,K+1), 1,
     $                           DWORK(K1), 1 )
                  ELSE
                     CALL DAXPY( KSIZE, -( R(K,K)*B(1,1) + R(K,K+1)
     $                           *B(1,2) ), S(1,K), 1, DWORK, 1 )
                     CALL DAXPY( KSIZE, -R(K+1,K+1)*B(1,2), S(1,K+1), 1,
     $                           DWORK, 1 )
                     CALL DAXPY( KSIZE, -( R(K,K)*B(2,1) + R(K,K+1)
     $                           *B(2,2) ), S(1,K), 1, DWORK(K1), 1 )
                     CALL DAXPY( KSIZE, -R(K+1,K+1)*B(2,2), S(1,K+1), 1,
     $                           DWORK(K1), 1 )
                  END IF
C
C                 SB03OR  solves the Sylvester equations. The solution
C                 is overwritten on DWORK.
C
                  CALL SB03OR( DISCR, LTRANS, KSIZE, 2, S, LDS, B, 2,
     $                         DWORK, KSIZE, SCALOC, INFO )
                  INFOM = MAX( INFO, INFOM )
                  IF( SCALOC.NE.ONE ) THEN
C
                     DO 110 J = 1, N
                        CALL DSCAL( J, SCALOC, R(1,J), 1 )
  110                CONTINUE
C
                     SCALE = SCALE*SCALOC
                  END IF
C
C                 Copy the solution into the next  2*( k - 1 ) elements
C                 of  DWORK.
C
                  CALL DCOPY( 2*KSIZE, DWORK, 1, DWORK(K2), 1 )
C
C                 Now form the matrix  Rhat  of equation corresponding
C                 to (6.4) or (10.22) (corrected version).
C
                  IF ( CONT ) THEN
C
C                    Swap the two columns of R with DWORK.
C
                     CALL DSWAP( KSIZE, DWORK, 1, R(1,K), 1 )
                     CALL DSWAP( KSIZE, DWORK(K1), 1, R(1,K+1), 1 )
C
C                    1st column:
C
                     CALL DAXPY( KSIZE, -A(1,1), DWORK(K2), 1, DWORK,
     $                           1 )
C
C                    2nd column:
C
                     CALL DAXPY( KSIZE, -A(1,2), DWORK(K2), 1,
     $                           DWORK(K1), 1 )
                     CALL DAXPY( KSIZE, -A(2,2), DWORK(K3), 1,
     $                           DWORK(K1), 1 )
                  ELSE
C
C                    Form  v = S1*u + s*u11, overwriting  v  on DWORK.
C
C                    Compute  S1*u,  first multiplying by the triangular
C                    part of  S1.
C
                     CALL DTRMM( 'Left', 'Upper', 'No transpose',
     $                           'Non-unit', KSIZE, 2, ONE, S, LDS,
     $                           DWORK, KSIZE )
C
C                    Then multiply by the subdiagonal of  S1  and add in
C                    to the above result.
C
                     J1 = K1
C
                     DO 120 J = 2, KSIZE
                        J1 = J1 + 1
                        IF ( S(J,J-1).NE.ZERO ) THEN
                           DWORK(J)  = S(J,J-1)*DWORK(K2+J-2) + DWORK(J)
                           DWORK(J1) = S(J,J-1)*DWORK(K3+J-2) +
     $                                 DWORK(J1)
                        END IF
  120                CONTINUE
C
C                    Add in s*u11.
C
                     CALL DAXPY( KSIZE, R(K,K), S(1,K), 1, DWORK, 1 )
                     CALL DAXPY( KSIZE, R(K,K+1), S(1,K), 1, DWORK(K1),
     $                           1 )
                     CALL DAXPY( KSIZE, R(K+1,K+1), S(1,K+1), 1,
     $                           DWORK(K1), 1 )
C
C                    Next recover r from R, swapping r with u.
C
                     CALL DSWAP( KSIZE, DWORK(K2), 1, R(1,K), 1 )
                     CALL DSWAP( KSIZE, DWORK(K3), 1, R(1,K+1), 1 )
C
C                    Now we perform the QL factorization.
C
C                    ( a' ) = Q*( t ),
C                    ( b' )
C
C                    and form
C
C                    ( p' ) = Q'*( r' ).
C                    ( y' )      ( v' )
C
C                    y  is then the correct vector to use in the
C                    relation corresponding to (10.22).
C                    Note that  a  is upper triangular and that  t  and
C                    p  are not required.
C
                     CALL DLARFG( 3, A(2,2), B(2,1), 2, TAU1 )
                     V1  = B(2,1)
                     T1  = TAU1*V1
                     V2  = B(2,2)
                     T2  = TAU1*V2
                     SUM = A(1,2) + V1*B(1,1) + V2*B(1,2)
                     B(1,1) = B(1,1) - SUM*T1
                     B(1,2) = B(1,2) - SUM*T2
                     CALL DLARFG( 3, A(1,1), B(1,1), 2, TAU2 )
                     V3 = B(1,1)
                     T3 = TAU2*V3
                     V4 = B(1,2)
                     T4 = TAU2*V4
                     J1 = K1
                     J2 = K2
                     J3 = K3
C
                     DO 130 J = 1, KSIZE
                        SUM = DWORK(J3) + V1*DWORK(J) + V2*DWORK(J1)
                        D1  = DWORK(J)  - SUM*T1
                        D2  = DWORK(J1) - SUM*T2
                        SUM = DWORK(J2) + V3*D1 + V4*D2
                        DWORK(J)  =  D1 - SUM*T3
                        DWORK(J1) =  D2 - SUM*T4
                        J1 = J1 + 1
                        J2 = J2 + 1
                        J3 = J3 + 1
  130                CONTINUE
C
                  END IF
C
C                 Now update  R1  to give  Rhat.
C
                  CALL MB04ND( 'Full', KSIZE, 0, 2, R, LDR, DWORK,
     $                         KSIZE, DWORK, 1, DWORK, 1, DWORK(K2),
     $                         DWORK(K3) )
               END IF
            ELSE
C
C              1-by-1 block.
C
C              Make sure S is stable or convergent and find u11 in
C              equation corresponding to (5.13) or (10.15).
C
               IF ( DISCR ) THEN
                  ABSSKK = ABS( S(K,K) )
                  IF ( ( ABSSKK - ONE ).GE.ZERO ) THEN
                     INFO = 2
                     RETURN
                  END IF
                  TEMP = SQRT( ( ONE - ABSSKK )*( ONE + ABSSKK ) )
               ELSE
                  IF ( S(K,K).GE.ZERO ) THEN
                     INFO = 2
                     RETURN
                  END IF
                  TEMP = SQRT( ABS( TWO*S(K,K) ) )
               END IF
C
               SCALOC = ONE
               IF( TEMP.LT.SMIN ) THEN
                  TEMP  = SMIN
                  INFOM = 1
               END IF
               DR = ABS( R(K,K) )
               IF( TEMP.LT.ONE .AND. DR.GT.ONE ) THEN
                  IF( DR.GT.BIGNUM*TEMP )
     $               SCALOC = ONE / DR
               END IF
               ALPHA = SIGN( TEMP, R(K,K) )
               R(K,K) = R(K,K)/ALPHA
               IF( SCALOC.NE.ONE ) THEN
C
                  DO 140 J = 1, N
                     CALL DSCAL( J, SCALOC, R(1,J), 1 )
  140             CONTINUE
C
                  SCALE = SCALE*SCALOC
               END IF
C
C              If we are not at the front of  S  then set up and solve
C              equation corresponding to (5.14) or (10.16).  ksize is
C              the order of the remainder leading part of  S.  k1 and k2
C              point to the start of vectors in  DWORK.
C
               IF ( KOUNT.GE.1 ) THEN
                  KSIZE = K - 1
                  K1 = KSIZE + 1
                  K2 = KSIZE + K1
C
C                 Form the right-hand side in DWORK( 1 ),...,
C                 DWORK( k - 1 ).
C
                  CALL DCOPY( KSIZE, R(1,K), 1, DWORK, 1 )
                  CALL DSCAL( KSIZE, -ALPHA, DWORK, 1 )
                  IF ( CONT ) THEN
                     CALL DAXPY( KSIZE, -R(K,K), S(1,K), 1, DWORK, 1 )
                  ELSE
                     CALL DAXPY( KSIZE, -S(K,K)*R(K,K), S(1,K), 1,
     $                          DWORK, 1 )
                  END IF
C
C                 SB03OR solves the Sylvester equations. The solution is
C                 overwritten on  DWORK.
C
                  CALL SB03OR( DISCR, LTRANS, KSIZE, 1, S, LDS, S(K,K),
     $                         1, DWORK, KSIZE, SCALOC, INFO )
                  INFOM = MAX( INFO, INFOM )
                  IF( SCALOC.NE.ONE ) THEN
C
                     DO 150 J = 1, N
                        CALL DSCAL( J, SCALOC, R(1,J), 1 )
  150                CONTINUE
C
                     SCALE = SCALE*SCALOC
                  END IF
C
C                 Copy the solution into the next  ( k - 1 ) elements
C                 of  DWORK,  copy the solution back into  R  and copy
C                 the column of  R  back into  DWORK.
C
                  CALL DCOPY( KSIZE, DWORK, 1, DWORK(K1), 1 )
                  CALL DSWAP( KSIZE, DWORK, 1, R(1,K), 1 )
C
C                 Now form the matrix  Rhat  of equation corresponding
C                 to (5.15) or (10.17), first computing  y  in  DWORK,
C                 and then updating  R1.
C
                  IF ( CONT ) THEN
                     CALL DAXPY( KSIZE, -ALPHA, DWORK(K1), 1, DWORK, 1 )
                  ELSE
C
C                    First form  lambda( 1 )*r  and then add in
C                    alpha*u11*s.
C
                     CALL DSCAL( KSIZE, -S(K,K), DWORK, 1 )
                     CALL DAXPY( KSIZE, ALPHA*R(K,K), S(1,K), 1, DWORK,
     $                           1 )
C
C                    Now form  alpha*S1*u,  first multiplying by the
C                    sub-diagonal of  S1  and then the triangular part
C                    of  S1,  and add the result in DWORK.
C
                     DO 160 J = 2, KSIZE
                        IF ( S(J,J-1).NE.ZERO ) DWORK(J)
     $                         = ALPHA*S(J,J-1)*DWORK(K1+J-2) + DWORK(J)
  160                CONTINUE
C
                     CALL DTRMV( 'Upper', 'No transpose', 'Non-unit',
     $                           KSIZE, S, LDS, DWORK(K1), 1 )
                     CALL DAXPY( KSIZE, ALPHA, DWORK(K1), 1, DWORK, 1 )
                  END IF
                  CALL MB04ND( 'Full', KSIZE, 0, 1, R, LDR, DWORK,
     $                         KSIZE, DWORK, 1, DWORK, 1, DWORK(K2),
     $                         DWORK(K1) )
               END IF
            END IF
            GO TO 90
         END IF
C        END WHILE 90
C
      END IF
      INFO = INFOM
      RETURN
C *** Last line of SB03OT ***
      END