!      SUBROUTINE RLLAV  (XY,IXY,N,M,IOPT,BETA,SUMRE,ITER,IRANK)
      SUBROUTINE RLLAV (XY,N,M,IOPT, BETA)

      INTEGER :: IXY,N,M,IOPT,ITER,IRANK,IWK(N),IER
!      Parameter (Nmax=30, Mmax=4)
      Real :: XY(Nmax, Mmax+5)
      REAL :: BETA(Mmax+1)
!      REAL               XY(IXY,1),BETA(1),SUMRE,WK(1)
      REAL SUMRE,WK(2*M+4)
!                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER            IN,IOUT,I,J,KL,KOUNT,KR,K,L,M1,M2,M3,M5,MM,N1,N2
      REAL               AMAX,AMIN,BIG,D,FMM,ONE,PIVOT,TOLER,TWO,ZERO
      LOGICAL            STAGE,TEST
      DOUBLE PRECISION   SUM
      DATA               TOLER/1.1921E-07/
      DATA               BIG/1.7E+38/
      DATA               ONE/1.0/,ZERO/0.0/,TWO/2.0/
!                                  FIRST EXECUTABLE STATEMENT
      MM = M
      N1 = N+1
      N2 = N+2
      M1 = M+1
      M2 = M+2
      M5 = M+5
!      IOPT = 0
!      print*,XY

      IF (IOPT.NE.0) GO TO 10
      MM = M+1
      DO 5 I=1,N
         XY(I,M2) = XY(I,M1)
         XY(I,M1) = ONE
    5 CONTINUE
      M1 = MM+1
      M2 = MM+2
   10 M3 = MM+3
      FMM = MM
      DO 15 J=1,MM
         WK(M1+J) = J
         BETA(J) = ZERO
   15 CONTINUE
      DO 30 I=1,N
         XY(I,M2) = MM+I
         XY(I,M3) = XY(I,M1)
         IF (XY(I,M3).GE.ZERO) GO TO 25
         DO 20 J=1,M2
            XY(I,J) = -XY(I,J)
   20    CONTINUE
   25    XY(I,M5) = ZERO
   30 CONTINUE
!                                  COMPUTE THE MARGINAL COSTS
      DO 40 J=1,M1
         SUM = 0.D0
         DO 35 I=1,N
            SUM = SUM+XY(I,J)
   35    CONTINUE
         WK(J) = SUM
   40 CONTINUE
!                                  STAGE I
!                                    DETERMINE THE VECTOR TO
!                                    ENTER THE BASIS
      STAGE = .TRUE.
      KOUNT = 0
      KR = 1
      KL = 1
   45 AMAX = -1
      DO 50 J=KR,MM
         IF (ABS(WK(M1+J)).GT.FMM) GO TO 50
         D = ABS(WK(J))
         IF (D.LE.AMAX) GO TO 50
         AMAX = D
         IN = J
   50 CONTINUE
      IF (WK(IN).GE.ZERO) GO TO 60
      DO 55 I=1,N
         XY(I,IN) = -XY(I,IN)
   55 CONTINUE
      WK(IN) = -WK(IN)
      WK(M1+IN) = -WK(M1+IN)
!                                  DETERMINE THE VECTOR TO LEAVE THE
!                                    BASIS
   60 K = 0
      DO 65 I=KL,N
         D = XY(I,IN)
         IF (D.LE.TOLER) GO TO 65
         K = K+1
         XY(K,M3) = XY(I,M1)/D
         IWK(K) = I
         TEST = .TRUE.
   65 CONTINUE
   70 IF (K.GT.0) GO TO 75
      TEST = .FALSE.
      GO TO 85
   75 AMIN = BIG
      DO 80 I=1,K
         IF (XY(I,M3).GE.AMIN) GO TO 80
         J = I
         AMIN = XY(I,M3)
         IOUT = IWK(I)
   80 CONTINUE
      XY(J,M3) = XY(K,M3)
      IWK(J) = IWK(K)
      K = K-1
!                                  CHECK FOR LINEAR DEPENDENCE IN
!                                    STAGE I
   85 IF (TEST.OR..NOT.STAGE) GO TO 95
      DO 90 I=1,N
         D = XY(I,KR)
         XY(I,KR) = XY(I,IN)
         XY(I,IN) = D
   90 CONTINUE
      D = WK(KR)
      WK(KR) = WK(IN)
      WK(IN) = D
      D = WK(M1+KR)
      WK(M1+KR) = WK(M1+IN)
      WK(M1+IN) = D
      KR = KR+1
      GO TO 145
   95 IF (TEST) GO TO 100
      IER = 129
      GO TO 9000
  100 PIVOT = XY(IOUT,IN)
      IF (WK(IN)-PIVOT-PIVOT.LE.TOLER) GO TO 110
      DO 105 J=KR,M1
         D = XY(IOUT,J)
         WK(J) = WK(J)-D-D
         XY(IOUT,J) = -D
  105 CONTINUE
      XY(IOUT,M2) = -XY(IOUT,M2)
      GO TO 70
!                                  PIVOT ON XY(IOUT,IN)
  110 DO 115 J=KR,M1
         IF (J.EQ.IN) GO TO 115
         XY(IOUT,J) = XY(IOUT,J)/PIVOT
  115 CONTINUE
      DO 125 I=1,N
         IF (I.EQ.IOUT) GO TO 125
         D = XY(I,IN)
         DO 120 J=KR,M1
            IF (J.EQ.IN) GO TO 120
            XY(I,J) = XY(I,J)-D*XY(IOUT,J)
  120    CONTINUE
  125 CONTINUE
      DO 130 J=KR,M1
         IF (J.EQ.IN) GO TO 130
         WK(J) = WK(J)-WK(IN)*XY(IOUT,J)
  130 CONTINUE
      DO 135 I=1,N
         IF (I.EQ.IOUT) GO TO 135
         XY(I,IN) = -XY(I,IN)/PIVOT
  135 CONTINUE
      WK(IN) = -WK(IN)/PIVOT
      XY(IOUT,IN) = ONE/PIVOT
      D = XY(IOUT,M2)
      XY(IOUT,M2) = WK(M1+IN)
      WK(M1+IN) = D
      KOUNT = KOUNT+1
      IF (.NOT.STAGE) GO TO 150
!                                  INTERCHANGE ROWS IN STAGE I
      KL = KL+1
      DO 140 J=KR,M2
         D = XY(IOUT,J)
         XY(IOUT,J) = XY(KOUNT,J)
         XY(KOUNT,J) = D
  140 CONTINUE
  145 IF (KOUNT+KR.NE.M1) GO TO 45
!                                  STAGE II
      STAGE = .FALSE.
!                                  DETERMINE THE VECTOR TO ENTER THE
!                                    BASIS
  150 AMAX = -BIG
      DO 160 J=KR,MM
         D = WK(J)
         IF (D.GE.ZERO) GO TO 155
         IF (D.GT.(-TWO)) GO TO 160
         D = -D-TWO
  155    IF (D.LE.AMAX) GO TO 160
         AMAX = D
         IN = J
  160 CONTINUE
      IF (AMAX.LE.TOLER) GO TO 170
      IF (WK(IN).GT.ZERO) GO TO 60
      DO 165 I=1,N
         XY(I,IN) = -XY(I,IN)
  165 CONTINUE
      WK(IN) = -WK(IN)-TWO
      WK(M1+IN) = -WK(M1+IN)
      GO TO 60
!                                  PREPARE OUTPUT
  170 L = KL-1
      DO 180 I=1,L
         IF (XY(I,M1).GE.ZERO) GO TO 180
         DO 175 J=KR,M2
            XY(I,J) = -XY(I,J)
  175    CONTINUE
  180 CONTINUE
      IER = 33
      IF (KR.NE.1) GO TO 190
      DO 185 J=1,MM
         D = ABS(WK(J))
         IF (D.LE.TOLER.OR.TWO-D.LE.TOLER) GO TO 190
  185 CONTINUE
      IER = 0
  190 DO 205 I=1,N
         K = XY(I,M2)
         D = XY(I,M1)
         IF (K.GT.0) GO TO 195
         K = -K
         D = -D
  195    IF (I.GE.KL) GO TO 200
         BETA(K) = D
!         print*,"K, BETA=",K,BETA(K)
         GO TO 205
  200    K = K-MM
         XY(K,M5) = D
  205 CONTINUE
      ITER = KOUNT
      IRANK = M1-KR
      SUM = 0.D0
      DO 210 I=KL,N
         SUM = SUM+XY(I,M1)
  210 CONTINUE
      SUMRE = SUM
      IF (IER.EQ.0) GO TO 9005
 9000 CONTINUE
!      CALL UERTST (IER,'RLLAV ')
      Print*,"Sai dau do trong RLLAV"
 9005 RETURN
      END SUBROUTINE RLLAV
!--------------------- End of this file ----------------------