!  Basic routines for simplex method
!  File:  simplex.f90		header:  simplex.h

!  Tableau:
!
!	LHS  |  Const.  X(v(1))	X(v(2))	. . .	X(v(N))
!    --------+--------------------------------------------
!    Objctv  |	A(0,0)	A(0,1)	A(0,2)	. . .	A(0,N)
!    X(u(1)) |	A(1,0)	A(1,1)	A(1,2)	. . .	A(1,N)
!    X(u(2)) |	A(2,0)	A(2,1)	A(2,2)	. . .	A(2,N)
!       :    |     :       :       :               :
!       :    |     :       :       :               :
!    X(u(M)) |	A(M,0)	A(M,1)	A(M,2)	. . .	A(M,N)
!
!  Requires A(1:M,0) >= 0 .
!
!  Each row is read as
!
!	LHS = A(i,0) + A(i,1)*X(v(1)) + A(i,2)*X(v(2)) +...+ A(i,N)*X(v(N))


!  Attempt maximization of the objective.  Leaves updated tableau in
!  A, U, V.  Outcome is reported in Iout as
!
!	Iout = 0	Maximal solution found
!	Iout = 1	Objective can be arbitrarily large
!	Iout = 2	Failed due to unmanageable degeneracy
!
!  Icyc reports number of cycles.

      subroutine simplex(M, N, A, U, V, Iout, Icyc)
      integer M, N, Iout, Icyc
      real,    dimension(0:M,0:N) :: A
      integer, dimension(0:M)     :: U
      integer, dimension(0:N)     :: V

      INCLUDE 'simprnt.h'

      integer, parameter :: dlim = 100	! Max consecutive degeneracies
      real,    parameter :: eps = 1E-20	! Tolerance for 0 test

      real,    dimension(0:M,0:N) :: Cn, Pr, Gn, Xi	! Working arrays
      integer, dimension(0:M,0:N) :: I, J	! Coordinates
      logical, dimension(0:M,0:N) :: Lc, Lp	! Flag coefficients, pivots
      integer, dimension(1:2)     :: maxl	! For maxloc

      integer dcount		! Counts steps with degenerate update
      integer itmp		! For swap
      integer p, q		! Pivot coordinates
      real    pv		! Pivot coefficient ( < 0.0 )
      real    mxobj             ! Maximum objective coefficient

      I      = SPREAD([0:M], 2, N+1)		! Coordinates
      J      = SPREAD([0:N], 1, M+1)
      Lc     = (I > 0) .and. (J > 0)
      dcount = 0
      Iout   = 0				! OK so far
      Icyc   = 0
!      print *, "M, N = ", M, N, "  Initial A:"
!      call prflt(M, N, A)
!      print *, "Lc:"
!      call prlog(M, N, Lc)
!
!  Repeat until optimum or failure
!  f90 compiler had trouble with maxval inside while test - 21 Jan 1997
      mxobj = maxval(A(0, 1:N))
      do while ( (mxobj > 0.0) .and. (dcount < dlim) )
!  Locate pivot
	Pr = SPREAD(A(0,:), 1, M+1)		! Prices	   
	Cn = SPREAD(A(:,0), 2, N+1)		! Constants
	Gn = -1.0
	Lp = .false.
	WHERE ( Pr > 0.0 .and. A < 0.0 .and. Lc )	! Potential pivots
	  Gn = -Cn/A				! Allowed increase, >= 0.0
	  Lp = .true.				! Potential pivot positions
	end WHERE
!	print *, "Lp:"
!	call prlog(M, N, Lp)
	if ( ANY(Pr(0,1:N) > 0.0 .and. &
	          .not. ANY(Lp(:,1:N), 1)) ) then
	  Iout = 1				! Unbounded
	  go to 99				! Exit
	end if
	Xi = SPREAD(MINVAL(Gn, 1, Lp), 1, M+1)	! Allowed for column
	Lp = Lp .and. Gn == Xi			! Legitimate pivots
	WHERE ( Lp )				! Potential pivots
	  Gn = Gn*Pr				! Gain in objective
	ELSEWHERE
	  Gn = -1.0				! Not a possible pivot
	end WHERE
	if ( maxval(array=Gn, mask=Lp) < eps ) &  ! Degenerate update
	then
	  dcount = dcount + 1			! To insure termination
	  maxl = maxloc(Pr, Lp)			! Pivot with max price
	else					! Normal update
	  dcount = 0				! Pivot with max gain
	  maxl = maxloc(Gn, Lp)			! Pivot location
	end if
	p  = maxl(1) - 1			! Pivot row (0-based)
	q  = maxl(2) - 1			! Pivot col (0-based)
	pv = abs(A(p,q))			! | pivot value |
!	print 31, p, q, pv
! 31	FORMAT (" Pivot:  (", I2, ", ", I1, ")  ", F12.6)
	A(p,:) = A(p,:)/pv			! Divide pivot row
	Pr = SPREAD(A(p,:), 1, M+1)		!   and spread 
	Cn = SPREAD(A(:,q), 2, N+1)		! Spread pivot col
	WHERE ( I <> p .and. J <> q ) &		! Except pivot row, col
	  A = A + Pr*Cn
	A(:,q) = -A(:,q)/pv			! Update pivot col.
	A(p,q) = -1.0/pv			! New pivot entry
	WHERE ( A(1:M,0) < 0.0 ) A(1:M,0) = 0.0	! For rounding errors
	itmp   = U(p)				! Swap variables
	U(p)   = V(q)
	V(q)   = itmp
	Icyc   = Icyc + 1			! One more cycle
!	print *, "Current A"
!	call prflt(M, N, A)
        mxobj = maxval(A(0, 1:N))
      end do ! while ( maxval > 0 .and. dcount < dlim )
      if ( dcount == dlim ) Iout = 2		! Degeneracy failure
 99   continue					! For exit
      end subroutine simplex
