function [xopt, fopt,Iopt]=lp(A,b,c)
% This function takes as input a linear program in standard form and
% identifies the basic feasible solution with the highest value of the
% objective function (i.e., the best vertex of the feasible region). 
% It creates the matrix of all combinations of {1,...,n} of size m. For
% each combination I it finds the corresponding basic solution (if the
% corresponding matrix B is invertible). It then checks if the basic
% solution is feasible. If it is, it compares the value of the objective
% function with that of the current maximum. If it is higher it is replaced
% by the current value. 
%Inputs : A : mxn matrix, with m<=n, rank(A)=m
%b : column vector of size m, b>=0
%c: column vector of size n
%
%Results: xopt=optimal solution, fopt=optimal value, Iopt=vector of indices
%

     Iopt=0;
     fopt=-10^12;
     xopt=-1;
     [m,n]=size(A);
     M=combnk(1:n,m);
     K=size(M,1);
     for k=1:K
         I=M(k,:);
         [x0,f0]=basicsol(A,b,c,I);
         if sum(x0<0)==0
             if f0>fopt
                 xopt=x0;
                 fopt=f0;
                 Iopt=I;
             end;
         end;
     end;
 end
 