      program fullsimplex

      include 'simplex.h'
      logical phaseII
c
c********************************************************************
c
c   written by Yong Li, SCCS/NPAC at Syracuse University
c   version date: 9/29/1991
c
c   The simplex method for solving linear programming problems.
c   There is a data header file called simplex.h
c
c*******************************************************************
c

      phaseII = .true.
      print *, ' '
      print *, ' FULL SIMPLEX METHOD FOR LINEAR PROGRAMMING '
      print *, '      ON THE DECmpp 12000 (MasPar) '
      print *, ' '

      call initial(a, basic, nonbasic, iteration)
c
      if (numGE .EQ. 0) then
        print *, ' There is no phase I'
       else
        print *, ' Phase I is working...'
        call simplex (a,basic,nonbasic,iteration,cont,bounded,s,1)
        print *, ' Phase I takes', iteration, ' iterations'
        if ((DABS(a(mm, nn)) .GT. small) .OR. cont) then
          phaseII = .false.
          print *, ' Phase I not successful'
          print *, ' Sum of the artifical variables is ', a(mm, nn)
         else
          print *, ' Phase I successful'
        end if
      end if
      if (phaseII) then
        print *, ' Phase II is working...'
        call simplex (a,basic,nonbasic,iteration,cont,bounded,s,2)
c
        print *, ' Phase II takes  ', iteration, ' iterations'
        if (bounded) then
          print *, ' Phase II successful '
          print *, ' Minmum is Z = ', -a(m1, nn)
         else
          print *, ' Variable ', s, ' is unbounded'
        end if
      end if


      stop
      end

      subroutine initial(a, basic, nonbasic, iteration)

      include 'simplex.h'
      integer i, j

c
c  get data
c
      print *, ' '
      call getdata(a)
      print *, ' '
      print *, m, ' constrants'
      print *, n, ' variables'
      print *, ' constrants (G  E  L) = (', numG, numE, numL, ' )'
      print *, ' '

      a(1:mm, n+1:nn-1) = zero
      do i = 1, numG
        a(i, n+i) = -one
      end do
      do i = 1, numL
        a(numGE+i, n+numG+i) = one
        basic(numGE+i) = n+numG+i
      end do
      do i = 1, numGE
        a(i, n2+i) = one
        basic(i) = n2+i
      end do
      nonbasic(1:n1) = .true.
      nonbasic(basic(1:m)) = .false.
      do j = 1, n2
        a(mm, j ) = -SUM(a(1:numGE, j))
      end do
      a(mm, nn) = -SUM(a(1:numGE, nn))
      iteration = 0
 
      return 
      end

      subroutine simplex (a,basic,nonbasic,iteration,cont,
     *                    bounded,s,phase)

      include 'simplex.h'
      integer phase

      if (phase .EQ. 1) then
        rindx = mm
        cindx = n1
       else
        rindx = m1
        cindx = n2
      end if

      call next(a,basic,nonbasic,iteration,cont,bounded,
     *          rindx,cindx,s)
      do while (cont .AND. bounded)
        call next(a,basic,nonbasic,iteration,cont,bounded,
     *            rindx,cindx,s)
      end do

      return
      end

      subroutine next(a,basic,nonbasic,iteration,cont,bounded,
     *               rindx,cindx,s)

      include 'simplex.h'
      real*8 pivot, tmp, tmpC(m), d(mm)
      integer iloc(1), jloc(1), i , j

      jloc = MINLOC(a(rindx, 1:cindx), nonbasic(1:cindx))
      s = jloc(1)
      cont = (a(rindx, s) .LT. -small)
      if (cont) then
        d = a(:, s)
        bounded = (MAXVAL(d(1:m)) .GT. small)
        if (bounded) then
          WHERE(d(1:m) .GT. small) tmpC(1:m) = a(1:m, nn)/d(1:m)
          iloc = MINLOC(tmpC(1:m), d(1:m) .GT. small)
          r = iloc(1)
          nonbasic(basic(r)) = .true.
          basic(r) = s
          nonbasic(s) = .false.
c
c  update the solution
c
          pivot = a(r, s)
          tmp = a(r, nn)/pivot
          a(:, nn) = a(:, nn) - tmp*a(:, s)
          a(r, nn) = tmp
          a(r, 1:cindx) = a(r, 1:cindx)/pivot
          row(1:cindx) = a(r, 1:cindx)
          tmpA = SPREAD(d, DIM=2, NCOPIES=nn) *
     &           SPREAD(a(r, :), DIM=1, NCOPIES=mm)
          a(:, 1:cindx) = a(:, 1:cindx) - tmpA(:, 1:cindx)
          a(r, 1:cindx) = row(1:cindx)
          iteration = iteration + 1
        end if
      end if

      return
      end

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

	subroutine GetData(a)

        include 'simplex.h'
        integer i, j, k, ii
	integer numRow, colCount
	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
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
c
        a(1:numRow+1, 1:nn) = 0.0
        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
        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)
              a(k, nn) = tmp1Data
              if (tmp4Lab .NE. '        ') then
                 k = SearchRow(tmp4Lab, rowLab, numRow)
                 a(k, nn) = 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


