function [cost,schedule] = trans_solver(C,s,d)
## Function [cost, schedule] = trans_solver(C, s, d) solves the 
## balanced transportation problem
##
##         -----   -----  
##         \       \
##    min   >       >     C(i,j)*X(i,j)
##         /       /
##         -----   -----
##           i       j
##
##           ----- 
##           \  
##    s.t.    >     X(i,j) = s(i)  for all ports i
##           /
##           -----
##             j
##
##           -----
##           \
##            >     X(i,j) = d(j)  for all markets j       
##           /
##           -----
##             i
##
## Input:
##
## C, s, d
##      The matrix of the transportation costs, the column vector 
##      of the supplies, and the column vector of the demands. 
##
## Output:
##
## cost, schedule
##      The minimal transportation cost and schedule.
##
## See also:  glpk.

## Some short-hands.
m1 = length(s);                                  # Number of ports.
m2 = length(d);                                  # Number of markets.
n = m1*m2;                                       # Number of decisions.
o = zeros(n,1);                                  # For glpk lower bounds.

## Build the technology matrix A
A = zeros(m1+m2,n);                              # Initialization.
for i = 1:m1                           
  A(i, (1:m2)+(i-1)*m2) = ones(1,m2);            # Supply side.
  A((1:m2)+m1, (1:m2)+(i-1)*m2) = eye(m2);       # Demand side.
endfor

## Set the sense of bounds (equality for balanced).
ctype = "";                           
for i=1:(m1+m2)                                    
  ctype = [ctype, "S"];               
endfor 

## Set the type of each decision variable as continuous.
vtype = "";                                  
for i=1:n                                   
  vtype = [vtype, "C"];                     
endfor 

## Solve the system by calling glpk.
[schedule, cost] = glpk(C'(:), A, [s; d], o, [], ctype, vtype, 1);

## Reshape the schedule vector to a matrix.
schedule = reshape(schedule, m2, m1)';
endfunction