      program cm_fullsimplex

      include 'simplex.h'
      logical phaseII

      real*8 a_out
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 CONNECTION MACHINE'
      print *, ' '

      call initial (a, nonbasic)
c
c  set cm timer
c
      call cm_timer_clear(0)
      call cm_timer_start(0)

      if (numGE .EQ. 0) then
        print *, ' There is no phase I'
       else
        print *, ' Phase I is working...'
        call simplex (1, a, nonbasic)
        print *, ' Phase I takes', iteration, ' iterations'
        a_out = a(mm,nn)
        if ((DABS(a_out) .GT. small) .OR. cont) then
          phaseII = .false.
          print *, ' Phase I not successful'
          print *, ' Sum of the artifical variables is ', a_out
         else
          print *, ' Phase I successful'
        end if
      end if
      if (phaseII) then
        print *, ' Phase II is working...'
        call simplex (2, a, nonbasic)
c
c  stop timer
c
        call cm_timer_stop(0)

        print *, ' Phase II takes  ', iteration, ' iterations'
        if (bounded) then
          print *, ' Phase II successful '
          a_out = a(m1, nn)
          print *, ' Minmum is Z = ', -a_out
         else
          print *, ' Variable ', s, ' is unbounded'
        end if
        call cm_timer_print(0)
      end if


      end

      subroutine initial (a, nonbasic)

      include 'simplex.h'
      integer i, j, k
      real*8 asum

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

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

c     FORALL (j=1:n2) a(mm, j) = -SUM(a(1:numGE, j))
!hpf$ independent, local_access
      do j = 1, n2 
         a(mm,j) = 0
         do k = 1, numGE
            a(mm,j) = a(mm,j) - a(k,j)
         end do
      end do

      asum = -SUM(a(1:numGE, nn))
      a(mm, nn) = asum
      iteration = 0
 
      end

      subroutine simplex (phase, a, nonbasic)

      include 'simplex.h'
      integer phase

      if (phase .EQ. 1) then
        rindex = mm
        cindex = n1
       else
        rindex = m1
        cindex = n2
      end if

      call next (a, nonbasic)
      do while (cont .AND. bounded)
        call next (a, nonbasic)
      end do

      end

      subroutine next (a, nonbasic)

      include 'simplex.h'
      real*8 pivot, tmp
      integer iloc, jloc, k, k1
      real*8 mina
      real*8 repa (m)
cmf$  layout repa (:replicated)

c     jloc = MINLOC(a(rindex, 1:cindex), nonbasic(1:cindex))
      jloc = 1
      mina = 100000.0
!hpf$ independent, local_access
      do k = 1, cindex 
         if (nonbasic(k)) then
            reduce (minval, mina, a(rindex,k), jloc, k)
         end if
      end do

      s = jloc
      cont = (mina .LT. -small)
      if (cont) then
        bounded = ANY(a(1:m, s) .GT. small)
        if (bounded) then
c         iloc = MINLOC(a(1:m, nn)/a(1:m, s), a(1:m, s) .GT. small)
          iloc = 1
          mina = 100000.0
          repa (1:m) = a(1:m,s)
!hpf$ independent, local_access
          do k = 1, m 
             if (repa(k) .GT. small) then
                reduce (minval, mina, a(k,nn)/repa(k), iloc, k)
             end if
          end do
          r = iloc
          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:cindex) = a(r, 1:cindex)/pivot
          row(1:cindex) = a(r, 1:cindex)
          forall (k=1:cindex,k1=1:mm)
            tmpA (k1,k) = a(k1,s) * a(r,k)
          end forall
c         tmpA = SPREAD(a(:, s), 2, cindex) *
c    &           SPREAD(a(r, 1:cindex), 1, mm)
          a(:, 1:cindex) = a(:, 1:cindex) - tmpA(:, 1:cindex)
          a(r, 1:cindex) = row(1:cindex)
          iteration = iteration + 1
        end if
      end if

      end

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

	subroutine GetData (a, nonbasic)

        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

	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
	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
		goto 100
	    endif
	end do
	print *, 'ERROR: unmatched row lable --- ', label
  100   continue
	end


