        SUBROUTINE AOLSYS(SYSFIL,ANLFIL,ERRCD,ERRMSG)
C
C       FUNCTION:
CF
CF
C       USAGE:
CU
CU
C       INPUTS:
CI
CI
C       OUTPUTS:
CO
CO      ERROR_CODE      integer
CO              value                   range           description
CO              INSYS code              [0..100]        see INSYS
CO              100 + rg code           [101..200]      see rg (eispack)
CO              200 - mvzero code       [201..300]      see mvzero
CO              300 + ahcon code        [301..400]      see ahcon (CC calc)
CO              400 + ahcon code        [401..500]      see ahcon (CO calc)
CO
C       ALGORITHM:
CA
CA
C       MACHINE DEPENDENCIES:
CM
CM
C       HISTORY:
CH
CH      written by:             J. Douglas Birdwell
CH      date:                   23-july-1986
CH      current version:        2.1
CH      modifications:          standardized,modified to accomodate
CH                              header.c bb: 8-86.
CH                              fixed problem with z2 and z1
CH                              bb and jdb :8-86.
CH				added dpcom: 7/16/88 jdb
CH
C       ROUTINES CALLED:
CC
CC
C       COMMON MEMORY USED:
CM
CM      DPCOM -- see dpcommon.f and dpcom.f
CM
C----------------------------------------------------------------------
C       written for:    The CASCADE Project
C                       Oak Ridge National Laboratory
C                       U.S. Department of Energy
C                       contract number DE-AC05-840R21400
C                       subcontract number ##########
C                       organization:   xxxxxx
C----------------------------------------------------------------------
C       THIS SOFTWARE IS IN THE PUBLIC DOMAIN
C       NO RESTRICTIONS ON ITS USE ARE IMPLIED
C----------------------------------------------------------------------
C
        INCLUDE         'Parameter.f'
C
        DOUBLE PRECISION        A(SIZE,SIZE)
        DOUBLE PRECISION        B(SIZE,SIZE)
        DOUBLE PRECISION        C(SIZE,SIZE)
        DOUBLE PRECISION        D(SIZE,SIZE)
        DOUBLE PRECISION        ATEMP(SIZE,SIZE)
        DOUBLE PRECISION        Z1(SIZE2,SIZE2)
        DOUBLE PRECISION        Z2(SIZE2,SIZE2)
        DOUBLE PRECISION        REZERO(SIZE)
        DOUBLE PRECISION        IMZERO(SIZE)
        DOUBLE PRECISION        REPOLE(SIZE)
        DOUBLE PRECISION        IMPOLE(SIZE)
        DOUBLE PRECISION        POLES(SIZE,2)
        DOUBLE PRECISION        ZEROS(SIZE,2)
        DOUBLE PRECISION        BETA(SIZE)
        DOUBLE PRECISION        SUM(SIZE)
        DOUBLE PRECISION        DUMMY(SIZE)
        DOUBLE PRECISION        SIGMA(SIZE)
        DOUBLE PRECISION        U
        DOUBLE PRECISION        V
        DOUBLE PRECISION        EIGVEC(SIZE,SIZE)
        DOUBLE PRECISION        EPS
        DOUBLE PRECISION        EPS1
        DOUBLE PRECISION        EPSLON
        DOUBLE PRECISION        WNZERO
        DOUBLE PRECISION        THZERO
        DOUBLE PRECISION        FREQ
        DOUBLE PRECISION        WNMIN
        DOUBLE PRECISION        WNMAX
C
        INTEGER                 NZEROS
        INTEGER                 NRANK
        INTEGER                 NSTATS
        INTEGER                 NINPS
        INTEGER                 NOUTS
        INTEGER                 ERRCD
        INTEGER                 IPVT(SIZE)
        INTEGER                 JPVT(SIZE)
        INTEGER                 NPORHP
        INTEGER                 NPJWAX
        INTEGER                 NPORIG
        INTEGER                 NPOLHP
        INTEGER                 NZORHP
        INTEGER                 NZJWAX
        INTEGER                 NZORIG
        INTEGER                 NZOLHP
        INTEGER                 NOBSMD
        INTEGER                 NCONMD
        INTEGER                 IDXSIG(SIZE)
        INTEGER                 NCON0(SIZE)
        INTEGER                 NOBS0(SIZE)
        INTEGER                 IDXPOR(SIZE)
        INTEGER                 ISEED
        INTEGER                 LENGTH
        INTEGER                 MERR
        INTEGER                 MAXLMN
C
        LOGICAL                 INSBSP
        LOGICAL                 CON(SIZE)
C
        CHARACTER*(*)           ANLFIL
        CHARACTER*(*)           SYSFIL
        CHARACTER*(*)           ERRMSG
        CHARACTER*30            STRING
C
	INCLUDE			'dpcom.f'
C
        ERRCD=0
C
C       READ IN THE SYSTEM FROM THE SYSTEM CONTAINER FILE
C
        CALL INSYS(SYSFIL,NINPS,NOUTS,NSTATS,
     2     SIZE,SIZE,SIZE,A,B,C,D,ERRCD)
          CLOSE (UNIT=UNIT1)
          IF (ERRCD .NE. 0) THEN
            ERRMSG = 'AOLSYS: Fatal error from INSYS '//
     1      'accessing '//SYSFIL
            RETURN
          END IF
C
C--calculate the poles of the system
C
        CALL SAVE (SIZE,SIZE,NSTATS,NSTATS,A,ATEMP)
        CALL RG (SIZE,NSTATS,ATEMP,REPOLE,IMPOLE,
     .          0,U,IPVT,DUMMY,ERRCD)
C
        IF (ERRCD .NE. 0) THEN
          ERRMSG='AOLSYS: Fatal error from '//
     2          'eispack (RG) ERRCD = 100 + return code: '
          ERRCD = 100 + ERRCD
          RETURN
        END IF
C
C--calculate multivariable zeros of system
C
        CALL MVZERO (NOUTS,NINPS,NSTATS,
     2          SIZE,SIZE,SIZE2,
     3          A,B,C,D,Z1,Z2,NZEROS,NRANK,
     4          REZERO,IMZERO,
     5          BETA,SUM,DUMMY)
C
        IF (NZEROS .LT. 0) THEN
          ERRMSG = ' AOLSYS: Zero calculation failed.'
          ERRCD = 200 - NZEROS
          RETURN
        END IF
C
C--find maximum natural frequency
C
          WNMAX = 0.D0
C
          DO 10, I = 1, NSTATS
C
            FREQ = SQRT(REPOLE(I)**2
     2          + IMPOLE(I)**2)
            IF (FREQ .GT. WNMAX) WNMAX = FREQ
C
 10       CONTINUE
C
C--compute machine precision
C
        EPS = 1.D0
20      CONTINUE
        EPS = EPS / 2.D0
        EPS1 = EPS + 1.D0
        IF (EPS1 .GT. 1.D0) GO TO 20
        EPS = 2.D0 * EPS
C
C--set width of "jw" axis
C
        EPSLON = EPS**0.3333333333333333D0 * MAX(WNMAX,1.D0)
C
C--find minimum natural frequency
C
          WNMIN = 1.D30
C
          DO 30, I = 1, NSTATS
C
            FREQ = SQRT(REPOLE(I)**2
     2          + IMPOLE(I)**2)
            IF (FREQ .LT. WNMIN .AND. FREQ .GT. EPSLON) WNMIN = FREQ
C
 30       CONTINUE
C
          IF (WNMIN .EQ. 1.D30) WNMIN = 0.001D0
C
C--count zeros in ORHP, on jw axis, at 0, and in OLHP.
C
        NZORHP = 0
        NZJWAX = 0
        NZORIG = 0
        NZOLHP = 0
C
        DO 40, I = 1, NZEROS
C
C       Save the zeros of the system
C
          ZEROS(I,1) = REZERO(I)
          ZEROS(I,2) = IMZERO(I)
C
          IF (REZERO(I) .LT. -EPSLON) THEN
            NZOLHP = NZOLHP + 1
          ELSE IF (REZERO(I) .GT. EPSLON) THEN
            NZORHP = NZORHP + 1
          ELSE IF (DABS(IMZERO(I)) .GT. EPSLON) THEN
            NZJWAX = NZJWAX + 1
          ELSE
            NZORIG = NZORIG + 1
          END IF
C
 40     CONTINUE
C
C--count poles in ORHP, on jw axis, at 0, and in OLHP.
C
        NPORHP = 0
        NPJWAX = 0
        NPORIG = 0
        NPOLHP = 0
C
        DO 50, I = 1, NSTATS
C
C    Save the poles of the system.
C
          POLES(I,1) = REPOLE(I)
          POLES(I,2) = IMPOLE(I)
C
          IF (REPOLE(I) .LT. -EPSLON) THEN
            NPOLHP = NPOLHP + 1
          ELSE IF (REPOLE(I) .GT. EPSLON) THEN
            NPORHP = NPORHP + 1
          ELSE IF (DABS(IMPOLE(I)) .GT. EPSLON) THEN
            NPJWAX = NPJWAX + 1
          ELSE
            NPORIG = NPORIG + 1
C
C--keep track of where the eigenvectors corresponding to zero eigenvalues
C--are for controllability calculations
C
            IDXPOR(NPORIG) = I
C
          END IF
C
 50       CONTINUE
C
C--create output file
C
        OPEN (UNIT=UNIT1,FILE=ANLFIL)
C
C--write pole data
C
          WRITE (UNIT1,*)  NPORHP,
     2                     NPJWAX,
     3                     NPORIG,
     4                     NPOLHP
C
C--write zero data
C
          WRITE (UNIT1,*)  NZORHP,
     2                     NZJWAX,
     3                     NZORIG,
     4                     NZOLHP
C
C--find natural frequency of smallest non-minimum phase zero
C
          IF (NZORHP .EQ. 0) THEN
            WNZERO = 0.D0
          ELSE
            WNZERO = 1.0D30
            DO 60, I = 1, NZEROS
              IF (REZERO(I) .GT. EPSLON) THEN
                THZERO = SQRT(REZERO(I)**2
     2                  + IMZERO(I)**2)
                WNZERO = MIN(WNZERO,THZERO)
              END IF
 60         CONTINUE
          END IF
C
          WRITE (UNIT1,*) WNZERO
C
C--begin controllability and observability calculations
C
        ISEED = 0
        CALL AHCON (SIZE,NSTATS,NINPS,A,B,REPOLE,
     .          IMPOLE,REZERO,IMZERO,
     .          DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD)
        IF (ERRCD .NE. 0) THEN
          ERRCD = 300 + ERRCD
          ERRMSG = ' AOLSYS: Controllability calculation failed.'//
     1             ' System is not controllable.'
          RETURN
        END IF
        NCONMD = 0
        DO 70, I = 1, NSTATS
          IF (CON(I)) NCONMD = NCONMD + 1
 70     CONTINUE
C
        CALL TRNATB (SIZE,SIZE2,NSTATS,NSTATS,A,Z1)
        CALL TRNATB (SIZE,SIZE2,NOUTS,NSTATS,C,Z2)
        CALL AHCON (SIZE2,NSTATS,NOUTS,Z1,Z2,REPOLE,
     .          IMPOLE,REZERO,IMZERO,
     .          DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD)
        IF (ERRCD .NE. 0) THEN
          ERRCD = 400 + ERRCD
          ERRMSG = ' AOLSYS: Observability calculation failed.'//
     1             ' System is not observable.'
          RETURN
        END IF
        NOBSMD = 0
        DO 80, I = 1, NSTATS
          IF (CON(I)) NOBSMD = NOBSMD + 1
 80     CONTINUE
C
C--now compute the number of controllable s=0 modes from each input
C
C--if there are no poles at s=0, skip the work
C
        IF (NPORIG .GT. 0) THEN
C
C--for each input, check controllability using ahcon
C
          DO 90, I = 1, NINPS
C
            CALL AHCON (SIZE,NSTATS,1,A,B(1,I),REPOLE,
     .          IMPOLE,REZERO,IMZERO,
     .          DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD)
            IF (ERRCD .NE. 0) THEN
              CALL ICNVRT(0,I,STRING,LENGTH,MERR)
              ERRCD = 500 + ERRCD
              ERRMSG=' AOLSYS: Controllability calculation '//
     .               'failed for input '//STRING(1:LENGTH)//'.'
              RETURN
            END IF
C
C--check for each zero eigenvalue
C
            NCON0(I) = 0
            DO 100, J = 1, NSTATS
              FREQ = SQRT(REPOLE(J)**2
     .          + IMPOLE(J)**2)
C
              IF (FREQ .LT. EPSLON) THEN
C
C--this is a zero eigenvalue
C--now, look for the eigenvalue closest to this eigenvalue
C
                DO 110, K = 1, NSTATS
                  IF (IPVT(K) .EQ. J) INDEX = K
 110            CONTINUE
C
C--and check to see if it is sufficiently far away from the 0 eigenvalue
C
                IF (CON(INDEX)) NCON0(I) = NCON0(I) + 1
              END IF
 100        CONTINUE
 90       CONTINUE
C
C--else, there are no s=0 modes and all this can be skipped
C
        ELSE
C
          DO 120, I = 1, NINPS
            NCON0(I) = 0
 120      CONTINUE
C
        END IF
C
C--now compute the number of observable s=0 modes from each output
C
C--if there are no poles at s=0, skip the work
C
        IF (NPORIG .GT. 0) THEN
C
C--remember that z1 contains a-transpose and z2 contains b-transpose
C
C--for each output, check observability using ahcon
C
          DO 130 I = 1, NOUTS
            CALL AHCON (SIZE2,NSTATS,1,Z1,Z2(1,I),REPOLE,
     .          IMPOLE,REZERO,IMZERO,
     .          DUMMY,SIGMA,IPVT,JPVT,CON,ATEMP,ISEED,ERRCD)
            IF (ERRCD .NE. 0) THEN
              CALL ICNVRT(0,I,STRING,LENGTH,MERR)
              ERRCD = 600 + ERRCD
              ERRMSG = ' AOLSYS: observability calculation '//
     .          'failed for output '//STRING(1:LENGTH)//'.'
              RETURN
            END IF
C
C--check for each zero eigenvalue
C
            NOBS0(I) = 0
            DO 140, J = 1, NSTATS
              FREQ = SQRT(REPOLE(J)**2
     .          + IMPOLE(J)**2)
C
              IF (FREQ .LT. EPSLON) THEN
C
C--this is a zero eigenvalue
C--now, look for the eigenvalue closest to this eigenvalue
C
                DO 150, K = 1, NSTATS
                  IF (IPVT(K) .EQ. J) INDEX = K
 150            CONTINUE
C
C--and check to see if it is sufficiently far away from the 0 eigenvalue
C
                IF (CON(INDEX)) NOBS0(I) = NOBS0(I) + 1
              END IF
 140        CONTINUE
 130      CONTINUE
C
C--else, there are no s=0 modes and all this can be skipped
C
        ELSE
C
          DO 160, I = 1, NOUTS
            NOBS0(I) = 0
 160      CONTINUE
C
        END IF
C
C--print out number of controllable modes of the system
C
          WRITE (UNIT1,*) NCONMD
C
C--print out number of observable modes of the system
C
          WRITE (UNIT1,*) NOBSMD
C
C--PRINT OUT VECTOR OF NUMBER OF CONTROLLABLE S=0 MODES OF
C--the system from input i, for i=1,..,NINPS
C
          DO 170, I = 1, NINPS
            WRITE (UNIT1,*) NCON0(I)
 170      CONTINUE
C
C--print out vector of number of observable s=0 modes of
C--the system from output i, for i =1,..,NOUTS
C
          DO 180, I = 1, NOUTS
            WRITE (UNIT1,*) NOBS0(I)
 180      CONTINUE
C
          WRITE (UNIT1,*) WNMIN, WNMAX
C
C Write out poles of system
C
          WRITE(UNIT1,185)
 185      FORMAT(1X,'POLES:')
C
          DO 190, I=1,NSTATS
             WRITE(UNIT1,*) POLES(I,1),POLES(I,2)
 190      CONTINUE
C
C  Write out zeros of system
C
          WRITE(UNIT1,195)
 195      FORMAT(1X,'ZEROS:')
C
          DO 200, I=1,NZEROS
             WRITE(UNIT1,*) ZEROS(I,1),ZEROS(I,2)
 200      CONTINUE
C
C--end of processing
C
        CLOSE (UNIT=UNIT1)
        RETURN
        END
