function [z_max, x_max, status] = stu_lp_solver(c, A, b)
## Function [z_max, x_max, status] = stu_lp_solver(c, A, b) solves the LP 
##
##     max  c'*x
##     s.t.  A*x <= b
##             x >= 0, 
##
## by using glpk.
##
## Input:
##
## c, A, b
##      The (column) vector defining the objective function, the 
##       technology matrix, and the (column) vector of the constraints.
##
## Output:
##
## z_max, x_max
##      The maximal value of the objective and an optimal decision.
##
## status
##      A string indicating the status of the solution: 
##	    "bounded"
##           A bounded solution was found. 
##      "unbounded"
##           The solution is known to be unbounded.
##      "infeasible"
##           Solutions are known not to exist.
##      "unknown"
##           The algorithm failed.
##
## See also:  glpk.

## Get m, the number of constraints (excluding the sign constraints), and
## n, is the number of decision variables, from the technology matrix A.
[m,n] = size(A);

## Some input-checking: 
if ( nargin!=3 )
  error("stu_lp_solver: The number of input arguments must be 3.\n");
elseif ( columns(c)>1 )
  error("stu_lp_solver: Objective c must be a column vector.\n"); 
elseif (columns(b)>1 )
  error("stu_lp_solver: Upper bounds b must be a column vector.\n");
elseif (rows(c)!=n || rows(b)!=m )
  error("stu_lp_solver: The dimensions of c, A, and b do not match.\n");
endif 

## Set the parameters for glpk:

## Set the decision-wise lower bounds all to zero, i.e. build a zero
## column vector with n zeros. 
o = zeros(n,1);

## Set the sense of each constraint as an upper bound.
ctype = "";                                 # Start with an empty string. 
for i=1:m                                   # Loop m times.
  ctype = [ctype, "U"];                     # Append "U" to the string.
endfor 

## Set the type of each variable as continuous.
vtype = "";                                 # Start with an empty string. 
for i=1:n                                   # Loop n times.
  vtype = [vtype, "C"];                     # Append "C" to the string.
endfor 

## Solve the system by calling glpk.
[x_max, z_max, STATUS] = glpk(c, A, b, o, [], ctype, vtype, -1);

## Set the STATUS code given by glpk to the appropriate string
status = "unknown";                         # Pessimism by default.
if ( STATUS==180 )                          # Everything went fine.
  status="bounded";
elseif ( STATUS==183 )                      # LP infeasible detected.
  status="infeasible";
elseif ( STATUS==184 )                      # LP unboundedness detected.
  status="unbounded";
endif
endfunction