      program revised
c
c*********************************************************************
c
c written by Yong Li, SCCS/NPAC at Syracuse University
c version dated: 9/29/1991
c version dated: 10/28/1991
c version dated: 2/8/92 for MasPar, DECmpp 12000
c
c The revised simplex method for solving linear programming problems.
c In this code, the coefficient matrix A is given in a matrix format.
c This code also has a data header file which is called revised.h
c
c This code and data header can be found in the SCCS/NPAC test suite.
c
c********************************************************************
c
      include 'revised.h'
      logical phaseII

      phaseII = .true.
      print *, ' REVISED SIMPLEX METHOD FOR SOLVING LINEAR PROGRAMMING '
      print *, '              ON THE DECmpp 12000 (MasPar)             '
      print *, ' This program finds the minimal if it exists'
      print *, '  '

      call initial(a, b, xx, basic, nonbasic, iteration)
c
c  set cm timer
c
c      call cm_timer_clear(0)
c      call cm_timer_start(0)
 
      if (numGE .EQ. 0) then
         print *, ' There is no phase I'
        else
         print *, ' Phase I is working...'
         call simplex (a, b, xx, basic, nonbasic, iteration, 
     &                 cont, bounded, s, 1)
         print *, ' Phase I takes ', iteration, ' iteration'
         if ((ABS(xx(mm)) .GT. smallvalue) .OR. cont) then
            phaseII = .false.
            print *, ' Phase I not successful'
            print *, ' The sum of artifical variables is ', xx(mm)
          else
            print *, ' Phase I successful'
         end if
      end if
      if (phaseII) then
         print *, ' Phase II is working...'
         call simplex (a, b, xx, basic, nonbasic, iteration, 
     &                 cont, bounded, s, 2)
c
c  stop timer for phase II
c
c         call cm_timer_stop(0)
         if (bounded) then
            print *, ' Phase II successful '
            print *, ' Minmum is Z = ', -xx(m1)
            print *, ' Total are ', iteration, ' iteration'
         else
            print *, ' Variable ', s, ' is unbounded'
         end if
c         call cm_timer_print(0)
      end if

      stop
      end

      subroutine initial(a, b, xx, basic, nonbasic, iteration)

      include 'revised.h'
      integer i, j
c 
c get data
c
      call getdata(a, xx)
      print *, ' ' 
      print *, m, ' constrants'
      print *, n, ' variables'
      print *, numG, ' greater or equal constrants'
      print *, numE, ' equality constrants'
      print *, numL, ' less or equal constrants'
      print *, ' '

      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
      nonbasic(1:n1) = .true.
      nonbasic(basic(1:m)) = .false.
      a(mm, n2+1:n1) = one 
      xx(mm) = -SUM(xx(1:numGE))
      b(1:mm, 1:m) = zero
      do i = 1, m
        b(i, i) = one
      end do
      b(mm, 1:numGE) = -one
      iteration = 0

      return 
      end

      subroutine simplex (a, b, xx, basic, nonbasic, iteration, 
     &                    cont, bounded, s, phase)

      include 'revised.h'
      integer rin, cin, iloc(1), jloc(1), phase, i
      real*8 c(n1), d(mm), row(m), pivot, tmp, tmpA(m, n1),
     &       tmpD(mm, m), tmpC(m)
cCMF$  LAYOUT c(), d(), row()
cCMF$  ALIGN tmpA(I, J) with A(I, J)
cCMF$  ALIGN tmpD(I, J) with b(I, J)

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

      cont = .true.
      bounded = .true.
c
c  routine next
c
      do while (cont .AND. bounded)
c
c***  c = b(rin, :) X a(1:m, 1:cin)
c
        tmpA = SPREAD(b(rin,1:m),DIM=2,NCOPIES=cin)*a(1:m,1:cin)
        c(1:cin) = SUM(tmpA, DIM=1)
        c(1:cin) = c(1:cin) + a(rin, 1:cin)
        jloc = MINLOC(c(1:cin), nonbasic(1:cin))
        s = jloc(1)
        cont = (c(s) .LT. -smallvalue)

        if (cont) then
c
c***  d = b X a(1:m, s)
c
          tmpD = b*SPREAD(a(1:m, s), DIM=1, NCOPIES=mm)
          d = SUM(tmpD, DIM=2)
          bounded = MAXVAL(d(1:m)) .GT. smallvalue 
          if (bounded) then
            WHERE (d(1:m) .GT. smallvalue) tmpC(1:m) = xx(1:m)/d(1:m)
            iloc = MINLOC(tmpC(1:m), d(1:m) .GT. smallvalue)
            r = iloc(1)
            nonbasic(basic(r)) = .true.
            nonbasic(s) = .false.
            basic(r) = s
c
c  routine trans
c
            d(m1:mm) = c(s)
            if (phase .EQ. 1) then
              d(m1) = a(m1, s) + SUM(b(m1, :)*a(1:m, s))
            end if
            pivot = d(r)
            tmp = xx(r)/pivot
            xx = xx - tmp*d
            xx(r) = tmp
            b(r, :) = b(r, :)/pivot
            row = b(r, :)
            b = b - SPREAD(d, DIM=2, NCOPIES=m) *
     *              SPREAD(b(r, :), DIM=1, NCOPIES=mm)
            b(r, :) = row
            iteration = iteration + 1
          end if       
        end if
      end do

      return
      end

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

	subroutine GetData(a, xx)

      include 'revised.h'
      integer i, j, k, ii, SearchRow, count, numRow, colCount
      character*1 rowType(d1)
      character*4 tmpConstant
      character*8 rowlab(d3), nameStr, tmp1Lab, tmp3Lab, tmp4Lab, tmpLab
      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
c
	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
c  initialize the solution vector and coefficient matrix A
c
        xx = zero
        a = zero
c
        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
c
        do i = 1, d2 
           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)
              xx(k) = tmp1Data
              if (tmp4Lab .NE. '        ') then
                 k = SearchRow(tmp4Lab, rowLab, numRow)
                 xx(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 'revised.h'

	integer i, j
        character*8 tmpLab, rowLab(d3)
        character*1 tmpType, rowType(d1)

	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


