      SUBROUTINE SIMPLX(A,M,N,MP,NP,M1,M2,M3,ICASE,IZROV,IPOSV)
      IMPLICIT NONE
C     PARAMETERS
      INTEGER MMAX
      DOUBLE PRECISION EPS
      PARAMETER(MMAX=100,EPS=1.D-6)
C     SCALAR ARGUMENTS
      INTEGER ICASE,M,M1,M2,M3,MP,N,NP
C     ARRAY ARGUMENTS
      INTEGER IZROV(N),IPOSV(M),L1(MMAX),L2(MMAX),L3(MMAX)
      DOUBLE PRECISION A(MP,NP)
C     LOCAL SCALARS
      INTEGER I,IP,IS,IR,K,KH,KP,M12,NL1,NL2
      DOUBLE PRECISION BMAX,Q1
      IF(M.NE.M1+M2+M3)PAUSE 'Bad input constraint counts.'
      NL1=N
      DO 11 K=1,N
        L1(K)=K
        IZROV(K)=K
11    CONTINUE
      NL2=M
      DO 12 I=1,M
        IF(A(I+1,1).LT.0.)PAUSE 'Bad input tableau.'
        L2(I)=I
        IPOSV(I)=N+I
12    CONTINUE
      DO 13 I=1,M2
        L3(I)=1
13    CONTINUE
      IR=0
      IF(M2+M3.EQ.0)GO TO 30
      IR=1
      DO 15 K=1,N+1
        Q1=0.D0
        DO 14 I=M1+1,M
          Q1=Q1+A(I+1,K)
14      CONTINUE
        A(M+2,K)=-Q1
15    CONTINUE
10    CALL SIMP1(A,MP,NP,M+1,L1,NL1,0,KP,BMAX)
      IF(BMAX.LE.EPS.AND.A(M+2,1).LT.-EPS)THEN
        ICASE=-1
        RETURN
      ELSE IF(BMAX.LE.EPS.AND.A(M+2,1).LE.EPS)THEN
        M12=M1+M2+1
        IF(M12.LE.M)THEN
          DO 16 IP=M12,M
            IF(IPOSV(IP).EQ.IP+N)THEN
              CALL SIMP1(A,MP,NP,IP,L1,NL1,1,KP,BMAX)
              IF(BMAX.GT.0.)GO TO 1
            ENDIF
16        CONTINUE
        ENDIF
        IR=0
        M12=M12-1
        IF(M1+1.GT.M12)GO TO 30
        DO 18 I=M1+1,M12
          IF(L3(I-M1).EQ.1)THEN
            DO 17 K=1,N+1
              A(I+1,K)=-A(I+1,K)
17          CONTINUE
          ENDIF
18      CONTINUE
        GO TO 30
      ENDIF
      CALL SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
      IF(IP.EQ.0)THEN
        ICASE=-1
        RETURN
      ENDIF
1     CALL SIMP3(A,MP,NP,M+1,N,IP,KP)
      IF(IPOSV(IP).GE.N+M1+M2+1)THEN
        DO 19 K=1,NL1
          IF(L1(K).EQ.KP)GO TO 2
19      CONTINUE
2       NL1=NL1-1
        DO 21 IS=K,NL1
          L1(IS)=L1(IS+1)
21      CONTINUE
      ELSE
        IF(IPOSV(IP).LT.N+M1+1)GO TO 20
        KH=IPOSV(IP)-M1-N
        IF(L3(KH).EQ.0)GO TO 20
        L3(KH)=0
      ENDIF
      A(M+2,KP+1)=A(M+2,KP+1)+1.D0
      DO 22 I=1,M+2
        A(I,KP+1)=-A(I,KP+1)
22    CONTINUE
20    IS=IZROV(KP)
      IZROV(KP)=IPOSV(IP)
      IPOSV(IP)=IS
      IF(IR.NE.0)GO TO 10
30    CALL SIMP1(A,MP,NP,0,L1,NL1,0,KP,BMAX)
      IF(BMAX.LE.0.)THEN
        ICASE=0
        RETURN
      ENDIF
      CALL SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
      IF(IP.EQ.0)THEN
        ICASE=1
        RETURN
      ENDIF
      CALL SIMP3(A,MP,NP,M,N,IP,KP)
      GO TO 20
      END

      SUBROUTINE SIMP1(A,MP,NP,MM,LL,NLL,IABF,KP,BMAX)
      IMPLICIT NONE
C     SCALAR ARGUMENTS
      INTEGER IABF,KP,MP,MM,NLL,NP
      DOUBLE PRECISION BMAX
C     ARRAY ARGUMENTS
      INTEGER LL(NP)
      DOUBLE PRECISION A(MP,NP)
C     LOCAL SCALARS
      INTEGER K
      DOUBLE PRECISION TEST
      KP=LL(1)
      BMAX=A(MM+1,KP+1)
      IF(NLL.LT.2)RETURN
      DO 11 K=2,NLL
        IF(IABF.EQ.0)THEN
          TEST=A(MM+1,LL(K)+1)-BMAX
        ELSE
          TEST=ABS(A(MM+1,LL(K)+1))-ABS(BMAX)
        ENDIF
        IF(TEST.GT.0.)THEN
          BMAX=A(MM+1,LL(K)+1)
          KP=LL(K)
        ENDIF
11    CONTINUE
      RETURN
      END

      SUBROUTINE SIMP2(A,M,N,MP,NP,L2,NL2,IP,KP,Q1)
      IMPLICIT NONE
C     PARAMETERS
      DOUBLE PRECISION EPS
      PARAMETER (EPS=1.D-6)
C     SCALAR ARGUMENTS
      INTEGER IP,KP,M,MP,N,NL2,NP
      DOUBLE PRECISION Q1
C     ARRAY ARGUMENTS
      INTEGER L2(MP)
      DOUBLE PRECISION A(MP,NP)
C     LOCAL SCALARS
      INTEGER I,II,K
      DOUBLE PRECISION Q,Q0,QP
      IP=0
      IF(NL2.LT.1)RETURN
      DO 11 I=1,NL2
        IF(A(L2(I)+1,KP+1).LT.-EPS)GO TO 2
11    CONTINUE
      RETURN
2     Q1=-A(L2(I)+1,1)/A(L2(I)+1,KP+1)
      IP=L2(I)
      IF(I+1.GT.NL2)RETURN
      DO 13 I=I+1,NL2
        II=L2(I)
        IF(A(II+1,KP+1).LT.-EPS)THEN
          Q=-A(II+1,1)/A(II+1,KP+1)
          IF(Q.LT.Q1)THEN
            IP=II
            Q1=Q
          ELSE IF (Q.EQ.Q1) THEN
            DO 12 K=1,N
              QP=-A(IP+1,K+1)/A(IP+1,KP+1)
              Q0=-A(II+1,K+1)/A(II+1,KP+1)
              IF(Q0.NE.QP)GO TO 6
12          CONTINUE
6           IF(Q0.LT.QP)IP=II
          ENDIF
        ENDIF
13    CONTINUE
      RETURN
      END

      SUBROUTINE SIMP3(A,MP,NP,I1,K1,IP,KP)
      IMPLICIT NONE
c     SCALAR ARGUMENTS
      INTEGER I1,IP,K1,KP,MP,NP
C     ARRAY ARGUMENTS
      DOUBLE PRECISION A(MP,NP)
C     LOCAL SCALARS
      INTEGER II,KK
      DOUBLE PRECISION PIV
      PIV=1.D0/A(IP+1,KP+1)
      IF(I1.GE.0)THEN
        DO 12 II=1,I1+1
          IF(II-1.NE.IP)THEN
            A(II,KP+1)=A(II,KP+1)*PIV
            DO 11 KK=1,K1+1
              IF(KK-1.NE.KP)THEN
                A(II,KK)=A(II,KK)-A(IP+1,KK)*A(II,KP+1)
              ENDIF
11          CONTINUE
          ENDIF
12      CONTINUE
      ENDIF
      DO 13 KK=1,K1+1
        IF(KK-1.NE.KP)A(IP+1,KK)=-A(IP+1,KK)*PIV
13    CONTINUE
      A(IP+1,KP+1)=PIV
      RETURN
      END


