!  Attempt to find basic, feasible solution and then optimum of
!  standard LP problem.
!
!  File:  lpsolve.f90		Header:  lpsolve.h

!  Problem form:
!
!  Maximize A(0,1:N).X(1:N)  subject to  X(1:N) >= 0 and
!
!	A(i,1:N).X(1:N) ## A(i,0)  for 1 <= i <= M
!
!  where, for each constraint, ## is determined from
!
!		R(i) | ##
!		-----+----
!		 -1  | >=
!		  0  | =
!		 +1  | <=
!
!  It is required that A(1:M,0) >= 0.  If A(0,0) <> 0 this will be incoporated
!  in the objective function.
!
!  Outcome is reported in Iout:
!
!	0 - optimal solution
!	1 - objective not bounded above
!	2 - optimization (phase 2) failed due to degeneracy
!	3 - infeasible
!	4 - phase 1 failed with unbounded feasibility
!	5 - phase 1 failed due to degeneracy
!
!  Ic1, Ic2 have the number of simplex cycles in phase1, phase2.
!
!  The final tableau is returned in A, U, V as
!
!	LHS  |  Const.  Z(V(1))	Z(V(2))	. . .	Z(V(N))
!    --------+--------------------------------------------
!    Objctv  =	A(0,0)	A(0,1)	A(0,2)	. . .	A(0,N)
!    Z(U(1)) =	A(1,0)	A(1,1)	A(1,2)	. . .	A(1,N)
!    Z(U(2)) =	A(2,0)	A(2,1)	A(2,2)	. . .	A(2,N)
!       :    :     :       :       :               :
!       :    :     :       :       :               :
!    Z(U(M)) =	A(M,0)	A(M,1)	A(M,2)	. . .	A(M,N)
!
!  where Z(U(i)) is X(U(i)) if U(i) > 0, the slack or artificial variable
!  associated with the original constraint -U(i) if U(i) < 0.  Similarly
!  for Z(V(i)).

      subroutine lpsolve(M, N, A, R, U, V, Iout, Ic1, Ic2)
      integer M, N, Iout, Ic1, Ic2
      real,    dimension(0:M,0:N) :: A
      integer, dimension(0:M)     :: R, U
      integer, dimension(0:N)     :: V

      real,    parameter :: eps = 1E-10		! Tolerance for 0 test

      INCLUDE 'simplex.h'
      INCLUDE 'simplex1.h'
      INCLUDE 'simprnt.h'

      integer, dimension(0:M,0:N) :: Lx
      logical, dimension(0:M,0:N) :: La, Lw
      real,    dimension(0:M)     :: C1
      real,    dimension(0:N)     :: C2		! The initial objective
      logical, dimension(0:N)     :: Vw

      A(1:M,1:N) = -A(1:M,1:N)			! Coeff's for phase 1
      C2   =  A(0,:)				! Remember objective
      R(0) =  1					! For convenience
      U    = -[0:M]				! Slacks and artificials
      where ( R == 0 ) U = U - M		! Equations, very negative code
      V    =  [0:N]				! Problem variables
      Lx   = SPREAD(R, 2, N+1)			! Rows with artificials
!      print *, "Lx:"
!      call printg(M, N, Lx)
      where ( Lx == -1 ) A = -A			! Neg. slacks -> artificial
!   Objective is -(Sum of infeasibilities)
      A(0,:) = SUM(A, 1, Lx == -1) - SUM(A, 1, Lx == 0)
      call simplex1(M, N, A, U, V, Iout, Ic1)	! Find feasible solution
      if ( Iout > 0 .or. A(0,0) < -eps) then	! Can't continue
	Iout = Iout + 3				! Error code
	Ic2 = 0					! Didn't try phase 2
	go to 99				! Exit
      end if
!      print *, "U(0:M):";  print *, U
!      print *, "V(0:N):";  print *, V
!  Parasoft f90 has trouble with the logical array expressions.  It's
!  not entirely clear that the following is compiled correctly.
      La = SPREAD(U < -M, 2, N+1)		! Artificials in basis
!      print *, "La:";  call prlog(M, N, La)
      Lw = SPREAD(ANY(La .and. A > eps, 1), 1, M+1)	! Var's = 0
      La = La .or. Lw
!      print *, "La:";  call prlog(M, N, La)
!  Separating the calculation of Vw from following SPREAD made it work.
!  Without that, it spread on axis 2.
      Vw = V < -M
!      print *, "Vw(0:N):";  print *, Vw
      Lw = SPREAD(Vw, 1, M+1)	                ! Non-basic artificials
!      print *, "Lw:";  call prlog(M, N, Lw)
      La = La .or. Lw
!      print *, "La:";  call prlog(M, N, La)
      where ( La ) A = 0.0			! Eliminate from problem
      A(0,:) = 0.0				! To form problem objective
      where ( V >= 0 ) A(0,:) = C2(V)		! Non-basic problem variables
      C1 = 0
      where ( U > 0 ) C1 = C2(U)		! Cost coeff for basic vars
      C1(0) = 1.0				! So we can include A(0,:)
      A(0,:) = SUM(A*SPREAD(C1, 2, N+1), 1)	! Objective in current basis
!      print *, "Tableau for phase 2:"
!      call prflt(M, N, A)
      call simplex(M, N, A, U, V, Iout, Ic2)	! Attempt optimization
 99   continue					! For exit
      where (V < -M) V = M + V			! Adjust codes for artificial
      where (U < -M) U = M + U			! variables from equations
      end subroutine lpsolve
