// Ecole Centrale Paris
// Numerical Optimization and Applications
// 
// Resolution of linear program
// By the  SIMPLEX METHOD
// Written by Tahar Z. Boulmezaoud



function [info, sol, val ] = simplex(c, A, b, iter)
  
  
//   This function solves the linear program
//       Max c'*x
//       x >= 0
//       A x <= b
//   where we assume that b > = 0 


// Input parameters
//    A  : a m*n rectangular matrix   
//    b  : a column or row vector of size "m" 
//    c  : a column or row vector of size "n"
//  iter : the maximum number of iteration 
 
// Output parameters
//    info : an integer which takes the following values
//            0 : if the problem has at least one solution  
//            1 : if the objective function is unbounded  
//            2 : the method do not converge after 
//               "iter" iteration  
//            3 : the sizes of A, b and c are not coherent
//                (error)          
//    sol  :  a column vector containing the solution when info = 0.
//            In the other cases (info > 0), sol is empty
//    val  : the optimal value of the objective function. 
//           (val = 0 if info > 0)



//    Sizes  of A, b and c
      sa = size(A);
      nc = sa(1);   
      nv = sa(2);

      sb = size(b); 
      sc = size(c); 
      
//    Rearrangment of b and c into column vectors
      rhs = b; 
      cf = c;
      if (sb(1) == 1) then, rhs = b', end; 
      if (sc(1) == 1) then,  cf  = c', end; 
      sb = size(rhs); 
      sc = size(cf);       
      
           
//   Checking the coherence of sizes 

      if ((sb(1) <> nc) | (sc(1) <> nv)) then
          info = 3;
          sol = cell();
          val = 0.0; 
          print (%io(2), 'Error : The sizes of A, b and c are not coherent. '); 
          [info, sol, val]= return; 
      end




// Treating the obvious case c <= 0


     if (max(cf) <= 0) then
       sol = zeros(1:nv); 
       info = 0;
       vaml = 0;
       print (%io(2), 'This is the trivial case: the solution is 0'); 
       [info, sol, val]= return; 
     end
//       ---------
//  The simplex algorithm
 
 
     M = [A eye(nc, nc) rhs; cf' zeros(1:nc+1)];
     print (%io(2), '      ----- Starting tableau ------       ');        
     print (%io(2), M); 
     print (%io(2),'            ----------             ');    
     
     
     niter = 0; 
     sol = cell();
     val = 0;
      
//   The initial corner
     for i=1:nc
       crn(i) = nv + i;
     end   
 
     niter = 0; 
     stp = 0; // The stopping parameter
     

 
     
     
     while (( stp <> 1) & (niter < iter)) 
       niter = niter + 1; 
       cf = M(nc+1,:); 
      
       if (max(cf) < = 0)  then
//        Stop: a solution is found
          info = 0; 
          print (%io(2), 'A solution is found'); 
          sol = zeros(1:nv)';
          for i=1:nc 
            if (crn(i) <= nv) then
              sol(crn(i)) = M(i, nc+nv+1)/M(i, crn(i)); 
            end; 
          end        
  
          val = - M(nc+1, nc+nv+1);
          stp = 1; 
       else  
         ord = order_posit(cf);
      
         k = 1; 
         choice = 0; 
         so = size(ord);
         while ((k <= so(1)) & choice == 0)

//         The choice of the column 
           col = ord(k);  
           if ( max(M(1:nc, col)) > 0  )  then
//           the column contains some positive values
             jpiv = col; 
             choice = 1; 
             theta = 0; 
             for i=1:nc
                if ((theta == 0) & (M(i, col) > 0) ) then
                    ipiv = i;
                    theta = M(i, nv+nc+1)/M(i, col);  
                elseif (M(i, nv+nc+1) < theta*M(i, col)) then
                    theta = M(i, nv+nc+1)/M(i, col);  
                    ipiv = i;
                end
             end           
           else
//         The column does not contain any positive value
             k = k + 1; 
           end 
       end // end of while 
       
       if (k >  so(1)) then
//         In this case, the objective function
//          is unbounded
           
           stp = 1;
           info = 1;
           print (%io(2), 'The objective function is unbounded. '); 
           return;
       else
//       making the elimination 
         M = eliminations(M, ipiv, jpiv); 
         crn(ipiv) = jpiv; // change the exterme point
       end
       
     end // end of if 

    end // end of while 
 
//   
     if (niter >= iter) then
//      The number of iterations is non sufficient or
//      there is a cycle
        info = 2; 
        print (%io(2), 'No convergence  '); 
     end
      
     return; 
      
endfunction   

    
//------------------------------------------------------------------------
       
function [ord] = order_posit(vec)

//   This function return an integer vector ord containing
//   the indices of strictly positive components of vec
//   orderder in the decreasing order
 
        sz = size(vec);
        nb = sz(1);
        if (nb == 1) then, nb = sz(2), end; 

       val = vec(1); 
       nbp = 0;
       for i = 1:nb
         if (vec(i) > 0) then
            nbp = nbp + 1; 
            ind(nbp) = i;
         end
       end 

       for k = 1:nbp
        for m=k+1:nbp
          if (vec(ind(m)) > vec(ind(k))) then 
              i = ind(k); 
              ind(k) = ind(m);
              ind(m) = i; 
          end
        end
       end
       ord = ind;

endfunction

//------------------------------------------------------------------------

     function [Mat] = eliminations(M, ip, jp)
//     Making eliminations in the matrix M, suing the pivot element (ip, jp)        

         Mat =  M; 
         sz = size(M);

         for j=1:sz(2)
              Mat(ip, j) = M(ip, j)/M(ip, jp); 
         end 
         
        for i = 1:sz(1)
         if (i <> ip) then 
           facto = M(i,jp)/M(ip, jp); 
           Mat(i,:) = M(i,:) - facto*M(ip, :);
         end
       end
       pivot = [ip jp]; 
       print (%io(2), pivot); 
       print (%io(2), Mat); 
       print (%io(2),'            ----------             '); 
    endfunction

   
       