!****************************************************************
!*           LINEAR PROGRAMMING: THE SIMPLEX METHOD             *
!* ------------------------------------------------------------ *
!* SAMPLE RUN:                                                  *
!* Maximize z = x1 + x2 + 3x3 -0.5x4 with conditions:           * 
!*          x1  + 2x3 <= 740                                    *
!*          2x2 - 7x4 <= 0                                      *
!*          x2  - x3 + 2x4 >= 0.5                               *
!*          x1 + x2 + x3 +x4 = 9                                *
!*          and all x's >=0.                                    *
!*                                                              *
!* Number of variables in E.F.: 4                               *
!* Number of <= inequalities..: 2                               *
!* Number of >= inequalities..: 1                               *
!* Number of = equalities.....: 1                               *
!* Input Economic Function:                                     *
!* Coefficient # 1: 1                                           *
!* Coefficient # 2: 1                                           *
!* Coefficient # 3: 3                                           *
!* Coefficient # 4: -0.5                                        *
!* Constant term..: 0                                           *
!* Input constraint # 1:                                        *
!* Coefficient # 1: 1                                           *
!* Coefficient # 2: 0                                           *
!* Coefficient # 3: 2                                           *
!* Coefficient # 4: 0                                           *
!* Constant term..: 740                                         *
!* Input constraint # 2:                                        *
!* Coefficient # 1: 0                                           *
!* Coefficient # 2: 2                                           *
!* Coefficient # 3: 0                                           *
!* Coefficient # 4: -7                                          *
!* Constant term..: 0                                           *
!* Input constraint # 3:                                        *
!* Coefficient # 1: 0                                           *
!* Coefficient # 2: 1                                           *
!* Coefficient # 3: -1                                          *
!* Coefficient # 4: 2                                           *
!* Constant term..: 0.5                                         *
!* Input constraint # 4:                                        *
!* Coefficient # 1: 1                                           *
!* Coefficient # 2: 1                                           *
!* Coefficient # 3: 1                                           *
!* Coefficient # 4: 1                                           *
!* Constant term..: 9                                           *
!*                                                              *
!*  Input Table:                                                *
!*    0.00    1.00    1.00    3.00   -0.50                      *
!*  740.00   -1.00    0.00   -2.00    0.00                      *
!*    0.00    0.00   -2.00    0.00    7.00                      *
!*    0.50    0.00   -1.00    1.00   -2.00                      *
!*    9.00   -1.00   -1.00   -1.00   -1.00                      *
!*                                                              *
!*  Maximum of E.F. =    17.02500                               *
!*  X 1 =    0.000000                                           *
!*  X 2 =    3.325000                                           *
!*  X 3 =    4.725000                                           *
!*  X 4 =    0.950000                                           *
!*                                                              *
!* ------------------------------------------------------------ *
!* Reference: "Numerical Recipes by W.H. Press, B. P. Flannery, *
!*             S.A. Teukolsky and W.T. Vetterling, Cambridge    *
!*             University Press, 1986".                         *
!*                                                              *
!*                       F90 Release 1.0 By J-P Moreau, Paris   *
!*                                (www.jpmoreau.fr)             *
!**************************************************************** 
PROGRAM TEST_SIMPLEX
implicit none
integer I,ICASE,J,M,N,M1,M2,M3
integer, parameter :: MMAX=25, NMAX=25
REAL A(MMAX,NMAX)
INTEGER IPOSV(MMAX), IZROV(NMAX)
REAL R

print *,' '
write(*,10,advance='no'); read *, N
write(*,20,advance='no'); read *, M1
write(*,30,advance='no'); read *, M2
write(*,40,advance='no'); read *, M3

M=M1+M2+M3  !Total number of constraints

A=0.
print *,'Input Economic Function:'
do i=2,N+1
  write(*,50,advance='no') i-1; read *,A(1,i)
end do
write(*,51,advance='no'); read *,A(1,1)
!input constraints    
do i=1, M
  write(*,60) i
  do j=2,N+1
    write(*,50,advance='no') j-1; read *, R
	A(i+1,j) = -R
  end do
  write(*,51,advance='no'); read *,A(i+1,1)   
end do

print *,' '
print *,' Input Table:'
write(*,100) ((A(I,J),J=1,N+1),I=1,M+1)

call simplx(A,M,N,MMAX,NMAX,M1,M2,M3,ICASE,IZROV,IPOSV)

print *,' '
print *,' Maximum of E.F. = ', A(1,1)

do I=1, N
  do J=1, M
    if (IPOSV(J).eq.I) then
	  write(*,110) I, A(J+1, 1)
	  goto 3
    end if
  end do
  write(*,110) I, 0.0 
3 end do
  print *,' '

10  format(' Number of variables in E.F.: ')
20  format(' Number of <= inequalities..: ')
30  format(' Number of >= inequalities..: ')
40  format(' Number of = equalities.....: ')
50  format(' Coefficient #',I2,': ')
51  format(' Constant term..: ')
60  format(' Input constraint #',I2,': ')

100 format(5F8.2)
110 format('  X',I2,' = ',F12.6)

END

SUBROUTINE simplx(a,m,n,mp,np,m1,m2,m3,icase,izrov,iposv)
implicit none
INTEGER icase,m,m1,m2,m3,mp,n,np,iposv(m),izrov(n) 
REAL a(mp,np) 
INTEGER, PARAMETER :: MMAX=100, NMAX=100
REAL, PARAMETER :: EPS=1.E-6
!---------------------------------------------------------------------------------------- 
! USES simp1,simp2,simp3 
!Simplex method for linear programming. Input parameters a, m, n, mp, np, m1, m2, and m3, 
!and output parameters a, icase, izrov, and iposv are described above. 
!Parameters: MMAX is the maximum number of constraints expected; NMAX is the maximum 
!number of variables expected; EPS is the absolute precision, which should be adjusted to 
!the scale of your variables. 
!----------------------------------------------------------------------------------------
INTEGER i,ip,ir,is,k,kh,kp,m12,nl1,nl2,l1(NMAX),l2(MMAX),l3(MMAX) 
REAL bmax,q1 
if(m.ne.m1+m2+m3) pause ' Bad input constraint counts in simplx.' 
nl1=n 
do k=1,n 
  l1(k)=k    !Initialize index list of columns admissible for exchange.
  izrov(k)=k !Initially make all variables right-hand. 
end do 
nl2=m
do i=1,m 
  if(a(i+1,1).lt.0.)pause ' Bad input tableau in simplx, Constants bi must be nonnegative.'
  l2(i)=i
  iposv(i)=n+i 
!-------------------------------------------------------------------------------------------------
!Initial left-hand variables. m1 type constraints are represented by having their slackv ariable 
!initially left-hand, with no artificial variable. m2 type constraints have their slack 
!variable initially left-hand, with a minus sign, and their artificial variable handled implicitly 
!during their first exchange. m3 type constraints have their artificial variable initially 
!left-hand.
!-------------------------------------------------------------------------------------------------
end do
do i=1,m2
  l3(i)=1
end do
ir=0
if(m2+m3.eq.0) goto 30 !The origin is a feasible starting solution. Go to phase two. 
ir=1
do k=1,n+1             !Compute the auxiliary objective function. 
  q1=0. 
  do i=m1+1,m
    q1=q1+a(i+1,k) 
  end do 
  a(m+2,k)=-q1 
end do 
10 call simp1(a,mp,np,m+1,l1,nl1,0,kp,bmax) !Find max. coeff. of auxiliary objective fn 
if(bmax.le.EPS.and.a(m+2,1).lt.-EPS)then 
  icase=-1        !Auxiliary objective function is still negative and can’t be improved, 
  return          !hence no feasible solution exists. 
else if(bmax.le.EPS.and.a(m+2,1).le.EPS)then 
!Auxiliary objective function is zero and can’t be improved; we have a feasible starting vector. 
!Clean out the artificial variables corresponding to any remaining equality constraints by 
!goto 1’s and then move on to phase two by goto 30.
  m12=m1+m2+1
  if (m12.le.m) then
    do ip=m12,m 
      if(iposv(ip).eq.ip+n)then !Found an artificial variable for an equalityconstraint. 
        call simp1(a,mp,np,ip,l1,nl1,1,kp,bmax) 
        if(bmax.gt.EPS) goto 1  !Exchange with column corresponding to maximum 
      end if                    !pivot element in row. 
    end do	                     
  end if
  ir=0
  m12=m12-1
  if (m1+1.gt.m12) goto 30 
  do i=m1+1,m1+m2               !Change sign of row for any m2 constraints 
                                !still present from the initial basis. 
    if(l3(i-m1).eq.1)then 
      do k=1,n+1 
        a(i+1,k)=-a(i+1,k) 
      end do 
    end if 
  end do 
  goto 30                        !Go to phase two. 
end if 
call simp2(a,m,n,mp,np,l2,nl2,ip,kp,q1) !Locate a pivot element (phase one). 
if(ip.eq.0)then                  !Maximum of auxiliary objective function is 
  icase=-1                       !unbounded, so no feasible solution exists. 
  return 
end if 

1 call simp3(a,mp,np,m+1,n,ip,kp) 

!Exchange a left- and a right-hand variable (phase one), then update lists. 
if(iposv(ip).ge.n+m1+m2+1)then   !Exchanged out an artificial variable for an 
                                 !equality constraint. Make sure it stays 
  do k=1,nl1                     !out by removing it from the l1 list. 
    if(l1(k).eq.kp) goto 2 
  end do 

2 nl1=nl1-1 
  do is=k,nl1 
    l1(is)=l1(is+1) 
  end do 
else
  if(iposv(ip).lt.n+m1+1) goto 20
  kh=iposv(ip)-m1-n 
  if(l3(kh).eq.0) goto 20      !Exchanged out an m2 type constraint. 
  l3(kh)=0                     !If it’s the first time, correct the pivot column 
                               !or the minus sign and the implicit 
end if                         !artificial variable. 

a(m+2,kp+1)=a(m+2,kp+1)+1. 
do i=1,m+2 
  a(i,kp+1)=-a(i,kp+1) 
end do 

20 is=izrov(kp)                !Update lists of left- and right-hand variables. 
izrov(kp)=iposv(ip) 
iposv(ip)=is 
if (ir.ne.0) goto 10           !if still in phase one, go back to 10. 
!End of phase one code for finding an initial feasible solution. Now, in phase two, optimize it. 
30 call simp1(a,mp,np,0,l1,nl1,0,kp,bmax) !Test the z-row for doneness. 
if(bmax.le.EPS)then            !Done. Solution found. Return with the good news. 
  icase=0 
  return 
end if

call simp2(a,m,n,mp,np,l2,nl2,ip,kp,q1)  !Locate a pivot element (phase two). 

if(ip.eq.0)then                !Objective function is unbounded. Report and return. 
  icase=1 
  return 
end if 

call simp3(a,mp,np,m,n,ip,kp)  !Exchange a left- and a right-hand variable (phase two), 
goto 20                        !update lists of left- and right-hand variables and 
END                            !return for another iteration.

!The preceding routine makes use of the following utility subroutines: 

SUBROUTINE simp1(a,mp,np,mm,ll,nll,iabf,kp,bmax)
implicit none
INTEGER iabf,kp,mm,mp,nll,np,ll(np) 
REAL bmax,a(mp,np) 
!Determines the maximum of those elements whose index is contained in the supplied list 
!ll, either with or without taking the absolute value, as flagged by iabf. 
INTEGER k 
REAL test 
  kp=ll(1) 
  bmax=a(mm+1,kp+1)
  if(nll.lt.2) return 
  do 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 
  end do 
return 
END 

SUBROUTINE simp2(a,m,n,mp,np,l2,nl2,ip,kp,q1)
implicit none
INTEGER ii,ip,kp,mp,l2(mp),m,n,nl2,np 
REAL a(mp,np),EPS 
PARAMETER (EPS=1.e-6) 
!Locate a pivot element, taking degeneracy into account. 
INTEGER i,k 
REAL q,q0,q1,qp 
ip=0 
if(nl2.lt.1) return
do i=1,nl2 
  if(a(i+1,kp+1).lt.-EPS) goto 2 
end do 
return  !No possible pivots. Return with message. 
2 q1=-a(l2(i)+1,1)/a(l2(i)+1,kp+1) 
ip=l2(i)
if(i+1.gt.nl2) return 
do 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 !We have a degeneracy.
      do 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)goto 6 
      end do 
6     if(q0.lt.qp) ip=ii 
    end if 
  end if 
end do 
return 
END 

SUBROUTINE simp3(a,mp,np,i1,k1,ip,kp)
implicit none
INTEGER i1,ip,k1,kp,mp,np 
REAL a(mp,np) 
!Matrix operations to exchange a left-hand and right-hand variable (see text). 
INTEGER ii,kk 
REAL piv 
piv=1./a(ip+1,kp+1)
if (i1.ge.0) then 
  do ii=1,i1+1 
    if(ii-1.ne.ip)then 
      a(ii,kp+1)=a(ii,kp+1)*piv 
      do kk=1,k1+1 
        if(kk-1.ne.kp)then 
          a(ii,kk)=a(ii,kk)-a(ip+1,kk)*a(ii,kp+1) 
        end if 
      end do 
    end if 
  end do 
end if
do kk=1,k1+1 
  if(kk-1.ne.kp) a(ip+1,kk)=-a(ip+1,kk)*piv 
end do 
a(ip+1,kp+1)=piv 
return 
END 

!end of file tsimplex.f90