AB13DD

L-infinity norm of a state space system

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To compute the L-infinity norm of a continuous-time or
  discrete-time system, either standard or in the descriptor form,

                                  -1
     G(lambda) = C*( lambda*E - A ) *B + D .

  The norm is finite if and only if the matrix pair (A,E) has no
  eigenvalue on the boundary of the stability domain, i.e., the
  imaginary axis, or the unit circle, respectively. It is assumed
  that the matrix E is nonsingular.

Specification
      SUBROUTINE AB13DD( DICO, JOBE, EQUIL, JOBD, N, M, P, FPEAK,
     $                   A, LDA, E, LDE, B, LDB, C, LDC, D, LDD, GPEAK,
     $                   TOL, IWORK, DWORK, LDWORK, CWORK, LCWORK,
     $                   INFO )
C     .. Scalar Arguments ..
      CHARACTER          DICO, EQUIL, JOBD, JOBE
      INTEGER            INFO, LCWORK, LDA, LDB, LDC, LDD, LDE, LDWORK,
     $                   M, N, P
      DOUBLE PRECISION   TOL
C     .. Array Arguments ..
      COMPLEX*16         CWORK(  * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   D( LDD, * ), DWORK(  * ), E( LDE, * ),
     $                   FPEAK(  2 ), GPEAK(  2 )
      INTEGER            IWORK(  * )

Arguments

Mode Parameters

  DICO    CHARACTER*1
          Specifies the type of the system, as follows:
          = 'C':  continuous-time system;
          = 'D':  discrete-time system.

  JOBE    CHARACTER*1
          Specifies whether E is a general square or an identity
          matrix, as follows:
          = 'G':  E is a general square matrix;
          = 'I':  E is the identity matrix.

  EQUIL   CHARACTER*1
          Specifies whether the user wishes to preliminarily
          equilibrate the system (A,E,B,C) or (A,B,C), as follows:
          = 'S':  perform equilibration (scaling);
          = 'N':  do not perform equilibration.

  JOBD    CHARACTER*1
          Specifies whether or not a non-zero matrix D appears in
          the given state space model:
          = 'D':  D is present;
          = 'Z':  D is assumed a zero matrix.

Input/Output Parameters
  N       (input) INTEGER
          The order of the system.  N >= 0.

  M       (input) INTEGER
          The column size of the matrix B.  M >= 0.

  P       (input) INTEGER
          The row size of the matrix C.  P >= 0.

  FPEAK   (input/output) DOUBLE PRECISION array, dimension (2)
          On entry, this parameter must contain an estimate of the
          frequency where the gain of the frequency response would
          achieve its peak value. Setting FPEAK(2) = 0 indicates an
          infinite frequency. An accurate estimate could reduce the
          number of iterations of the iterative algorithm. If no
          estimate is available, set FPEAK(1) = 0, and FPEAK(2) = 1.
          FPEAK(1) >= 0, FPEAK(2) >= 0.
          On exit, if INFO = 0, this array contains the frequency
          OMEGA, where the gain of the frequency response achieves
          its peak value GPEAK, i.e.,

              || G ( j*OMEGA ) || = GPEAK ,  if DICO = 'C', or

                      j*OMEGA
              || G ( e       ) || = GPEAK ,  if DICO = 'D',

          where OMEGA = FPEAK(1), if FPEAK(2) > 0, and OMEGA is
          infinite, if FPEAK(2) = 0.

  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
          The leading N-by-N part of this array must contain the
          state dynamics matrix A.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).

  E       (input) DOUBLE PRECISION array, dimension (LDE,N)
          If JOBE = 'G', the leading N-by-N part of this array must
          contain the descriptor matrix E of the system.
          If JOBE = 'I', then E is assumed to be the identity
          matrix and is not referenced.

  LDE     INTEGER
          The leading dimension of the array E.
          LDE >= MAX(1,N), if JOBE = 'G';
          LDE >= 1,        if JOBE = 'I'.

  B       (input) DOUBLE PRECISION array, dimension (LDB,M)
          The leading N-by-M part of this array must contain the
          system input matrix B.

  LDB     INTEGER
          The leading dimension of the array B.  LDB >= max(1,N).

  C       (input) DOUBLE PRECISION array, dimension (LDC,N)
          The leading P-by-N part of this array must contain the
          system output matrix C.

  LDC     INTEGER
          The leading dimension of the array C.  LDC >= max(1,P).

  D       (input) DOUBLE PRECISION array, dimension (LDD,M)
          If JOBD = 'D', the leading P-by-M part of this array must
          contain the direct transmission matrix D.
          The array D is not referenced if JOBD = 'Z'.

  LDD     INTEGER
          The leading dimension of array D.
          LDD >= MAX(1,P), if JOBD = 'D';
          LDD >= 1,        if JOBD = 'Z'.

  GPEAK   (output) DOUBLE PRECISION array, dimension (2)
          The L-infinity norm of the system, i.e., the peak gain
          of the frequency response (as measured by the largest
          singular value in the MIMO case), coded in the same way
          as FPEAK.

Tolerances
  TOL     DOUBLE PRECISION
          Tolerance used to set the accuracy in determining the
          norm.  0 <= TOL < 1.

Workspace
  IWORK   INTEGER array, dimension (N)

  DWORK   DOUBLE PRECISION array, dimension (LDWORK)
          On exit, if INFO = 0, DWORK(1) contains the optimal value
          of LDWORK.

  LDWORK  INTEGER
          The dimension of the array DWORK.
          LDWORK >= K, where K can be computed using the following
          pseudo-code (or the Fortran code included in the routine)

             d = 6*MIN(P,M);
             c = MAX( 4*MIN(P,M) + MAX(P,M), d );
             if ( MIN(P,M) = 0 ) then
                K = 1;
             else if( N = 0 or B = 0 or C = 0 ) then
                if( JOBD = 'D' ) then
                   K = P*M + c;
                else
                   K = 1;
                end
             else
                if ( DICO = 'D' ) then
                   b = 0;  e = d;
                else
                   b = N*(N+M);  e = c;
                   if ( JOBD = Z' ) then  b = b + P*M;  end
                end
                if ( JOBD = 'D' ) then
                   r = P*M;
                   if ( JOBE = 'I', DICO = 'C',
                        N > 0, B <> 0, C <> 0 ) then
                      K = P*P + M*M;
                      r = r + N*(P+M);
                   else
                      K = 0;
                   end
                   K = K + r + c;  r = r + MIN(P,M);
                else
                   r = 0;  K = 0;
                end
                r = r + N*(N+P+M);
                if ( JOBE = 'G' ) then
                   r = r + N*N;
                   if ( EQUIL = 'S' ) then
                      K = MAX( K, r + 9*N );
                   end
                   K = MAX( K, r + 4*N + MAX( M, 2*N*N, N+b+e ) );
                else
                   K = MAX( K, r + N +
                               MAX( M, P, N*N+2*N, 3*N+b+e ) );
                end
                w = 0;
                if ( JOBE = 'I', DICO = 'C' ) then
                   w = r + 4*N*N + 11*N;
                   if ( JOBD = 'D' ) then
                      w = w + MAX(M,P) + N*(P+M);
                   end
                end
                if ( JOBE = 'E' or DICO = 'D' or JOBD = 'D' ) then
                   w = MAX( w, r + 6*N + (2*N+P+M)*(2*N+P+M) +
                            MAX( 2*(N+P+M), 8*N*N + 16*N ) );
                end
                K = MAX( 1, K, w, r + 2*N + e );
             end

          For good performance, LDWORK must generally be larger.

          An easily computable upper bound is

          K = MAX( 1, 15*N*N + P*P + M*M + (6*N+3)*(P+M) + 4*P*M +
                      N*M + 22*N + 7*MIN(P,M) ).

          The smallest workspace is obtained for DICO = 'C',
          JOBE = 'I', and JOBD = 'Z', namely

          K = MAX( 1, N*N + N*P + N*M + N +
                      MAX( N*N + N*M + P*M + 3*N + c,
                           4*N*N + 10*N ) ).

          for which an upper bound is

          K = MAX( 1, 6*N*N + N*P + 2*N*M + P*M + 11*N + MAX(P,M) +
                      6*MIN(P,M) ).

  CWORK   COMPLEX*16 array, dimension (LCWORK)
          On exit, if INFO = 0, CWORK(1) contains the optimal
          LCWORK.

  LCWORK  INTEGER
          The dimension of the array CWORK.
          LCWORK >= 1,  if N = 0, or B = 0, or C = 0;
          LCWORK >= MAX(1, (N+M)*(N+P) + 2*MIN(P,M) + MAX(P,M)),
                        otherwise.
          For good performance, LCWORK must generally be larger.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  the matrix E is (numerically) singular;
          = 2:  the (periodic) QR (or QZ) algorithm for computing
                eigenvalues did not converge;
          = 3:  the SVD algorithm for computing singular values did
                not converge;
          = 4:  the tolerance is too small and the algorithm did
                not converge.

Method
  The routine implements the method presented in [1], with
  extensions and refinements for improving numerical robustness and
  efficiency. Structure-exploiting eigenvalue computations for
  Hamiltonian matrices are used if JOBE = 'I', DICO = 'C', and the
  symmetric matrices to be implicitly inverted are not too ill-
  conditioned. Otherwise, generalized eigenvalue computations are
  used in the iterative algorithm of [1].

References
  [1] Bruinsma, N.A. and Steinbuch, M.
      A fast algorithm to compute the Hinfinity-norm of a transfer
      function matrix.
      Systems & Control Letters, vol. 14, pp. 287-293, 1990.

Numerical Aspects
  If the algorithm does not converge in MAXIT = 30 iterations
  (INFO = 4), the tolerance must be increased.

Further Comments
  If the matrix E is singular, other SLICOT Library routines
  could be used before calling AB13DD, for removing the singular
  part of the system.

Example

Program Text

*     AB13DD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, MMAX, PMAX
      PARAMETER        ( NMAX = 10, MMAX = 10, PMAX = 10 )
      INTEGER          LDA, LDB, LDC, LDD, LDE
      PARAMETER        ( LDA = NMAX, LDB = NMAX, LDC = PMAX,
     $                   LDD = PMAX, LDE = NMAX )
      INTEGER          LIWORK
      PARAMETER        ( LIWORK = NMAX )
      INTEGER          LCWORK
      PARAMETER        ( LCWORK = ( NMAX + MMAX )*( NMAX + PMAX ) +
     $                             2*MIN( PMAX, MMAX ) +
     $                             MAX( PMAX, MMAX ) )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = 15*NMAX*NMAX + PMAX*PMAX + MMAX*MMAX +
     $                            ( 6*NMAX + 3 )*( PMAX + MMAX ) +
     $                            4*PMAX*MMAX + NMAX*MMAX + 22*NMAX +
     $                            7*MIN( PMAX, MMAX ) )
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )

*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, INFO, J, M, N, P
      CHARACTER        DICO, EQUIL, JOBD, JOBE
*     .. Local Arrays ..
      INTEGER          IWORK( LIWORK )
      DOUBLE PRECISION A( LDA, NMAX ), B( LDB, MMAX ),  C( LDC, NMAX ),
     $                 D( LDD, MMAX ), DWORK( LDWORK ), E( LDE, NMAX ),
     $                 FPEAK( 2 ), GPEAK( 2 )
      COMPLEX*16       CWORK( LCWORK )
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         AB13DD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, M, P, FPEAK, TOL, DICO, JOBE, EQUIL, JOBD
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99994 ) N
      ELSE IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) M
      ELSE IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
         WRITE ( NOUT, FMT = 99992 ) P
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         IF ( LSAME( JOBE, 'G' ) )
     $      READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,N )
         READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
         READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
         IF ( LSAME( JOBD, 'D' ) )
     $      READ ( NIN, FMT = * ) ( ( D(I,J), J = 1,M ), I = 1,P )
*        Computing the Linf norm.
         CALL AB13DD( DICO, JOBE, EQUIL, JOBD, N, M, P, FPEAK, A, LDA,
     $                E, LDE, B, LDB, C, LDC, D, LDD, GPEAK, TOL, IWORK,
     $                DWORK, LDWORK, CWORK, LCWORK, INFO )
*
         IF ( INFO.EQ.0 ) THEN
            IF ( GPEAK( 2 ).EQ.ZERO ) THEN
               WRITE ( NOUT, FMT = 99991 )
            ELSE
               WRITE ( NOUT, FMT = 99997 )
               WRITE ( NOUT, FMT = 99995 ) GPEAK( 1 )
            END IF
            IF ( FPEAK( 2 ).EQ.ZERO ) THEN
               WRITE ( NOUT, FMT = 99990 )
            ELSE
               WRITE ( NOUT, FMT = 99996 )
               WRITE ( NOUT, FMT = 99995 ) FPEAK( 1 )
            END IF
         ELSE
            WRITE( NOUT, FMT = 99998 ) INFO
         END IF
      END IF
      STOP
*
99999 FORMAT (' AB13DD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (/' INFO on exit from AB13DD =',I2)
99997 FORMAT (/' The L_infty norm of the system is'/)
99996 FORMAT (/' The peak frequency is'/)
99995 FORMAT (D17.10)
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' M is out of range.',/' M = ',I5)
99992 FORMAT (/' P is out of range.',/' P = ',I5)
99991 FORMAT (/' The L_infty norm of the system is infinite')
99990 FORMAT (/' The peak frequency is infinite'/)
      END
Program Data
 AB13CD EXAMPLE PROGRAM DATA
   6     1     1     0.0     1.0   0.000000001     C     I     N     D   
   0.0  1.0     0.0   0.0      0.0  0.0
  -0.5 -0.0002  0.0   0.0      0.0  0.0
   0.0  0.0     0.0   1.0      0.0  0.0
   0.0  0.0    -1.0  -0.00002  0.0  0.0
   0.0  0.0     0.0   0.0      0.0  1.0
   0.0  0.0     0.0   0.0     -2.0 -0.000002
   1.0 
   0.0  
   1.0 
   0.0 
   1.0 
   0.0
   1.0  0.0  1.0  0.0  1.0  0.0
   0.0

Program Results
 AB13DD EXAMPLE PROGRAM RESULTS


 The L_infty norm of the system is

 0.5000000001D+06

 The peak frequency is

 0.1414213562D+01

Return to index