      program simplex_node
      implicit none

      include '/usr/include/fcube.h'
      include 'simcube.h'
      include 'simnode.h'
   
c---------------------------------------------------------------------------
c   
c  This code perfomes the simplex method to solve
c  linear programming problems on hypercube iPSC/2
c  The matrix A is distributed by using block partition
c 
c  This is a node program
c
c version date: 4/16/1991
c
c-------------------------------------------------------------------------

c  local variables
c
      integer i, j, k

      thisNode = mynode()
      hostNode = myhost()
      numNode = numnodes()
      numDim = nodeDim()
      iteration = 0

      call initial
c
c  code for node 0
c
      if (thisNode .eq. 0) then
        if (numGE .eq. 0) then
          print *, ' There is no phase I'
         else
          print *, 'w = ', x(numEqu+1)
          print *, ' Starting phase I....'
        end if
        timeBegin = mclock()
      end if
c
c  code for all node
c
      if ( numGE .ne. 0 ) then
        maxRow = numEqu + 1
        call simplex(1)
      end if
c
c  code for node 0
c
      if ((numGE .ne. 0) .and. (thisNode .eq. 0)) then
        print *, 'w = ', x(numEqu+1)
        if ( .not. bounded) then
          print *, ' It is not bounded'
        end if
        if (cont) then
          print *, ' Phase I was not successful, please killcube'
         else
          print *, ' Phase I was successful'
          print *, ' Starting phase II'
        end if
      end if
c
c  send matrix a back to hostNode
c
       do i = 1, numEqu
         do j = 1, numCol
           realBuf(j) = a(i, j)
         end do
         call csend(gtype+4+i,realBuf,
     *                numCol*doubleSize,hostNode,npid)
       end do
c
c  redistribution coefficient matrx a in phase II
c
      if (numGE .ne. 0) then
        if (numNode .gt. 1) then
          minCol = allColRow(3, thisNode+1)
          maxCol = allColRow(4, thisNode+1)
          numCol = maxCol - minCol + 1
          k = minCol
          do j = 1, numCol
            nonbasic_node(j) = nonbasic(k)
            k = k + 1
          end do
          do i = 1, numEqu
            call crecv(gtype+2,realBuf,numCol*doubleSize)
            do j = 1, numCol
              a(i, j) = realBuf(j)
            end do
          end do 
         else
          if ( maxCol .eq. numVarI) then
            maxCol = numVarII
            numCol = maxCol - minCol + 1
          end if
        end if
      end if
c
c code for all node
c
      maxRow = numEqu
      call simplex(2)
c
c  code for node 0
c
      if (thisNode .eq. 0) then
        timeEnd = mclock()
        if (bounded) then
          print *, ' Phase II was successful'
         else
          print *, ' Phase II was not successful, plsase killcube'
          print *, ' Variable', pivotCol, ' is not bounded'
        end if
        print *, ' Elapsed time: ', timeEnd-timeBegin
        print *, ' At node, Min =', x(numEqu)
        call csend(gtype,x(numEqu),doubleSize,hostNode,NPID)
        print *, ' It takes ', iteration, ' iterations'
      end if
c
c  code for node 0
c
      if (thisNode .eq. 0) then
        call csend(gtype, x(numEqu), doubleSize, hostNode,npid)
      end if

      stop
      end

c-----------------------------------------------------------------

      subroutine initial
      implicit none

      include 'simcube.h'
      include 'simnode.h'
c
c  local variables
c
      integer i, j, k

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  receiving info numVar, numEqu, numCon, numG, numE, numL
c
      call crecv(gtype, info, 6*intSize)

      numVar = info(1)
      numEqu = info(2)
      numCon = info(3)
      numG = info(4)
      numE = info(5)
      numL = info(6)
      numGL = numG + numL
      numGE = numG + numE
      numVarI = numVar + numGE + numGL
      numVarII = numVar + numGL
c
c  receiving allColRow 
c
      call crecv(gtype-1,allColRow,4*numNode*intSize)
c
c  receiving the right hand side x
c
      call crecv(gtype-2, x, (numEqu+1)*doubleSize)
c
c  receiving coefficient matrix A at node
c
      minCol = allColRow(1, thisNode+1)
      maxCol = allColRow(2, thisNode+1)
      numCol = maxCol - minCol + 1
      do i = 1, numEqu+1
        call crecv(gtype+1,realBuf,numCol*doubleSize)
        do j = 1, numCol
          a(i, j) = realBuf(j)
        end do
      end do
c
c  code for all node
c
      do j = 1, numVarI
         nonbasic_node(j) = .true.
         nonbasic(j) = .true.
      end do
      do i = numVar+numG+1, numVarII
         basic(i-numVar+numE) = i
         nonbasic(i) = .false.
         if ((i .ge. minCol) .and. (i .le. maxCol)) then
           nonbasic_node(i-minCol+1) = .false.
         end if
      end do 
      do i = numVarII+1, numVarI
        basic(i-numVarII) = i
        nonbasic(i) = .false.
         if ((i .ge. minCol) .and. (i .le. maxCol)) then
           nonbasic_node(i-minCol+1) = .false.
         end if
      end do

      return
      end

c-----------------------------------------------------------------

      subroutine simplex(phase)
      implicit none

      include 'simcube.h'
      include 'simnode.h'
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      call next(phase)
      do while (cont .and. bounded)
         call trans(phase)
         call next(phase)
      end do

      return
      end

c------------------------------------------------------------------
      subroutine next(phase)
      implicit none
 
      include 'simcube.h'
      include 'simnode.h'
c
c  local variables
c
      double precision large, small
      parameter (large = 1.0d+7)
      parameter (small = 1.0d-7)

      double precision val, sum
      integer i, j, k
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   Find the leaving column
c
c  code for all nodes
c
      pivotCol = -1
      maxval = -large
      do j = 1, numCol
        if (nonbasic_node(j)) then
          if ((a(maxRow, j) .gt. maxval) .or.
     *                        (pivotCol .eq. -1)) then
            maxval = a(maxRow, j)
            pivotCol = j
          end if
        end if 
      end do
c
c  reduct on max over all nodes to get one pivotCol and maxval
c
      realBuf(1) = maxval
      realBuf(2) = minCol + pivotCol -1
      call GOPF(realBuf, 2*doubleSize, work2, MaxReduct)
      maxval = realBuf(1)
      pivotCol = realBuf(2)
      cont = (maxval .gt. small)
c
c  pivot column is in this node
c
      kIndxCol = (pivotCol .le. maxCol) .and. (pivotCol .ge. minCol)

      if (cont) then
        if (kIndxCol) then
c
c  check if this column is bounded
c
          bounded = .false.
          k = pivotCol - minCol + 1
          do i = 1, numCon
            if ((.not. bounded) .and. (a(i, k) .gt. small)) then
              bounded = .true.
            end if
          end do
        end if 
        call GLOR(bounded, 1, worklogical)
c
c  find the pivot row
c
        minval = large
        if (bounded) then
          if (kIndxCol) then
            do i = 1, numCon
              if (a(i, k) .gt. small) then
                val = x(i)/a(i, k)
                if (val .lt. minval) then
                  minval = val
                  pivotRow = i
                end if
              end if
            end do
          end if
          if (kIndxCol) then
            nonbasic_node(k) = .false.
c
c  leaving column of matrix a is in this node, this column needs
c  to send to all nodes, k is defined as above
c
            do i = 1, maxRow
              realBuf(i) = a(i, k)
            end do
            realBuf(maxRow+1) = minval
            realBuf(maxRow+2) = pivotRow
            if (numNode .gt. 1) then
              call csend(stype+3*iteration+1,realBuf,
     *                  (maxRow+2)*doubleSize,allNode,npid)
            end if
           else if (pivotCol .ne. -1) then
            call crecv(stype+3*iteration+1,realBuf,
     *                (maxRow+2)*doubleSize)
            pivotRow = realBuf(maxRow+2)
          end if
c
c   code for a special node
c
          if ((basic(pivotRow) .ge. minCol) .and.
     *        (basic(pivotRow) .le. maxCol)) then
            nonbasic_node(basic(pivotRow)-minCol+1) = .true.
          end if
          nonbasic(basic(pivotRow)) = .true.
c
c  code for all node
c
          basic(pivotRow) = pivotCol
          nonbasic(pivotCol) = .false.
        end if       
      end if

      return
      end

c-------------------------------------------------------------------

      subroutine trans(phase)
      implicit none

      include 'simcube.h'
      include 'simnode.h'
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c
c  local variables
c
      integer i, j
c
c  code for all nodes
c
      minval = realBuf(maxRow+1)
      pivot = realBuf(pivotRow)
      do j = 1, numCol 
        a(pivotRow, j) = a(pivotRow, j)/pivot
      end do
c
c  code for all nodes
c
      do i = 1, maxRow
        if (i .ne. pivotRow) then
          x(i) = x(i) - realBuf(i)*minval
          do j = 1, numCol
            a(i, j) = a(i, j) - realBuf(i)*a(pivotRow, j)
          end do
        end if
      end do
      x(pivotRow) = minval
      iteration = iteration + 1

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			MaxReduct				C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        subroutine MaxReduct(mine, work)
        implicit none

        double precision mine(2), work(2)

        if ((mine(2) .LT. 0.0) .or. ((work(2) .gt. 0.0) .and. 
     &	     ((mine(1) .LT. work(1)) .or. ((mine(1) .eq. work(1))
     &      .and. (mine(2) .GT. work(2)))))) then
          mine(1) = work(1)
          mine(2) = work(2)
        end if
        return
        end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C			MinReduct				C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        subroutine MinReduct(mine, work)
        implicit none

        double precision mine(2), work(2)

        if ((mine(2) .LT. 0.0) .or. ((work(2) .gt. 0.0) .and. 
     &	    ((mine(1) .GT. work(1)) .OR. ((mine(1) .eq. work(1))
     &      .AND. (mine(2) .GT. work(2)))))) then
          mine(1) = work(1)
          mine(2) = work(2)
        end if
        return
        end

