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

      SUBROUTINE AB01OD( STAGES, JOBU, JOBV, N, M, A, LDA, B, LDB, U, 
     $                   LDU, V, LDV, NCONT, INDCON, KSTAIR, TOL, IWORK,
     $                   DWORK, LDWORK, INFO )
C
C     RELEASE 4.0, WGS COPYRIGHT 1999.
C
C     PURPOSE
C
C     To reduce the matrices A and B using (and optionally accumulating)
C     state-space and input-space transformations U and V respectively,
C     such that the pair of matrices
C 
C        Ac = U' * A * U,    Bc = U' * B * V 
C
C     are in upper "staircase" form. Specifically, 
C
C             [ Acont     *    ]         [ Bcont ] 
C        Ac = [                ],   Bc = [       ], 
C             [   0    Auncont ]         [   0   ] 
C
C        and
C
C                [ A11 A12  . . .  A1,p-1 A1p ]         [ B1 ] 
C                [ A21 A22  . . .  A2,p-1 A2p ]         [ 0  ] 
C                [  0  A32  . . .  A3,p-1 A3p ]         [ 0  ] 
C        Acont = [  .   .   . . .    .     .  ],   Bc = [ .  ], 
C                [  .   .     . .    .     .  ]         [ .  ] 
C                [  .   .       .    .     .  ]         [ .  ] 
C                [  0   0   . . .  Ap,p-1 App ]         [ 0  ] 
C
C     where the blocks  B1, A21, ..., Ap,p-1  have full row ranks and 
C     p is the controllability index of the pair.  The size of the
C     block Auncont is equal to the dimension of the uncontrollable
C     subspace of the pair (A, B).  The first stage of the reduction,
C     the "forward" stage, accomplishes the reduction to the orthogonal 
C     canonical form (see SLICOT library routine AB01ND). The blocks 
C     B1, A21, ..., Ap,p-1 are further reduced in a second, "backward" 
C     stage to upper triangular form using RQ factorization. Each of
C     these stages is optional.
C
C     ARGUMENTS 
C
C     Mode Parameters
C
C     STAGES  CHARACTER*1
C             Specifies the reduction stages to be performed as follows:
C             = 'F':  Perform the forward stage only;
C             = 'B':  Perform the backward stage only;
C             = 'A':  Perform both (all) stages.
C
C     JOBU    CHARACTER*1
C             Indicates whether the user wishes to accumulate in a
C             matrix U the state-space transformations as follows:
C             = 'N':  Do not form U; 
C             = 'I':  U is internally initialized to the unit matrix (if
C                     STAGES <> 'B'), or updated (if STAGES = 'B'), and
C                     the orthogonal transformation matrix U is 
C                     returned.
C
C     JOBV    CHARACTER*1
C             Indicates whether the user wishes to accumulate in a 
C             matrix V the input-space transformations as follows:
C             = 'N':  Do not form V; 
C             = 'I':  V is initialized to the unit matrix and the
C                     orthogonal transformation matrix V is returned.
C             JOBV is not referenced if STAGES = 'F'. 
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The actual state dimension, i.e. the order of the
C             matrix A.  N >= 0. 
C
C     M       (input) INTEGER
C             The actual input dimension.  M >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must 
C             contain the state transition matrix A to be transformed.
C             If STAGES = 'B', A should be in the orthogonal canonical  
C             form, as returned by SLICOT library routine AB01ND.
C             On exit, the leading N-by-N part of this array contains 
C             the transformed state transition matrix U' * A * U.
C             The leading NCONT-by-NCONT part contains the upper block 
C             Hessenberg state matrix Acont in Ac, given by U' * A * U, 
C             of a controllable realization for the original system.
C             The elements below the first block-subdiagonal are set to 
C             zero.  If STAGES <> 'F', the subdiagonal blocks of A are 
C             triangularized by RQ factorization, and the annihilated 
C             elements are explicitly zeroed.  
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
C             On entry, the leading N-by-M part of this array must 
C             contain the input matrix B to be transformed.
C             If STAGES = 'B', B should be in the orthogonal canonical 
C             form, as returned by SLICOT library routine AB01ND.
C             On exit with STAGES = 'F', the leading N-by-M part of 
C             this array contains the transformed input matrix U' * B, 
C             with all elements but the first block set to zero.
C             On exit with STAGES <> 'F', the leading N-by-M part of 
C             this array contains the transformed input matrix 
C             U' * B * V, with all elements but the first block set to
C             zero and the first block in upper triangular form.
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     U       (input/output) DOUBLE PRECISION array, dimension (LDU,N)
C             If STAGES <> 'B' or JOBU = 'N', then U need not be set
C             on entry.
C             If STAGES = 'B' and JOBU = 'I', then, on entry, the
C             leading N-by-N part of this array must contain the 
C             transformation matrix U that reduced the pair to the 
C             orthogonal canonical form.
C             On exit, if JOBU = 'I', the leading N-by-N part of this 
C             array contains the transformation matrix U that performed
C             the specified reduction.
C             If JOBU = 'N', the array U is not referenced and can be 
C             supplied as a dummy array (i.e. set parameter LDU = 1 and 
C             declare this array to be U(1,1) in the calling program).
C
C     LDU     INTEGER
C             The leading dimension of array U. 
C             If JOBU = 'I', LDU >= MAX(1,N);  if JOBU = 'N', LDU >= 1.
C             
C     V       (output) DOUBLE PRECISION array, dimension (LDV,M)
C             If JOBV = 'I', then the leading M-by-M part of this array
C             contains the transformation matrix V.
C             If STAGES = 'F', or JOBV = 'N', the array V is not 
C             referenced and can be supplied as a dummy array (i.e. set 
C             parameter  LDV = 1 and declare this array to be V(1,1) in 
C             the calling program).
C
C     LDV     INTEGER
C             The leading dimension of array V. 
C             If STAGES <> 'F' and JOBV = 'I', LDV >= MAX(1,M);  
C             if STAGES = 'F' or JOBV = 'N', LDV >= 1.
C
C     NCONT   (input/output) INTEGER
C             The order of the controllable state-space representation.
C             NCONT is input only if STAGES = 'B'.
C
C     INDCON  (input/output) INTEGER
C             The number of stairs in the staircase form (also, the
C             controllability index of the controllable part of the
C             system representation).
C             INDCON is input only if STAGES = 'B'.
C
C     KSTAIR  (input/output) INTEGER array, dimension (N)
C             The leading INDCON elements of this array contain the
C             dimensions of the stairs, or, also, the orders of the
C             diagonal blocks of Acont.
C             KSTAIR is input if STAGES = 'B', and output otherwise.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used in rank determination when 
C             transforming (A, B). If the user sets TOL > 0, then
C             the given value of TOL is used as a lower bound for the
C             reciprocal condition number (see the description of the
C             argument RCOND in the SLICOT routine MB03OD);  a 
C             (sub)matrix whose estimated condition number is less than
C             1/TOL is considered to be of full rank.  If the user sets
C             TOL <= 0, then an implicitly computed, default tolerance,
C             defined by  TOLDEF = N*N*EPS,  is used instead, where EPS 
C             is the machine precision (see LAPACK Library routine 
C             DLAMCH).
C             TOL is not referenced if STAGES = 'B'.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (M)
C             IWORK is not referenced if STAGES = 'B'.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK. 
C             If STAGES <> 'B', LDWORK >= MAX(1, N + MAX(N,3*M));
C             If STAGES =  'B', LDWORK >= MAX(1, M + MAX(N,M)).
C             For optimum performance LDWORK should be larger.
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
C     METHOD
C
C     Staircase reduction of the pencil [B|sI - A] is used. Orthogonal
C     transformations U and V are constructed such that
C
C
C                        |B |sI-A      *  . . .  *      *       |
C                        | 1|    11       .      .      .       |
C                        |  |  A    sI-A    .    .      .       |
C                        |  |   21      22    .  .      .       |
C                        |  |        .     .     *      *       |
C     [U'BV|sI - U'AU] = |0 |     0    .     .                  |
C                        |  |            A     sI-A     *       |
C                        |  |             p,p-1    pp           |
C                        |  |                                   |
C                        |0 |         0          0   sI-A       |
C                        |  |                            p+1,p+1|
C
C
C     where the i-th diagonal block of U'AU has dimension KSTAIR(i),
C     for i = 1,...,p. The value of p is returned in INDCON. The last
C     block contains the uncontrollable modes of the (A,B)-pair which
C     are also the generalized eigenvalues of the above pencil.
C
C     The complete reduction is performed in two stages. The first,
C     forward stage accomplishes the reduction to the orthogonal 
C     canonical form. The second, backward stage consists in further
C     reduction to triangular form by applying left and right orthogonal
C     transformations.
C
C     REFERENCES
C
C     [1] Van Dooren, P.
C         The generalized eigenvalue problem in linear system theory.
C         IEEE Trans. Auto. Contr., AC-26, pp. 111-129, 1981.
C
C     [2] Miminis, G. and Paige, C.
C         An algorithm for pole assignment of time-invariant multi-input
C         linear systems.
C         Proc. 21st IEEE CDC, Orlando, Florida, 1, pp. 62-67, 1982.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires O((N + M) x N**2) operations and is
C     backward stable (see [1]).
C
C     FURTHER COMMENTS
C
C     If the system matrices A and B are badly scaled, it would be
C     useful to scale them with SLICOT routine TB01ID, before calling 
C     the routine.
C 
C     CONTRIBUTOR
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996.
C     Supersedes Release 2.0 routine AB01CD by M. Vanbegin, and
C     P. Van Dooren, Philips Research Laboratory, Brussels, Belgium.
C
C     REVISIONS
C
C     January 14, 1997, February 12, 1998.
C
C     KEYWORDS
C
C     Controllability, generalized eigenvalue problem, orthogonal
C     transformation, staircase form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER*1       JOBU, JOBV, STAGES
      INTEGER           INFO, INDCON, LDA, LDB, LDU, LDV, LDWORK, M, N,
     $                  NCONT
      DOUBLE PRECISION  TOL
C     .. Array Arguments ..
      INTEGER           IWORK(*), KSTAIR(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), U(LDU,*), V(LDV,*)
C     .. Local Scalars ..
      LOGICAL           LJOBUI, LJOBVI, LSTAGB, LSTGAB
      INTEGER           I, I0, IBSTEP, ITAU, J0, JINI, JWORK, MCRT, MM, 
     $                  NCRT, WRKOPT
C     .. External Functions ..
      LOGICAL           LSAME
      EXTERNAL          LSAME
C     .. External Subroutines ..
      EXTERNAL          AB01ND, DGERQF, DLACPY, DLASET, DORGRQ, DORMRQ,
     $                  DSWAP, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         INT, MAX, MIN
C     .. Executable Statements ..
C
      INFO = 0
      LJOBUI = LSAME( JOBU, 'I' )
C
      LSTAGB = LSAME( STAGES, 'B' )
      LSTGAB = LSAME( STAGES, 'A' ).OR.LSTAGB
C
      IF ( LSTGAB ) THEN
         LJOBVI = LSAME( JOBV, 'I' )
      END IF
C
C     Test the input scalar arguments.
C
      IF( .NOT.LSTGAB .AND. .NOT.LSAME( STAGES, 'F' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.LJOBUI .AND. .NOT.LSAME( JOBU, 'N' ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( M.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( .NOT.LJOBUI .AND. LDU.LT.1 .OR.
     $              LJOBUI .AND. LDU.LT.MAX( 1, N ) ) THEN
         INFO = -11
      ELSE IF( .NOT.LSTAGB .AND. LDWORK.LT.MAX( 1, N + MAX( N, 3*M ) )
     $         .OR. LSTAGB .AND. LDWORK.LT.MAX( 1, M + MAX( N, M ) ) )
     $         THEN 
         INFO = -20
      ELSE IF( LSTGAB ) THEN
         IF( .NOT.LJOBVI .AND. .NOT.LSAME( JOBV, 'N' ) ) THEN
            INFO = -3
         ELSE IF( .NOT.LJOBVI .AND. LDV.LT.1 .OR.
     $                 LJOBVI .AND. LDV.LT.MAX( 1, M ) ) THEN
            INFO = -13
         END IF
      ELSE IF( .NOT.LSTAGB .AND. (TOL.LT.ZERO .OR.  TOL.GT.ONE) ) THEN
C added by S. STEER (see mb03oy)
         INFO = -17
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'AB01OD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( MIN( N, M ).EQ.0 ) THEN
         NCONT  = 0
         INDCON = 0
         RETURN
      END IF
C
C     (Note: Comments in the code beginning "Workspace:" describe the
C     minimal amount of real workspace needed at that point in the
C     code, as well as the preferred amount for good performance.
C     NB refers to the optimal block size for the immediately
C     following subroutine, as returned by ILAENV.)
C
      ITAU = 1
      WRKOPT = 1
C 
      IF ( .NOT.LSTAGB ) THEN
C
C        Perform the forward stage computations of the staircase
C        algorithm on B and A: reduce the (A, B) pair to orthogonal
C        canonical form.
C
C        Workspace: N + MAX(N,3*M).
C        
         JWORK = N + 1 
         CALL AB01ND( JOBU, N, M, A, LDA, B, LDB, NCONT, INDCON,
     $                KSTAIR, U, LDU, DWORK(ITAU), TOL, IWORK,
     $                DWORK(JWORK), LDWORK-JWORK+1, INFO )
         IF(INFO.LT.0) RETURN
C
         WRKOPT = INT( DWORK(JWORK) ) + JWORK - 1
      END IF
C
C     Exit if no further reduction to triangularize B1 and subdiagonal
C     blocks of A is required, or if the order of the controllable part
C     is 0.
C
      IF ( .NOT.LSTGAB ) THEN 
         RETURN
      ELSE IF ( NCONT.EQ.0 .OR. INDCON.EQ.0 ) THEN
         IF( LJOBVI ) 
     $      CALL DLASET( 'F', M, M, ZERO, ONE, V, LDV )
         RETURN
      END IF
C
C     Now perform the backward steps except the last one.
C
      MCRT = KSTAIR(INDCON)
      I0 = NCONT - MCRT + 1
      JWORK = M + 1 
C     
      DO 10 IBSTEP = INDCON, 2, -1
         NCRT = KSTAIR(IBSTEP-1)
         J0 = I0 - NCRT
         MM = MIN( NCRT, MCRT )
C     
C        Compute the RQ factorization of the current subdiagonal block 
C        of A, Ai,i-1 = R*Q (where i is IBSTEP), of dimension
C        MCRT-by-NCRT, starting in position (I0,J0).
C        The matrix Q' should postmultiply U, if required.
C        Workspace: need   M + MCRT;  
C                   prefer M + MCRT*NB.
C
         CALL DGERQF( MCRT, NCRT, A(I0,J0), LDA, DWORK(ITAU),
     $                DWORK(JWORK), LDWORK-JWORK+1, INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
C        Set JINI to the first column number in A where the current 
C        transformation Q is to be applied, taking the block Hessenberg
C        form into account.
C
         IF ( IBSTEP.GT.2 ) THEN 
            JINI = J0 - KSTAIR(IBSTEP-2)
         ELSE
            JINI = 1
C
C           Premultiply the first block row (B1) of B by Q.
C           Workspace: need   2*M;
C                      prefer M + M*NB.
C     
            CALL DORMRQ( 'Left', 'No transpose', NCRT, M, MM, A(I0,J0),
     $                   LDA, DWORK(ITAU), B, LDB, DWORK(JWORK),
     $                   LDWORK-JWORK+1, INFO )
            WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
         END IF
C
C        Premultiply the appropriate block row of A by Q.
C        Workspace: need   M + N;
C                   prefer M + N*NB.
C     
         CALL DORMRQ( 'Left', 'No transpose', NCRT, N-JINI+1, MM, 
     $                A(I0,J0), LDA, DWORK(ITAU), A(J0,JINI), LDA, 
     $                DWORK(JWORK), LDWORK-JWORK+1, INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
C        Postmultiply the appropriate block column of A by Q'.
C        Workspace: need   M +  I0-1;
C                   prefer M + (I0-1)*NB.
C        
         CALL DORMRQ( 'Right', 'Transpose', I0-1, NCRT, MM, A(I0,J0),
     $                LDA, DWORK(ITAU), A(1,J0), LDA, DWORK(JWORK),
     $                LDWORK-JWORK+1, INFO )
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
         IF ( LJOBUI ) THEN
C     
C           Update U, postmultiplying it by Q'.
C           Workspace: need   M + N;
C                      prefer M + N*NB.
C     
            CALL DORMRQ( 'Right', 'Transpose', N, NCRT, MM, A(I0,J0), 
     $                   LDA, DWORK(ITAU), U(1,J0), LDU, DWORK(JWORK),
     $                   LDWORK-JWORK+1, INFO )
            WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
         END IF
C
C        Zero the subdiagonal elements of the current subdiagonal block 
C        of A.
C     
         CALL DLASET( 'F', MCRT, NCRT-MCRT, ZERO, ZERO, A(I0,J0), LDA ) 
         IF ( I0.LT.N )
     $      CALL DLASET( 'L', MCRT-1, MCRT-1, ZERO, ZERO,
     $                   A(I0+1,I0-MCRT), LDA ) 
C     
         MCRT = NCRT
         I0 = J0
C
   10 CONTINUE
C     
C     Now perform the last backward step on B, V = Qb'.
C     
C     Compute the RQ factorization of the first block of B, B1 = R*Qb.
C     Workspace: need   M + MCRT;
C                prefer M + MCRT*NB.
C
      CALL DGERQF( MCRT, M, B, LDB, DWORK(ITAU), DWORK(JWORK),
     $             LDWORK-JWORK+1, INFO )
      WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
C
      IF ( LJOBVI ) THEN
C
C        Accumulate the input-space transformations V.
C        Workspace: need 2*M;  prefer M + M*NB.
C     
         CALL DLACPY( 'F', MCRT, M-MCRT, B, LDB, V(M-MCRT+1,1), LDV )
         IF ( MCRT.GT.1 )
     $      CALL DLACPY( 'L', MCRT-1, MCRT-1, B(2,M-MCRT+1), LDB,
     $                   V(M-MCRT+2,M-MCRT+1), LDV )
         CALL DORGRQ( M, M, MCRT, V, LDV, DWORK(ITAU), DWORK(JWORK),
     $                LDWORK-JWORK+1, INFO )
C
         DO 20 I = 2, M
            CALL DSWAP( I-1, V(I, 1), LDV, V(1,I), 1 )
 20      CONTINUE    
C
         WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
      END IF
C     
C     Zero the subdiagonal elements of the submatrix B1.
C     
      CALL DLASET( 'F', MCRT, M-MCRT, ZERO, ZERO, B, LDB ) 
      IF ( MCRT.GT.1 )
     $   CALL DLASET( 'L', MCRT-1, MCRT-1, ZERO, ZERO, B(2,M-MCRT+1),
     $                LDB )
C
C     Set optimal workspace dimension.
C
      DWORK(1) = WRKOPT
      RETURN
C *** Last line of AB01OD ***
      END