      program fullsimplex
c   
c*********************************************************************
c
c   written by Yong Li, SCCS/NPAC at Syracuse University
c   version dated: 9/24/1991
c
c   The simplex method for solving linear programming problems.  In
c   this code, the coefficient matrix is given in a matrix format.
c   There is data header file, called simplex.h for this code.
c
c   This code and data header can be found in SCCS/NPSC test suite.
c
c*********************************************************************
c 
      include 'simplex.h'
      integer numGE
      real tarray(2), etime, elapsed, elapsed1
      logical phaseII

      phaseII = .true.
      print *, ' FULL SIMPLEX METHOD FOR SOLVING LINEAR PROGRAMMING '
      print *, ' ON SEQUENTIAL MACHINE (SUN SPARC WORKSTATION) '
      print *, ' This program finds the minmum, if it exists.'
      print *, ' '

      call initial(numGE)
      if (numGE .EQ. 0) then
         print *, ' There is no phase I'
        else
         print *, ' Phase I is working...'
         elapsed = etime(tarray)
         call simplex (1)
         print *, ' Phase I takes ', iteration, ' iterations'
         if ((dabs(x(mm)) .GT. smallvalue) .or. cont) then
            phaseII = .false.
            print *, ' Phase I was not successful'
	    print *, ' The sum of artificals is ', x(mm)
          else
            print *, ' Phase I is successful'
         end if
      end if
      if (phaseII) then
         print *, ' Phase II is working ...'
         call simplex (2)
         elapsed1 = etime(tarray) - elapsed
         print *, ' Total takes ', elapsed1, ' seconds'
         print *, ' Total are ', iteration, ' iterations'
         if (bounded) then
            print *, ' Phase II is successful '
            print *, ' Minmum is Z = ', -x(m1)
         else
            print *, ' Variable ', s, ' is unbounded'
         end if
      end if

      stop
      end

      subroutine initial(numGE)

      include 'simplex.h'
      integer i, j, numG, numE, numL, numGE, numGL
      real*8 sum

      call GetData(numG, numE, numL)
      print *, ' '
      print *, m, ' constrants'
      print *, n, ' variables'
      print *, numG, ' greater or equal constrants'
      print *, numE, ' equality constrants'
      print *, numL, ' less or equal constrants'
      print *, ' ' 
      numGE = numG + numE
      numGL = numG + numL
      n1 = n + numGE + numGL
      n2 = n + numGL
      m1 = m + 1
      mm = m + 2
      nn = n1 + 1
      iteration = 0

      do i = 1, numG
         a(i, n+i) = -one
      end do
      do i = 1, numL
         a(numGE+i, n+numG+i) = one
      end do
      do i = 1, numGE
         a(i, n2+i) = one
      end do
      do j = 1, numGE
         basic(j) = n2 + j 
      end do
      do j = 1, numL
         basic(numGE+j) = n + numG + j
      end do 
      do j = 1, n1
         nonbasic(j) = .true.
      end do
      do i = 1, m
         nonbasic(basic(i)) = .false.
      end do
      do j = 1, n2
         sum = zero
         do i = 1, numGE
            sum = sum + a(i, j)
         end do
         a(mm, j) = -sum
      end do

      sum = zero
      do i = 1, numGE
         sum = sum + x(i)
      end do
      x(mm) = -sum
 
      return 
      end

      subroutine simplex (phase)

      include 'simplex.h'
      integer phase

      if (phase .EQ. 1) then
        rin = mm
        cin = n1
       else
        rin = m1
        cin = n2
      end if

      call next
      do while (cont .AND. bounded)
        call next
      end do

      return
      end

      subroutine next

      include 'simplex.h'
      integer i, j
      real*8 min, pivot, tmp, col(d1)

c
c  determining the variable to enter the basis
c 
      min = largevalue
      do j = 1, cin
        if (nonbasic(j)) then
          if (a(rin, j) .LT. min) then
            min = a(rin, j)
            s = j
          end if
        end if
      end do
c  
c  x(s) will enter the basis
c
      cont = (min .LT. -smallvalue)
c
c  check if the problem is bounded
c
      if (cont) then
c
c  check feasibility
c
        do i = 1, rin
          col(i) = a(i, s)
        end do
        i = 1
        bounded = (col(i) .GT. smallvalue)
        do while ((.NOT. bounded) .AND. (i .LT. m))
          i = i + 1
          bounded = (col(i) .GT. smallvalue)
        end do 
        if (bounded) then
          min = largevalue
          do i = 1, m
            if (col(i) .GT. smallvalue) then
              tmp = x(i)/col(i)
              if (tmp .LT. min) then
                min = tmp
                r = i
              end if
            end if
          end do
c
c  xxs[basic(r)] will leave the basis
c
          nonbasic(basic(r)) = .true.
          basic(r) = s
          nonbasic(s) = .false.
c
c  update the solution
c
          pivot = a(r, s)
          x(r) = min
          do j = 1, cin 
             a(r, j) = a(r, j)/pivot
          end do
          do i = 1, rin
            if (i .NE. r) then
    	      x(i) = x(i) - min*col(i)
              do j = 1, cin 
                a(i, j) = a(i, j) - col(i)*a(r, j)
              end do
            end if
          end do
          iteration = iteration + 1
        end if
      end if

      return
      end

C
C  Version Dated:  January 24, 1991
C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			GetData					C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	subroutine GetData(numG,numE,numL)

        include 'simplex.h'
        integer i, j, k, ii
	integer numRow, colCount, numG, numE, numL
	character*8 rowLab(d3)
	integer SearchRow, count
	character*4 tmpConstant
	character*8 nameStr, tmp1Lab, tmp3Lab, tmp4Lab, tmpLab
	character*1 rowType(d3)
	real*8 tmp1Data, tmp2Data
	
	READ(*,'(A4,T15,A8)') tmpConstant, nameStr
	if (tmpConstant .NE. 'NAME') then
	    print *, 'Input format error, missing name field'
	    stop
	endif
        print *, ' The test data file is ', nameStr
	READ(*,'(A4)') tmpConstant
	if (tmpConstant .NE. 'ROWS') then
	    print *, 'Input format error, missing ROWS field'
	    stop
	endif
	do i = 1, d1 
	    READ(*,'(A4,T2,A1,T5,A8)') tmpConstant, rowType(i), rowLab(i)
	    if (tmpConstant .EQ. 'COLU') then
		numRow = i - 1
	        goto 10
	    endif
	end do
	if (tmpConstant .NE. 'COLU') then
	    print *, 'Error: too many rows'
	    stop
	endif
10      continue
	m = numRow -1
	numG = 0
	numE = 0
	numL = 0
C
C  Looking for the row with row type N, and exchange it with the last row
C
        do i = 1, numRow
           if (rowType(i) .EQ. 'N') then
              call Exchange(rowType, rowLab, i, numRow)
              goto 20
           end if
        end do
20      continue
C
C  Rearrange the rows according to the order of row type G, E, L
C
	i = 1
	j = numRow-1
	do ii = 1, numRow-1
	   do k = i, numRow-1
	      if (rowType(k) .NE. 'G') then
	         i = k
                 goto 30
	      endif
	   end do
30        continue
	  do k = j, 1, -1
	     if (rowType(k) .EQ. 'G') then
		j = k
		goto 40
	     endif
	  end do
40        continue
	  if (i .GT. j) then
             goto 50
	  else 
	     call Exchange(rowType,rowLab,i,j)
	  endif
	end do
50      continue

	i = 1
	j = numRow-1
	do ii = 1, numRow-1
	  do k = i, numRow-1
	     if ((rowType(k) .NE. 'E') .AND. (rowType(k) .NE. 'G')) then
		i = k
		goto 60
	     endif
	  end do
60        continue
	  do k = j, 1, -1
	     if (rowType(k) .NE. 'L') then
		j = k
		goto 70
	     endif
	  end do
70        continue
	  if (i .GT. j) then
             goto 80
	  else 
	     call Exchange(rowType,rowLab,i,j)
	  endif
	end do
80      continue
	do k = 1, numRow-1
	  if (rowType(k) .EQ. 'G') then
	    numG = numG + 1
	  else
	    if (rowType(k) .EQ. 'E') then
	      numE = numE + 1
	    else
	      if (rowType(k) .EQ. 'L') then
		numL = numL + 1
	      endif
	    endif
	  endif
	end do

        do i = 1, numRow+1 
           x(i) = 0.0
           do j = 1, d2 
              a(i, j) = 0.0
           end do
        end do
        Count = 0
        colCount = 1
        do i = 1, d3 
           READ(*, '(A3,T5,A8,T15,A8,T25,F12.3,T40,A8,T50,F12.3)')
     &            tmpConstant, tmp1Lab, tmp3Lab, tmp1Data,
     &            tmp4Lab, tmp2Data
 
           if (Count .EQ. 0) tmpLab = tmp1Lab
           Count = Count + 1
           if (tmpConstant .EQ. 'RHS') then
              goto 90
            else
              if (tmp1Lab .NE. tmpLab) then
                 tmpLab = tmp1Lab
                 colCount = colCount + 1
              end if
              k = SearchRow(tmp3Lab, rowLab, numRow)
              a(k, colCount) = tmp1Data
              if (tmp4Lab .NE. '        ') then
                 k = SearchRow(tmp4Lab, rowLab, numRow)
                 a(k, colCount) = tmp2Data
              end if
           endif
        end do

        if (tmpConstant .NE. 'RHS') then
           print *, 'Error:  too many columns'
           stop
        end if

90      continue
        n = colCount
        do i = 1, d3 
           READ(*, '(A4,T5,A8,T15,A8,T25,F12.3,T40,A8,T50,F12.3)')
     &            tmpConstant, tmp1Lab, tmp3Lab, tmp1Data,
     &            tmp4Lab, tmp2Data

           if (tmpConstant .EQ. 'ENDA') then
              goto 100
            else
              k = SearchRow(tmp3Lab, rowLab, numRow)
              x(k) = tmp1Data
              if (tmp4Lab .NE. '        ') then
                 k = SearchRow(tmp4Lab, rowLab, numRow)
                 x(k) = tmp2Data
              end if
           endif
        end do

        if (tmpConstant .NE. 'ENDA') then
           print *, 'Error:  too many RHS'
           stop
        end if

100     continue
C	end of GetData

	return
	end

	subroutine Exchange(rowType,rowLab,i,j)

        include 'simplex.h'
        integer i, j
	character*8 rowLab(d1), tmpLab
	character*1 rowType(d1), tmpType

	tmpLab = rowLab(i)
	tmpType = rowType(i)
	rowLab(i) = rowLab(j)
	rowType(i) = rowType(j)
	rowLab(j) = tmpLab
	rowType(j) = tmpType
	return
	end

	integer function SearchRow(label,rowLab,num)

	implicit none
	integer i, num
	character*8 label, rowLab(num)

	do i = 1, num
	    if (label .EQ. rowLab(i)) then
		SearchRow = i
		return
	    endif
	end do
	print *, 'ERROR: unmatched row lable --- ', label
	stop
	end


