text:       SUBROUTINE SB01FY( DISCR, N, M, A, LDA, B, LDB, F, LDF, V, LDV,
       $                   INFO )
  C
  C     SLICOT RELEASE 5.5.
  C
  C     Copyright (c) 2002-2012 NICONET e.V.
  C
  C     PURPOSE
  C
  C     To compute the inner denominator of a right-coprime factorization
  C     of a system of order N, where N is either 1 or 2. Specifically,
  C     given the N-by-N unstable system state matrix A and the N-by-M
  C     system input matrix B, an M-by-N state-feedback matrix F and
  C     an M-by-M matrix V are constructed, such that the system
  C     (A + B*F, B*V, F, V) is inner.
  C
  C     ARGUMENTS
  C
  C     Mode Parameters
  C
  C     DISCR   LOGICAL
  C             Specifies the type of system as follows:
  C             = .FALSE.:  continuous-time system;
  C             = .TRUE. :  discrete-time system.
  C
  C     Input/Output Parameters
  C
  C     N       (input) INTEGER
  C             The order of the matrix A and also the number of rows of
  C             the matrix B and the number of columns of the matrix F.
  C             N is either 1 or 2.
  C
  C     M       (input) INTEGER
  C             The number of columns of the matrices B and V, and also
  C             the number of rows of the matrix F.  M >= 0.
  C
  C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
  C             The leading N-by-N part of this array must contain the
  C             system state matrix A whose eigenvalues must have positive
  C             real parts if DISCR = .FALSE. or moduli greater than unity
  C             if DISCR = .TRUE..
  C
  C     LDA     INTEGER
  C             The leading dimension of array A.  LDA >= N.
  C
  C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
  C             The leading N-by-M part of this array must contain the
  C             system input matrix B.
  C
  C     LDB     INTEGER
  C             The leading dimension of array B.  LDB >= N.
  C
  C     F       (output) DOUBLE PRECISION array, dimension (LDF,N)
  C             The leading M-by-N part of this array contains the state-
  C             feedback matrix F which assigns one eigenvalue (if N = 1)
  C             or two eigenvalues (if N = 2) of the matrix A + B*F in
  C             symmetric positions with respect to the imaginary axis
  C             (if DISCR = .FALSE.) or the unit circle (if
  C             DISCR = .TRUE.).
  C
  C     LDF     INTEGER
  C             The leading dimension of array F.  LDF >= MAX(1,M).
  C
  C     V       (output) DOUBLE PRECISION array, dimension (LDV,M)
  C             The leading M-by-M upper triangular part of this array
  C             contains the input/output matrix V of the resulting inner
  C             system in upper triangular form.
  C             If DISCR = .FALSE., the resulting V is an identity matrix.
  C
  C     LDV     INTEGER
  C             The leading dimension of array V.  LDF >= MAX(1,M).
  C
  C     Error Indicator
  C
  C     INFO    INTEGER
  C             = 0:  successful exit;
  C             = 1:  if uncontrollability of the pair (A,B) is detected;
  C             = 2:  if A is stable or at the stability limit;
  C             = 3:  if N = 2 and A has a pair of real eigenvalues.
  C
  C     CONTRIBUTOR
  C
  C     A. Varga, German Aerospace Center,
  C     DLR Oberpfaffenhofen, July 1998.
  C     Based on the RASP routine RCFID2.
  C
  C     REVISIONS
  C
  C     Nov. 1998, V. Sima, Research Institute for Informatics, Bucharest.
  C     Dec. 1998, V. Sima, Katholieke Univ. Leuven, Leuven.
  C     Feb. 1999, A. Varga, DLR Oberpfaffenhofen.
  C
  C     KEYWORDS
  C
  C     Coprime factorization, eigenvalue, eigenvalue assignment,
  C     feedback control, pole placement, state-space model.
  C
  C     ******************************************************************
  C
  C     .. Parameters ..
        DOUBLE PRECISION  ONE, TWO, ZERO
        PARAMETER         ( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0 )
  C     .. Scalar Arguments ..
        LOGICAL           DISCR
        INTEGER           INFO, LDA, LDB, LDF, LDV, M, N
  C     .. Array Arguments ..
        DOUBLE PRECISION  A(LDA,*), B(LDB,*), F(LDF,*), V(LDV,*)
  C     .. Local Scalars ..
        INTEGER           I
        DOUBLE PRECISION  CS, R11, R12, R22, SCALE, SN, TEMP
  C     .. Local Arrays ..
        DOUBLE PRECISION  AT(2,2), DUMMY(2,2), U(2,2)
  C     .. External Functions ..
        DOUBLE PRECISION  DLAPY2, DLAPY3
        EXTERNAL          DLAPY2, DLAPY3
  C     .. External Subroutines ..
        EXTERNAL          DLARFG, DLASET, DLATZM, DROTG, DTRTRI, MA02AD,
       $                  MB04OX, SB03OY
  C     .. Intrinsic Functions ..
        INTRINSIC         ABS, SQRT
  C     .. Executable Statements ..
  C
  C     For efficiency reasons, the parameters are not checked.
  C
        INFO = 0
  C
  C     Compute an N-by-N upper triangular R such that R'*R = B*B' and
  C     find an upper triangular matrix U in the equation
  C
  C     A'*U'*U + U'*U*A = R'*R if DISCR = .FALSE. or
  C     A'*U'*U*A - U'*U = R'*R if DISCR = .TRUE. .
  C
        CALL MA02AD( 'Full', N, M, B, LDB, F, LDF )
  C
        IF( N.EQ.1 ) THEN
  C
  C        The N = 1 case.
  C
           IF( M.GT.1 )
       $      CALL DLARFG( M, F(1,1), F(2,1), 1, TEMP )
           R11 = ABS( F(1,1) )
  C
  C        Make sure A is unstable or divergent and find U.
  C
           IF( DISCR ) THEN
              TEMP = ABS( A(1,1) )
              IF( TEMP.LE.ONE ) THEN
                 INFO = 2
                 RETURN
              ELSE
                 TEMP = R11 / SQRT( ( TEMP - ONE )*( TEMP + ONE ) )
              END IF
           ELSE
              IF( A(1,1).LE.ZERO ) THEN
                 INFO = 2
                 RETURN
              ELSE
                 TEMP = R11 / SQRT( ABS( TWO*A(1,1) ) )
              END IF
           END IF
           U(1,1) = TEMP
           SCALE  = ONE
        ELSE
  C
  C        The N = 2 case.
  C
           IF( M.GT.1 ) THEN
              CALL DLARFG( M, F(1,1), F(2,1), 1, TEMP )
              CALL DLATZM( 'Left', M, N-1, F(2,1), 1, TEMP, F(1,2),
       $                   F(2,2), LDF, V )
           END IF
           R11 = F(1,1)
           R12 = F(1,2)
           IF( M.GT.2 )
       $      CALL DLARFG( M-1, F(2,2), F(3,2), 1, TEMP )
           IF( M.EQ.1 ) THEN
              R22 = ZERO
           ELSE
              R22 = F(2,2)
           END IF
           AT(1,1) = A(1,1)
           AT(1,2) = A(2,1)
           AT(2,1) = A(1,2)
           AT(2,2) = A(2,2)
           U(1,1)  = R11
           U(1,2)  = R12
           U(2,2)  = R22
           CALL SB03OY( DISCR, .FALSE., -1, AT, 2, U, 2, DUMMY, 2,
       $                SCALE, INFO )
           IF( INFO.NE.0 ) THEN
              IF( INFO.NE.4 ) THEN
                 INFO = 2
              ELSE
                 INFO = 3
              END IF
              RETURN
           END IF
        END IF
  C
  C     Check the controllability of the pair (A,B).
  C
  C     Warning. Only an exact controllability check is performed.
  C              If the pair (A,B) is nearly uncontrollable, then
  C              the computed results may be inaccurate.
  C
        DO 10 I = 1, N
           IF( U(I,I).EQ.ZERO ) THEN
              INFO = 1
              RETURN
           END IF
     10 CONTINUE
  C
  C     Set V = I.
  C
        CALL DLASET( 'Upper', M, M, ZERO, ONE, V, LDV )
  C
        IF( DISCR ) THEN
  C
  C        Compute an upper triangular matrix V such that
  C                                 -1
  C        V*V' = (I+B'*inv(U'*U)*B)  .
  C
  C        First compute F = B'*inv(U) and the Cholesky factorization
  C        of I + F*F'.
  C
           DO 20 I = 1, M
              F(I,1) = B(1,I)/U(1,1)*SCALE
     20    CONTINUE
           IF( N.EQ.2 ) THEN
              DO 30 I = 1, M
                 F(I,2) = ( B(2,I) - F(I,1)*U(1,2) )/U(2,2)*SCALE
     30       CONTINUE
              CALL MB04OX( M, V, LDV, F(1,2), 1 )
           END IF
           CALL MB04OX( M, V, LDV, F(1,1), 1 )
           CALL DTRTRI( 'Upper', 'NonUnit', M, V, LDV, INFO )
        END IF
  C
  C     Compute the feedback matrix F as:
  C
  C     1)   If DISCR = .FALSE.
  C
  C             F = -B'*inv(U'*U);
  C
  C     2)   If DISCR = .TRUE.
  C                                -1
  C             F = -B'*(U'*U+B*B')  *A.
  C
        IF( N.EQ.1 ) THEN
           IF( DISCR ) THEN
              TEMP = -A(1,1)
              R11  = DLAPY2( U(1,1), R11 )
              DO 40 I = 1, M
                 F(I,1) = ( ( B(1,I)/R11 )/R11 )*TEMP
     40       CONTINUE
           ELSE
              R11 = U(1,1)
              DO 50 I = 1, M
                 F(I,1) = -( ( B(1,I)/R11 )/R11 )
     50       CONTINUE
           END IF
        ELSE
  C
  C        Set R = U  if DISCR = .FALSE. or compute the Cholesky
  C        factorization of R'*R = U'*U+B*B' if DISCR = .TRUE..
  C
           IF( DISCR ) THEN
              TEMP = U(1,1)
              CALL DROTG( R11, TEMP, CS, SN )
              TEMP = -SN*R12 + CS*U(1,2)
              R12  =  CS*R12 + SN*U(1,2)
              R22  = DLAPY3( R22, TEMP, U(2,2) )
           ELSE
              R11 = U(1,1)
              R12 = U(1,2)
              R22 = U(2,2)
           END IF
  C
  C        Compute F = -B'*inv(R'*R).
  C
           DO 60 I = 1, M
              F(I,1) = -B(1,I)/R11
              F(I,2) = -( B(2,I) + F(I,1)*R12 )/R22
              F(I,2) =  F(I,2)/R22
              F(I,1) =  ( F(I,1) - F(I,2)*R12 )/R11
     60    CONTINUE
           IF( DISCR ) THEN
  C
  C           Compute F <-- F*A.
  C
              DO 70 I = 1, M
                 TEMP   = F(I,1)*A(1,1) + F(I,2)*A(2,1)
                 F(I,2) = F(I,1)*A(1,2) + F(I,2)*A(2,2)
                 F(I,1) = TEMP
     70       CONTINUE
           END IF
        END IF
  C
        RETURN
  C *** Last line of SB01FY ***
        END
  