SUBROUTINE  LQR(SYSFIL,BALFIL,OUTFIL,RHO,ERRCD,ERRMSG)
C
C       FUNCTION:
CF
CF
C       USAGE:
CU
CU
C       INPUTS:
CI
CI
C       OUTPUTS:
CO
CO
C       ALGORITHM:
CA
CA
C       MACHINE DEPENDENCIES:
CM
CM
C       HISTORY:
CH
CH      written by:  Bobby Bodenheimer
CH      date:        12-86
CH      current version:  1.2
CH      modifications:    declared undeclared variables - jdb - 5/11/88
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 37B-07685C S13
C                       organization:   The University of Tennessee
C----------------------------------------------------------------------
C       THIS SOFTWARE IS IN THE PUBLIC DOMAIN
C       NO RESTRICTIONS ON ITS USE ARE IMPLIED
C----------------------------------------------------------------------
C
C
      INCLUDE               'Parameter.f'
C
      DOUBLE PRECISION      A(SIZE2,SIZE2)
      DOUBLE PRECISION      ACL(SIZE2,SIZE2)
      DOUBLE PRECISION      B(SIZE2,SIZE2)
      DOUBLE PRECISION      C(SIZE2,SIZE2)
      DOUBLE PRECISION      D(SIZE2,SIZE2)
      DOUBLE PRECISION      CQC(SIZE2,SIZE2)
      DOUBLE PRECISION      E(SIZE2,SIZE2)
      DOUBLE PRECISION      F(SIZE2,SIZE2)
      DOUBLE PRECISION      G(SIZE2,SIZE2)
      DOUBLE PRECISION      Q(SIZE2,SIZE2)
      DOUBLE PRECISION      R(SIZE2,SIZE2)
      DOUBLE PRECISION      RI(SIZE2,SIZE2)
      DOUBLE PRECISION      RS(SIZE2,SIZE2)
      DOUBLE PRECISION      S(SIZE2,SIZE2)
      DOUBLE PRECISION      U(SIZE2,SIZE2)
      DOUBLE PRECISION      WK(SIZE2,SIZE2)
      DOUBLE PRECISION      Z(SIZE2,SIZE2)
C
      DOUBLE COMPLEX        P(SIZE2,SIZE2)
C
      DOUBLE PRECISION      ALFI(SIZE2)
      DOUBLE PRECISION      ALFR(SIZE2)
      DOUBLE PRECISION      BETA(SIZE2)
      DOUBLE PRECISION      CPERM(SIZE2)
      DOUBLE PRECISION      CSCALE(SIZE2)
C
      DOUBLE PRECISION      RHO
      DOUBLE PRECISION      RTOL
      DOUBLE PRECISION      RSD
C
      INTEGER               I
      INTEGER               J
      INTEGER               ERRCD
      INTEGER               RSTRUC
      INTEGER               IBAL
      INTEGER               IERR
      INTEGER               IND(SIZE2)
      INTEGER               MAXIT
      INTEGER               NINPS
      INTEGER               NOUTS
      INTEGER               NSTATS
C
      CHARACTER*(*)         SYSFIL
      CHARACTER*(*)         BALFIL
      CHARACTER*(*)         OUTFIL
      CHARACTER*(*)         ERRMSG
C
      INCLUDE			'dpcom.f'
C
C  Begin. Initialize IERR, MAXIT, RTOL, IBAL and RSTRUC.
C
      ERRCD = 0
      IERR = 0
      MAXIT = 5
      RTOL = 0.0D0
      IF (RHO.EQ.1.0D0) THEN
          RSTRUC = 0
      ELSE
          RSTRUC = 1
      END IF
C
C  Read in the system matrices.
C
      CALL INSYS(SYSFIL,NINPS,NOUTS,NSTATS,
     1           SIZE2,SIZE2,SIZE2,A,B,C,D,ERRCD)
C
       CLOSE(UNIT=UNIT1)
C
       IF (ERRCD.NE.0) THEN
          ERRMSG = 'LQR: Fatal Error from INSYS '//
     1             'accessing '//SYSFIL
          RETURN
       END IF
C
C  For the LQR solution of the Riccatti Equation, we
C  want the output matrix C to be the identity.
C
       DO 20 I=1,NSTATS
          DO 10 J=1,NSTATS
             C(J,I) = 0.0D0
 10       CONTINUE
          C(I,I) = 1.0D0
 20    CONTINUE
C
C   The R matrix is a diagonal matrix equal to RHO*I
C
       DO 40 I=1,NINPS
          DO 30 J=1,NINPS
             R(J,I) = 0.0D0
 30       CONTINUE
          R(I,I) = RHO
 40    CONTINUE
C
C   Read in the balancing matrix. P isn't used.
C
       OPEN(UNIT=UNIT1,FILE=BALFIL,STATUS='OLD',ERR=9990)
C
       DO 50 I=1,NINPS
          READ(UNIT1,*,ERR=9990) (P(I,J),J=1,NSTATS)
 50    CONTINUE
C
       DO 60 I=1,NSTATS
          READ(UNIT1,*,ERR=9990) (Q(I,J),J=1,NSTATS)
 60    CONTINUE
C
       CLOSE(UNIT=UNIT1,ERR=9990)
C
C   Compute the regulator solution.
C
       CALL CREG(SIZE2,NSTATS,NINPS,NSTATS,A,B,C,R,Q,G,F,Z,ACL,
     1          CQC,E,RI,RS,S,U,WK,ALFI,ALFR,BETA,CPERM,
     2          CSCALE,RSD,RTOL,MAXIT,
     3          RSTRUC,IBAL,IND,IERR,ERRMSG)
C
       IF (IERR.GT.0) THEN
          ERRCD = IERR
          RETURN
       END IF
C
C   Output the solution. Z is the optimal gain matrix.
C
       CALL OUTSYS(OUTFIL,NINPS,NINPS,NSTATS,SIZE2,
     1             SIZE2,SIZE2,A,B,Z,D,ERRCD)
C
       IF (ERRCD.NE.0) THEN
          ERRMSG = ' LQR: Fatal Error from OUTSYS.'
          RETURN
       END IF
C
C  Append the closed-loop matrix to the output file.
C
       DO 70 I=1,NSTATS
          WRITE(UNIT1,*,ERR=9991) (ACL(I,J),J=1,NSTATS)
 70    CONTINUE
C
C  Append the solution to the ARE to the output file.
C
       DO 80 I=1,NSTATS
          WRITE(UNIT1,*,ERR=9991) (F(I,J),J=1,NSTATS)
 80    CONTINUE
C
C  Append the closed-loop eigenvalues to the output file.
C
       DO 90 I=1,NSTATS
          WRITE(UNIT1,*,ERR=9991) ALFR(I),ALFI(I)
 90    CONTINUE
C
C  Close the output file.
C
       CLOSE(UNIT=UNIT1,ERR=9991)
C
C  End of program.
C
       RETURN
C
C  Error traps.
C
9990   ERRCD = 100 + ERRCD
       ERRMSG = 'LQR: Error reading Balance File '//BALFIL
       RETURN
9991   ERRCD = 200 + ERRCD
       ERRMSG = 'LQR: Error writing Output File '//OUTFIL
       RETURN
C
       END