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

      SUBROUTINE IB01PX( JOB, NOBR, N, M, L, UF, LDUF, UN, LDUN, UL,
     $                   LDUL, PGAL, LDPGAL, K, LDK, R, LDR, X, B, LDB,
     $                   D, LDD, TOL, IWORK, DWORK, LDWORK, IWARN,
     $                   INFO )
C
C     RELEASE 4.0, WGS COPYRIGHT 2000.
C
C     PURPOSE
C
C     To build and solve the least squares problem  T*X = Kv,  and 
C     estimate the matrices B and D of a linear time-invariant (LTI)
C     state space model, using the solution  X,  and the singular
C     value decomposition information and other intermediate results,
C     provided by other routines. 
C
C     The matrix  T  is computed as a sum of Kronecker products,
C
C        T = T + kron(Uf(:,(i-1)*m+1:i*m),N_i),  for i = 1 : s,
C
C     (with  T  initialized by zero), where  Uf  is the triangular
C     factor of the QR factorization of the future input part (see 
C     SLICOT Library routine IB01ND),  N_i  is given by the i-th block
C     row of the matrix
C
C            [ Q_11  Q_12  ...  Q_1,s-2  Q_1,s-1  Q_1s ]   [ I_L  0  ] 
C            [ Q_12  Q_13  ...  Q_1,s-1    Q_1s    0   ]   [         ]
C        N = [ Q_13  Q_14  ...    Q_1s      0      0   ] * [         ],
C            [  :     :            :        :      :   ]   [         ]
C            [ Q_1s   0    ...     0        0      0   ]   [  0  GaL ]
C
C     and where
C
C               [   -L_1|1    ]          [ M_i-1 - L_1|i ]
C        Q_11 = [             ],  Q_1i = [               ],  i = 2:s,
C               [ I_L - L_2|1 ]          [     -L_2|i    ]
C
C     are  (n+L)-by-L  matrices, and  GaL  is built from the first  n  
C     relevant singular vectors,  GaL = Un(1:L(s-1),1:n),  computed 
C     by IB01ND.
C
C     The vector  Kv  is vec(K), with the matrix  K  defined by
C
C        K = [ K_1  K_2  K_3  ...  K_s ],
C
C     where  K_i = K(:,(i-1)*m+1:i*m),  i = 1:s,  is  (n+L)-by-m.
C     The given matrices are  Uf,  GaL,  and  
C
C            [ L_1|1  ...  L_1|s ]
C        L = [                   ],   (n+L)-by-L*s,
C            [ L_2|1  ...  L_2|s ]
C
C        M = [ M_1  ...  M_s-1 ],      n-by-L*(s-1),  and
C        K,                            (n+L)-by-m*s.
C
C     Matrix  M  is the pseudoinverse of the matrix  GaL,  computed by
C     SLICOT Library routine IB01PD.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             Specifies which of the matrices B and D should be
C             computed, as follows:
C             = 'B':  compute the matrix B, but not the matrix D;
C             = 'D':  compute both matrices B and D.
C
C     Input/Output Parameters
C
C     NOBR    (input) INTEGER
C             The number of block rows,  s,  in the input and output
C             Hankel matrices processed by other routines.  NOBR > 1.
C
C     N       (input) INTEGER
C             The order of the system.  NOBR > N > 0.
C
C     M       (input) INTEGER
C             The number of system inputs.  M >= 0.
C
C     L       (input) INTEGER
C             The number of system outputs.  L > 0.
C
C     UF      (input/output) DOUBLE PRECISION array, dimension 
C             ( LDUF,M*NOBR ) 
C             On entry, the leading  M*NOBR-by-M*NOBR  upper triangular
C             part of this array must contain the upper triangular 
C             factor of the QR factorization of the future input part,
C             as computed by SLICOT Library routine IB01ND.
C             The strict lower triangle need not be set to zero. 
C             On exit, the leading  M*NOBR-by-M*NOBR  upper triangular
C             part of this array is unchanged, and the strict lower 
C             triangle is set to zero.
C
C     LDUF    INTEGER
C             The leading dimension of the array  UF.
C             LDUF >= MAX( 1, M*NOBR ).
C
C     UN      (input) DOUBLE PRECISION array, dimension ( LDUN,N ) 
C             The leading  L*(NOBR-1)-by-N  part of this array must
C             contain the matrix  GaL,  i.e., the leading part of the
C             first  N  columns of the matrix  Un  of relevant singular
C             vectors.
C
C     LDUN    INTEGER
C             The leading dimension of the array  UN.
C             LDUN >= L*(NOBR-1).
C
C     UL      (input/output) DOUBLE PRECISION array, dimension 
C             ( LDUL,L*NOBR )
C             On entry, the leading  (N+L)-by-L*NOBR  part of this array
C             must contain the given matrix  L.
C             On exit, if  M > 0,  the leading  (N+L)-by-L*NOBR  part of
C             this array is overwritten by the matrix  
C             [ Q_11  Q_12  ...  Q_1,s-2  Q_1,s-1  Q_1s ]. 
C
C     LDUL    INTEGER
C             The leading dimension of the array  UL.  LDUL >= N+L.
C
C     PGAL    (input) DOUBLE PRECISION array, dimension 
C             ( LDPGAL,L*(NOBR-1) )
C             The leading  N-by-L*(NOBR-1)  part of this array must
C             contain the pseudoinverse of the matrix  GaL,  computed by
C             SLICOT Library routine IB01PD.
C
C     LDPGAL  INTEGER
C             The leading dimension of the array  PGAL.  LDPGAL >= N.
C
C     K       (input) DOUBLE PRECISION array, dimension ( LDK,M*NOBR )
C             The leading  (N+L)-by-M*NOBR  part of this array must
C             contain the given matrix  K.
C            
C     LDK     INTEGER
C             The leading dimension of the array  K.  LDK >= N+L.
C
C     R       (output) DOUBLE PRECISION array, dimension ( LDR,M*(N+L) )
C             The leading  (N+L)*M*NOBR-by-M*(N+L)  part of this array
C             contains details of the complete orthogonal factorization
C             of the coefficient matrix  T  of the least squares problem
C             which is solved for getting the system matrices B and D.
C            
C     LDR     INTEGER
C             The leading dimension of the array  R.
C             LDR >= MAX( 1, (N+L)*M*NOBR ).
C            
C     X       (output) DOUBLE PRECISION array, dimension 
C             ( (N+L)*M*NOBR ) 
C             The leading  M*(N+L)  elements of this array contain the
C             least squares solution of the system  T*X = Kv.
C             The remaining elements are used as workspace (to store the
C             corresponding part of the vector Kv = vec(K)).
C            
C     B       (output) DOUBLE PRECISION array, dimension ( LDB,M )
C             The leading N-by-M part of this array contains the system
C             input matrix.  
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= N.
C
C     D       (output) DOUBLE PRECISION array, dimension ( LDD,M )
C             If  JOB = 'D',  the leading L-by-M part of this array 
C             contains the system input-output matrix.  
C             If  JOB = 'B',  this array is not referenced.  
C
C     LDD     INTEGER
C             The leading dimension of the array D.
C             LDD >= L, if  JOB = 'D';
C             LDD >= 1, if  JOB = 'B'.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used for estimating the rank of
C             matrices. If the user sets  TOL > 0,  then the given value
C             of  TOL  is used as a lower bound for the reciprocal
C             condition number;  an m-by-n matrix whose estimated
C             condition number is less than  1/TOL  is considered to  
C             be of full rank.  If the user sets  TOL <= 0,  then an 
C             implicitly computed, default tolerance, defined by 
C             TOLDEF = m*n*EPS,  is used instead, where  EPS  is the 
C             relative machine precision (see LAPACK Library routine 
C             DLAMCH).
C
C     Workspace
C
C     IWORK   INTEGER array, dimension ( M*(N+L) )
C
C     DWORK   DOUBLE PRECISION array, dimension ( LDWORK )
C             On exit, if  INFO = 0,  DWORK(1) returns the optimal value
C             of LDWORK,  and, if  M > 0,  DWORK(2)  contains the  
C             reciprocal condition number of the triangular factor of
C             the matrix  T.
C             On exit, if  INFO = -26,  DWORK(1)  returns the minimum
C             value of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= MAX( (N+L)*(N+L), 4*M*(N+L)+1 ).
C             For good performance,  LDWORK  should be larger.
C
C     Warning Indicator
C
C     IWARN   INTEGER
C             = 0:  no warning;
C             = 4:  the least squares problem to be solved has a
C                   rank-deficient coefficient matrix. 
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     The matrix  T  is computed, evaluating the sum of Kronecker
C     products, and then the linear system  T*X = Kv  is solved in a
C     least squares sense. The matrices  B  and  D  are then directly
C     obtained from the least squares solution.
C
C     REFERENCES
C
C     [1] Verhaegen M., and Dewilde, P.
C         Subspace Model Identification. Part 1: The output-error 
C         state-space model identification class of algorithms.
C         Int. J. Control, 56, pp. 1187-1210, 1992.
C
C     [2] Van Overschee, P., and De Moor, B.
C         N4SID: Two Subspace Algorithms for the Identification
C         of Combined Deterministic-Stochastic Systems.
C         Automatica, Vol.30, No.1, pp. 75-93, 1994.
C
C     [3] Van Overschee, P.
C         Subspace Identification : Theory - Implementation -
C         Applications.
C         Ph. D. Thesis, Department of Electrical Engineering,
C         Katholieke Universiteit Leuven, Belgium, Feb. 1995.
C
C     NUMERICAL ASPECTS
C
C     The implemented method is numerically stable.
C
C     CONTRIBUTOR
C
C     V. Sima, Katholieke Universiteit Leuven, Feb. 2000.
C
C     REVISIONS
C
C     -
C
C     KEYWORDS
C
C     Identification methods; least squares solutions; multivariable
C     systems; QR decomposition; singular value decomposition.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO  = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      DOUBLE PRECISION   TOL
      INTEGER            INFO, IWARN, L, LDB, LDD, LDK, LDPGAL, LDR,
     $                   LDUF, LDUL, LDUN, LDWORK, M, N, NOBR
      CHARACTER          JOB
C     .. Array Arguments ..
      DOUBLE PRECISION   B(LDB, *), D(LDD, *), DWORK(*), K(LDK, *),  
     $                   PGAL(LDPGAL, *), R(LDR, *), UF(LDUF, *), 
     $                   UL(LDUL, *), UN(LDUN, *), X(*)
      INTEGER            IWORK( * )
C     .. Local Scalars ..
      DOUBLE PRECISION   RCOND, TOLL
      INTEGER            I, IERR, J, JWORK, LDUN2, LNOBR, LP1, MAXWRK,
     $                   MINWRK, MKRON, MNOBR, NKRON, NP1, NPL, RANK
      LOGICAL            WITHB, WITHD
C     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH, LSAME
C     .. External Subroutines ..
      EXTERNAL           DGELSY, DGEMM, DLACPY, DLASET, DTRCON, MB01VD,
     $                   XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      WITHD = LSAME( JOB, 'D' )
      WITHB = LSAME( JOB, 'B' ) .OR. WITHD
      MNOBR = M*NOBR
      LNOBR = L*NOBR
      LDUN2 = LNOBR - L
      LP1   = L + 1
      NP1   = N + 1
      NPL   = N + L
      IWARN = 0
      INFO  = 0
C
C     Check the scalar input parameters.
C
      IF( .NOT.WITHB ) THEN
         INFO = -1
      ELSE IF( NOBR.LE.1 ) THEN
         INFO = -2
      ELSE IF( N.GE.NOBR .OR. N.LE.0 ) THEN
         INFO = -3
      ELSE IF( M.LT.0 ) THEN
         INFO = -4
      ELSE IF( L.LE.0 ) THEN
         INFO = -5
      ELSE IF( LDUF.LT.MAX( 1, MNOBR ) ) THEN
         INFO = -7
      ELSE IF( LDUN.LT.LDUN2 ) THEN
         INFO = -9
      ELSE IF( LDUL.LT.NPL ) THEN
         INFO = -11
      ELSE IF( LDPGAL.LT.N ) THEN
         INFO = -13
      ELSE IF( LDK.LT.NPL ) THEN
         INFO = -15
      ELSE IF( LDR.LT.MAX( 1, MNOBR*NPL ) ) THEN
         INFO = -17
      ELSE IF( LDB.LT.N ) THEN
         INFO = -20
      ELSE IF( LDD.LT.1 .OR. ( WITHD .AND. LDD.LT.L ) ) THEN
         INFO = -22
      ELSE IF( LDWORK.GE.1 ) THEN
C
C        Compute workspace.
C        (Note: Comments in the code beginning "Workspace:" describe the
C         minimal amount of workspace needed at that point in the code,
C         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
         MINWRK = MAX( NPL*NPL, 4*M*NPL + 1 )
C
         IF ( LDWORK.LT.MINWRK ) THEN
            INFO = -26
            DWORK( 1 ) = MINWRK
         END IF
      END IF
C
C     Return if there are illegal arguments.
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'IB01PX', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF( M.EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Construct the matrix  [ Q_11  Q_12  ...  Q_1,s-1  Q_1s ]  in  UL.
C
      DO 20 J = 1, L
C
         DO 10 I = 1, NPL
            UL(I,J) = -UL(I,J)
   10    CONTINUE
C
         UL(N+J,J) = ONE + UL(N+J,J)
   20 CONTINUE
C
      DO 50 J = LP1, LNOBR
C
         DO 30 I = 1, N
            UL(I,J) = PGAL(I,J-L) - UL(I,J)
   30    CONTINUE
C
         DO 40 I = NP1, NPL
            UL(I,J) = -UL(I,J)
   40    CONTINUE
C
   50 CONTINUE
C
C     Compute the coefficient matrix T using Kronecker products.
C     Workspace: (N+L)*(N+L).
C     In the same loop, vectorize K in X.
C
      CALL DLASET( 'Full', MNOBR*NPL, M*NPL, ZERO, ZERO, R, LDR )
      CALL DLASET( 'Lower', MNOBR-1, MNOBR-1, ZERO, ZERO, UF(2,1), 
     $             LDUF )
      JWORK = NPL*L + 1
C
      DO 60 I = 1, NOBR
         CALL DLACPY( 'Full', NPL, L, UL(1,(I-1)*L+1), LDUL, DWORK,
     $                NPL )
         IF ( I.LT.NOBR ) THEN
            CALL DGEMM ( 'NoTranspose', 'NoTranspose', NPL, N,
     $                   L*(NOBR-I), ONE, UL(1,I*L+1), LDUL, UN, LDUN,
     $                   ZERO, DWORK(JWORK), NPL )
         ELSE
            CALL DLASET( 'Full', NPL, N, ZERO, ZERO, DWORK(JWORK), NPL )
         END IF
         CALL MB01VD( 'NoTranspose', 'NoTranspose', MNOBR, M, NPL, 
     $                NPL, ONE, ONE, UF(1,(I-1)*M+1), LDUF, DWORK, 
     $                NPL, R, LDR, MKRON, NKRON, IERR )
         CALL DLACPY( 'Full', NPL, M, K(1,(I-1)*M+1), LDK, 
     $                X((I-1)*NKRON+1), NPL )
   60 CONTINUE
C
C     Compute the tolerance.
C
      TOLL = TOL
      IF( TOLL.LE.ZERO ) 
     $   TOLL = MKRON*NKRON*DLAMCH( 'Precision' )
C
C     Solve the least square problem  T*X = vec(K).
C     Workspace:  need   4*M*(N+L)+1;
C                 prefer 3*M*(N+L)+(M*(N+L)+1)*NB. 
C
      DO 70 I = 1, NKRON
         IWORK(I) = 0
   70 CONTINUE
C
      CALL DGELSY( MKRON, NKRON, 1, R, LDR, X, MKRON, IWORK, TOLL, RANK,
     $             DWORK, LDWORK, IERR )
      MAXWRK = DWORK(1)
C
C     Compute the reciprocal of the condition number of the triangular 
C     factor  R  of  T.
C     Workspace: need 3*M*(N+L).
C
      CALL DTRCON( '1-norm', 'Upper', 'NonUnit', NKRON, R, LDR, RCOND,
     $             DWORK, IWORK, IERR )
C
      IF ( RANK.LT.NKRON ) THEN
C     
C        The least squares problem is rank-deficient.
C     
         IWARN = 4
      END IF
C
C     Construct the matrix  D,  if needed.
C
      IF ( WITHD ) 
     $   CALL DLACPY( 'Full', L, M, X, NPL, D, LDD )
C
C     Construct the matrix  B.
C
      CALL DLACPY( 'Full', N, M, X(LP1), NPL, B, LDB )
C
C     Return optimal workspace in DWORK(1) and reciprocal condition
C     number in  DWORK(2).
C
      DWORK(1) = MAX( MINWRK, MAXWRK )
      DWORK(2) = RCOND
C
      RETURN
C
C *** Last line of IB01PX ***
      END