program main

c*********************************************************************72
c
c  remarks
c  since this program is complete in all respects, it can be
c  run as it is without any additional modifications or
c  instructions.  in such case, follow the input format as given.
c
c  program for solving linear and quadratic programming
c  problems in the form w = m * z + q, w * z = 0, w and z nonnegative
c  by lemke's algorithm.
c
c  main program which calls the six subroutines - matrix,
c  initia, newbas, sort, pivot and pprint in proper order.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer ip
      integer ir
      integer l1
      integer mbasis(100)
      integer n 
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      integer no
      real q(50)
      real w(50)
      real z(50)

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z
c
c  description of parameters in common.
c
c    am     a two dimensional array containing the
c           elements of matrix m.
c
c    q      a singly subscripted array containing the
c           elements of vector q.
c
c    l1     an integer variable indicating the number of
c           iterations taken for each problem.
c
c    a      a two dimensional array containing the
c           elements of the inverse of the current basis.
c
c    w      a singly subscripted array containing the values
c           of w variables in each solution.
c
c    z      a singly subscripted array containing the values
c           of z variables in each solution.
c
c    nl1    an integer variable taking value 1 or 2 depend-
c           ing on whether variable w or z leaves the basis.
c
c    ne1    similar to nl1 but indicates variable entering.
c
c    nl2    an integer variable indicating what component
c           of w or z variable leaves the basis.
c
c    ne2    similar to nl2 but indicates variable entering.
c
c    a      a singly subscripted array containing the
c           elements of the transformed column that is
c           entering the basis.
c
c    ir     an integer variable denoting the pivot row at
c           each iteration.  also used to indicate termina-
c           tion of a problem by giving it a value of 1000.
c
c    mbasis a single subscripted array-indicator for the
c           basis variables.  two indicators are used for
c           each basic variable - one indicating whether
c           it is a w or z and another indicating what
c           component of w or z.
c
c  read in the value of variable ip indicating the
c  number of problems to be solved.
c
      read ( 5, 3 ) ip
c
c  variable no indicates the current problem being solved.
c
      no = 0
1     no = no + 1
      if ( no .gt. ip ) go to 5
      write ( 6, 2 ) no
2     format ( 1h1, 10x, 11hproblem no., i2 )
c
c  read in the size of the matrix m.
c
      read ( 5, 3 ) n
3     format ( i2 )
c
c  program calling sequence.
c
      call matrix ( n )
c
c  parameter n indicates the problem size.
c
      call initia ( n )
c
c  since for any problem termination can occur in initia,
c  newbas or sort subroutine, the value of ir is matched with
c  1000 to check whether to continue or to to next problem.
c
      if ( ir .eq. 1000 ) go to 1
4     call newbas ( n )
      if ( ir .eq. 1000 ) go to 1
      call sort ( n )
      if ( ir .eq. 1000 ) go to 1
      call pivot ( n )
      go to 4
5     stop
      end
      subroutine matrix ( n )

c*********************************************************************72
c
cc matrix() initializes and reads in the various input data.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer i
      integer ir
      integer j
      integer l1
      integer mbasis(100)
      integer n
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      real q(50)
      real w(50)
      real z(50)

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z
c
c  read the elements of m matrix column by column.
c
      do j = 1, n
        read ( 5, 2 ) ( am(i,j), i = 1, n )
      end do
2     format ( 7f10.5 )
c
c  read the elements of q vector.
c
      read ( 5, 2 ) ( q(i), i = 1, n )
c
c  in iteration 1, basis inverse is an identity matrix.
c
      do j = 1, n
        do i = 1, n
          if ( i .ne. j ) then
            b(i,j) = 0.0e+00
          else
            b(i,j) = 1.0e+00
          end if
        end do
      end do

      return
      end
      subroutine initia ( n )

c*********************************************************************72
c
c  purpose - to find the initial almost complementary solution
c  by adding an artificial variable z0.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer i
      integer ir
      integer j
      integer l
      integer l1
      integer mbasis(100)
      integer n
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      real q(50)
      real t1
      real w(50)
      real z(50)
      real z0

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z
c
c  set z0 equal to the most negative q(i).
c
      i = 1
      j = 2
1     if ( q(i) .le. q(j) ) go to 2
      i = j
2     j = j + 1
      if ( j .le. n ) go to 1
c
c  update q vector
c
      ir = i
      t1 = -q(ir)
      if ( t1 .le. 0.0e+00 ) go to 9

      do i = 1, n
        q(i) = q(i) + t1
      end do

      q(ir) = t1
c
c  update basis inverse and indicator vector
c  of basic variables.
c
      do j = 1, n
        b(j,ir) = -1.0e+00
        w(j) = q(j)
        z(j) = 0.0e+00
        mbasis(j) = 1
        l = n + j
        mbasis(l) = j
      end do

      nl1 = 1
      l = n + ir
      nl2 = ir
      mbasis(ir) = 3
      mbasis(l) = 0
      w(ir) = 0.0e+00
      z0 = q(ir)
      l1 = 1
c
c  print the initial almost complementary solution.
c
      write ( 6, 5 )
5     format ( 3(/), 5x, 29hinitial almost complementary ,
     &  8hsolution )
      do i = 1, n
        write ( 6, 6 ) i, w(i)
6       format ( 10x, 2hw(, i4, 2h)=, f15.5 )
      end do
      write ( 6, 8 ) z0
8     format ( 10x, 3hz0=, f15.5 )
      return
9     write ( 6, 10 )
10    format ( 5x, 36hproblem has a trivial complementary ,
     &  23hsolution with w=q, z=0. )
      ir = 1000
      return
      end
      subroutine newbas ( n )

c*********************************************************************72
c
c  purpose - to find the new basis column to enter in
c  terms of the current basis.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer i
      integer ir
      integer j
      integer l1
      integer mbasis(100)
      integer n
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      real q(50)
      real ti
      real w(50)
      real z(50)

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z
c
c  if nl1 is neither 1 nor 2 then the variable zo leaves the
c  basis indicating termination with a complementary solution.
c
      if ( nl1 .eq. 1 ) go to 2
      if ( nl1 .eq. 2 ) go to 5
      write ( 6, 1 )
1     format ( 5x, 22hcomplementary solution )
      call pprint ( n )
      ir = 1000
      return
2     ne1 = 2
      ne2 = nl2
c
c  update new basis column by multiplying by basis inverse.
c
      do i = 1, n
        ti = 0.0e+00
        do j = 1, n
          ti = ti - b(i,j) * am(j,ne2)
        end do
        a(i) = ti
      end do

      return

5     ne1 = 1
      ne2 = nl2
      do i = 1, n
        a(i) = b(i,ne2)
      end do

      return
      end
      subroutine sort ( n )

c*********************************************************************72
c
c  purpose - to find the pivot row for the next iteration by the
c  use of (simplex-type) minimum ratio rule.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer i
      integer ir
      integer l1
      integer mbasis(100)
      integer n
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      real t1
      real t2
      real q(50)
      real w(50)
      real z(50)

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z

      i = 1
1     if ( a(i) .gt. 0.0e+00 ) go to 2
      i = i + 1
      if ( i .gt. n ) go to 6
      go to 1

2     t1 = q(i) / a(i)
      ir = i

3     i = i + 1
      if ( i .gt. n ) go to 5
      if ( a(i) .gt. 0.0e+00 ) go to 4
      go to 3

4     t2 = q(i) / a(i)
      if ( t2 .ge. t1 ) go to 3
      ir = i
      t1 = t2
      go to 3

5     return
c
c  failure of the ratio rule indicates termination with
c  no complementary solution.
c
6     write ( 6, 7 )
7     format ( 5x, 37hproblem has no complementary solution )
      write ( 6, 8 ) l1
8     format ( 10x, 13hiteration no., i4 )
      ir = 1000
      return
      end
      subroutine pivot ( n )

c*********************************************************************72
c
cc pivot() performs the pivot operation by updating the
c  inverse of the basis and q vector.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer i
      integer ir
      integer j
      integer l
      integer l1
      integer mbasis(100)
      integer n
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      real q(50)
      real w(50)
      real z(50)

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z

      do i = 1, n
        b(ir,i) = b(ir,i) / a(ir)
      end do

      q(ir) = q(ir) / a(ir)

      do i = 1, n
        if ( i .ne. ir ) then
          q(i) = q(i) - q(ir) * a(i)
          do j = 1, n
            b(i,j) = b(i,j) - b(ir,j) * a(i)
          end do
        end if
      end do
c
c  update the indicator vector of basic variables.
c
      nl1 = mbasis(ir)
      l = n + ir
      nl2 = mbasis(l)
      mbasis(ir) = ne1
      mbasis(l) = ne2
      l1 = l1 + 1

      return
      end
      subroutine pprint ( n )

c*********************************************************************72
c
c  purpose - to print the current solution to complementary
c  problem and the iteration number.
c
      implicit none

      real a(50)
      real am(50,50)
      real b(50,50)
      integer i
      integer ir
      integer j
      integer k1
      integer k2
      integer l1
      integer mbasis(100)
      integer n
      integer ne1
      integer ne2
      integer nl1
      integer nl2
      real q(50)
      real w(50)
      real z(50)

      common am,q,l1,b,nl1,nl2,a,ne1,ne2,ir,mbasis,w,z

      write ( 6, 1 ) l1
1     format ( 10x, 13hiteration no., i4 )
      i = n + 1
      j = 1

2     k1 = mbasis(i)
      k2 = mbasis(j)
      if ( q(j) .ge. 0.0e+00 ) go to 3
      q(j) = 0.0e+00

3     if ( k2 .eq. 1 ) go to 5
      write ( 6, 4 ) k1, q(j)
4     format ( 10x, 2hz(, i4, 2h)=, f15.5 )
      go to 7

5     write ( 6, 6 ) k1, q(j)
6     format ( 10x, 2hw(, i4, 2h)=, f15.5 )

7     i = i + 1
      j = j + 1
      if ( j .le. n ) go to 2

      return
      end