home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / adaptor.zip / adapt.zip / adaptor / examples / hpf_gen / simplex / simplex.f < prev    next >
Text File  |  1993-06-25  |  9KB  |  388 lines
  1.       program cm_fullsimplex
  2.  
  3.       include 'simplex.h'
  4.       logical phaseII
  5.  
  6.       real*8 a_out
  7. c
  8. c********************************************************************
  9. c
  10. c   written by Yong Li, SCCS/NPAC at Syracuse University
  11. c   version date: 9/29/1991
  12. c
  13. c   The simplex method for solving linear programming problems.
  14. c   There is a data header file called simplex.h
  15. c
  16. c*******************************************************************
  17. c
  18.  
  19.       phaseII = .true.
  20.       print *, ' '
  21.       print *, ' FULL SIMPLEX METHOD FOR LINEAR PROGRAMMING '
  22.       print *, '      ON THE CONNECTION MACHINE'
  23.       print *, ' '
  24.  
  25.       call initial (a, nonbasic)
  26. c
  27. c  set cm timer
  28. c
  29.       call cm_timer_clear(0)
  30.       call cm_timer_start(0)
  31.  
  32.       if (numGE .EQ. 0) then
  33.         print *, ' There is no phase I'
  34.        else
  35.         print *, ' Phase I is working...'
  36.         call simplex (1, a, nonbasic)
  37.         print *, ' Phase I takes', iteration, ' iterations'
  38.         a_out = a(mm,nn)
  39.         if ((DABS(a_out) .GT. small) .OR. cont) then
  40.           phaseII = .false.
  41.           print *, ' Phase I not successful'
  42.           print *, ' Sum of the artifical variables is ', a_out
  43.          else
  44.           print *, ' Phase I successful'
  45.         end if
  46.       end if
  47.       if (phaseII) then
  48.         print *, ' Phase II is working...'
  49.         call simplex (2, a, nonbasic)
  50. c
  51. c  stop timer
  52. c
  53.         call cm_timer_stop(0)
  54.  
  55.         print *, ' Phase II takes  ', iteration, ' iterations'
  56.         if (bounded) then
  57.           print *, ' Phase II successful '
  58.           a_out = a(m1, nn)
  59.           print *, ' Minmum is Z = ', -a_out
  60.          else
  61.           print *, ' Variable ', s, ' is unbounded'
  62.         end if
  63.         call cm_timer_print(0)
  64.       end if
  65.  
  66.  
  67.       end
  68.  
  69.       subroutine initial (a, nonbasic)
  70.  
  71.       include 'simplex.h'
  72.       integer i, j, k
  73.       real*8 asum
  74.  
  75. c
  76. c  get data
  77. c
  78.       print *, ' '
  79.       call getdata (a, nonbasic)
  80.       print *, ' '
  81.       print *, m, ' constrants'
  82.       print *, n, ' variables'
  83.       print *, ' constrants (G  E  L) = (', numG, numE, numL, ' )'
  84.       print *, ' '
  85.  
  86.       FORALL (i=1:mm, j=n+1:nn-1) a(i, j) = zero
  87.       FORALL (i=n+1:n+numG) a(i-n,i) = -one      
  88.       FORALL (i=n+numG+1:n+numG+numL) a(numGE+i-n-numG, i) = one
  89.       FORALL (i=n2+1:n2+numGE) a(i-n2, i) = one
  90.       FORALL (j=1:numGE) basic(j) = n2 + j
  91.       FORALL (j=1:numL) basic(numGE+j) = n + numG + j
  92.       nonbasic(1:n1) = .true.
  93.       do j = 1, m
  94.          nonbasic(basic(j)) = .false.
  95.       end do
  96.  
  97. c     FORALL (j=1:n2) a(mm, j) = -SUM(a(1:numGE, j))
  98. !hpf$ independent, local_access
  99.       do j = 1, n2 
  100.          a(mm,j) = 0
  101.          do k = 1, numGE
  102.             a(mm,j) = a(mm,j) - a(k,j)
  103.          end do
  104.       end do
  105.  
  106.       asum = -SUM(a(1:numGE, nn))
  107.       a(mm, nn) = asum
  108.       iteration = 0
  109.  
  110.       end
  111.  
  112.       subroutine simplex (phase, a, nonbasic)
  113.  
  114.       include 'simplex.h'
  115.       integer phase
  116.  
  117.       if (phase .EQ. 1) then
  118.         rindex = mm
  119.         cindex = n1
  120.        else
  121.         rindex = m1
  122.         cindex = n2
  123.       end if
  124.  
  125.       call next (a, nonbasic)
  126.       do while (cont .AND. bounded)
  127.         call next (a, nonbasic)
  128.       end do
  129.  
  130.       end
  131.  
  132.       subroutine next (a, nonbasic)
  133.  
  134.       include 'simplex.h'
  135.       real*8 pivot, tmp
  136.       integer iloc, jloc, k, k1
  137.       real*8 mina
  138.       real*8 repa (m)
  139. cmf$  layout repa (:replicated)
  140.  
  141. c     jloc = MINLOC(a(rindex, 1:cindex), nonbasic(1:cindex))
  142.       jloc = 1
  143.       mina = 100000.0
  144. !hpf$ independent, local_access
  145.       do k = 1, cindex 
  146.          if (nonbasic(k)) then
  147.             reduce (minval, mina, a(rindex,k), jloc, k)
  148.          end if
  149.       end do
  150.  
  151.       s = jloc
  152.       cont = (mina .LT. -small)
  153.       if (cont) then
  154.         bounded = ANY(a(1:m, s) .GT. small)
  155.         if (bounded) then
  156. c         iloc = MINLOC(a(1:m, nn)/a(1:m, s), a(1:m, s) .GT. small)
  157.           iloc = 1
  158.           mina = 100000.0
  159.           repa (1:m) = a(1:m,s)
  160. !hpf$ independent, local_access
  161.           do k = 1, m 
  162.              if (repa(k) .GT. small) then
  163.                 reduce (minval, mina, a(k,nn)/repa(k), iloc, k)
  164.              end if
  165.           end do
  166.           r = iloc
  167.           nonbasic(basic(r)) = .true.
  168.           basic(r) = s
  169.           nonbasic(s) = .false.
  170. c
  171. c  update the solution
  172. c
  173.           pivot = a(r, s)
  174.           tmp = a(r, nn)/pivot
  175.           a(:, nn) = a(:, nn) - tmp*a(:, s)
  176.           a(r, nn) = tmp
  177.           a(r, 1:cindex) = a(r, 1:cindex)/pivot
  178.           row(1:cindex) = a(r, 1:cindex)
  179.           forall (k=1:cindex,k1=1:mm)
  180.             tmpA (k1,k) = a(k1,s) * a(r,k)
  181.           end forall
  182. c         tmpA = SPREAD(a(:, s), 2, cindex) *
  183. c    &           SPREAD(a(r, 1:cindex), 1, mm)
  184.           a(:, 1:cindex) = a(:, 1:cindex) - tmpA(:, 1:cindex)
  185.           a(r, 1:cindex) = row(1:cindex)
  186.           iteration = iteration + 1
  187.         end if
  188.       end if
  189.  
  190.       end
  191.  
  192. C
  193. C  Version Dated:  January 24, 1991
  194. C 
  195. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  196. C            GetData                    C
  197. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  198.  
  199.     subroutine GetData (a, nonbasic)
  200.  
  201.         include 'simplex.h'
  202.         integer i, j, k, ii
  203.     integer numRow, colCount
  204.     character*8 rowLab(d3)
  205.     integer SearchRow, count
  206.     character*4 tmpConstant
  207.     character*8 nameStr, tmp1Lab, tmp3Lab, tmp4Lab, tmpLab
  208.     character*1 rowType(d3)
  209.     real*8 tmp1Data, tmp2Data
  210.     
  211.     READ(*,'(A4,T15,A8)') tmpConstant, nameStr
  212.     if (tmpConstant .NE. 'NAME') then
  213.         print *, 'Input format error, missing name field'
  214.         stop
  215.     endif
  216.         print *, ' The test data file is ', nameStr
  217.     READ(*,'(A4)') tmpConstant
  218.     if (tmpConstant .NE. 'ROWS') then
  219.         print *, 'Input format error, missing ROWS field'
  220.         stop
  221.     endif
  222.     do i = 1, d1 
  223.         READ(*,'(A4,T2,A1,T5,A8)') tmpConstant, rowType(i), rowLab(i)
  224.         if (tmpConstant .EQ. 'COLU') then
  225.         numRow = i - 1
  226.             goto 10
  227.         endif
  228.     end do
  229.     if (tmpConstant .NE. 'COLU') then
  230.         print *, 'Error: too many rows'
  231.         stop
  232.     endif
  233. 10      continue
  234. C
  235. C  Looking for the row with row type N, and exchange it with the last row
  236. C
  237.         do i = 1, numRow
  238.            if (rowType(i) .EQ. 'N') then
  239.               call Exchange(rowType, rowLab, i, numRow)
  240.               goto 20
  241.            end if
  242.         end do
  243. 20      continue
  244. C
  245. C  Rearrange the rows according to the order of row type G, E, L
  246. C
  247.     i = 1
  248.     j = numRow-1
  249.     do ii = 1, numRow-1
  250.        do k = i, numRow-1
  251.           if (rowType(k) .NE. 'G') then
  252.              i = k
  253.                  goto 30
  254.           endif
  255.        end do
  256. 30        continue
  257.       do k = j, 1, -1
  258.          if (rowType(k) .EQ. 'G') then
  259.         j = k
  260.         goto 40
  261.          endif
  262.       end do
  263. 40        continue
  264.       if (i .GT. j) then
  265.              goto 50
  266.       else 
  267.          call Exchange(rowType,rowLab,i,j)
  268.       endif
  269.     end do
  270. 50      continue
  271.  
  272.     i = 1
  273.     j = numRow-1
  274.     do ii = 1, numRow-1
  275.       do k = i, numRow-1
  276.          if ((rowType(k) .NE. 'E') .AND. (rowType(k) .NE. 'G')) then
  277.         i = k
  278.         goto 60
  279.          endif
  280.       end do
  281. 60        continue
  282.       do k = j, 1, -1
  283.          if (rowType(k) .NE. 'L') then
  284.         j = k
  285.         goto 70
  286.          endif
  287.       end do
  288. 70        continue
  289.       if (i .GT. j) then
  290.              goto 80
  291.       else 
  292.          call Exchange(rowType,rowLab,i,j)
  293.       endif
  294.     end do
  295. 80      continue
  296. c
  297.         a(1:numRow+1, 1:nn) = 0.0
  298.         Count = 0
  299.         colCount = 1
  300.         do i = 1, d3 
  301.            READ(*, '(A3,T5,A8,T15,A8,T25,F12.3,T40,A8,T50,F12.3)')
  302.      &            tmpConstant, tmp1Lab, tmp3Lab, tmp1Data,
  303.      &            tmp4Lab, tmp2Data
  304.  
  305.            if (Count .EQ. 0) tmpLab = tmp1Lab
  306.            Count = Count + 1
  307.            if (tmpConstant .EQ. 'RHS') then
  308.               goto 90
  309.             else
  310.               if (tmp1Lab .NE. tmpLab) then
  311.                  tmpLab = tmp1Lab
  312.                  colCount = colCount + 1
  313.               end if
  314.               k = SearchRow(tmp3Lab, rowLab, numRow)
  315.               a(k, colCount) = tmp1Data
  316.               if (tmp4Lab .NE. '        ') then
  317.                  k = SearchRow(tmp4Lab, rowLab, numRow)
  318.                  a(k, colCount) = tmp2Data
  319.               end if
  320.            endif
  321.         end do
  322.  
  323.         if (tmpConstant .NE. 'RHS') then
  324.            print *, 'Error:  too many columns'
  325.            stop
  326.         end if
  327.  
  328. 90      continue
  329.         do i = 1, d3 
  330.            READ(*, '(A4,T5,A8,T15,A8,T25,F12.3,T40,A8,T50,F12.3)')
  331.      &            tmpConstant, tmp1Lab, tmp3Lab, tmp1Data,
  332.      &            tmp4Lab, tmp2Data
  333.  
  334.            if (tmpConstant .EQ. 'ENDA') then
  335.               goto 100
  336.             else
  337.               k = SearchRow(tmp3Lab, rowLab, numRow)
  338.               a(k, nn) = tmp1Data
  339.               if (tmp4Lab .NE. '        ') then
  340.                  k = SearchRow(tmp4Lab, rowLab, numRow)
  341.                  a(k, nn) = tmp2Data
  342.               end if
  343.            endif
  344.         end do
  345.  
  346.         if (tmpConstant .NE. 'ENDA') then
  347.            print *, 'Error:  too many RHS'
  348.            stop
  349.         end if
  350.  
  351. 100     continue
  352. C    end of GetData
  353.  
  354.     end
  355.  
  356.     subroutine Exchange(rowType,rowLab,i,j)
  357.  
  358.         include 'simplex.h'
  359.         integer i, j
  360.     character*8 rowLab(d1), tmpLab
  361.     character*1 rowType(d1), tmpType
  362.  
  363.     tmpLab = rowLab(i)
  364.     tmpType = rowType(i)
  365.     rowLab(i) = rowLab(j)
  366.     rowType(i) = rowType(j)
  367.     rowLab(j) = tmpLab
  368.     rowType(j) = tmpType
  369.     end
  370.  
  371.     integer function SearchRow(label,rowLab,num)
  372.  
  373.     implicit none
  374.     integer i, num
  375.     character*8 label, rowLab(num)
  376.  
  377.     do i = 1, num
  378.         if (label .EQ. rowLab(i)) then
  379.         SearchRow = i
  380.         goto 100
  381.         endif
  382.     end do
  383.     print *, 'ERROR: unmatched row lable --- ', label
  384.   100   continue
  385.     end
  386.  
  387.  
  388.