clear all

% read in the raw data for the trnsport model
trdata

% now do the modeling work in Matlab - set up canonical A,b,c
M = pN + dN;
N = pN * dN;
NNZ = (pN * dN) * 2;
c = zeros(N,1);
b = zeros(M,1);
A = zeros(NNZ,3);
k = 0;
for p = 1:pN
    for d = 1:dN
        j = (p-1)*dN + d;
        c(j,1) = f * distance(p,d) / 1000;
        k = k + 1;
	A(k,1) = p;   % row index
	A(k,2) = j;   % col index 
	A(k,3) = 1;   % coef value
        k = k + 1;
	A(k,1) = pN + d;  % row index
	A(k,2) = j;       % col index 
	A(k,3) = -1;      % coef value
   end;
end;
i = 0;
for p = 1:pN
    i = i + 1;
    b(i,1) = cap(p);
end;
for d = 1:dN
    i = i + 1;
    b(i,1) = -dem(d);
end;
Is.name = 'I';
Is.val  = [ linspace(1,M,M) ]';
Is.type = 'set';
Is.load = 'initialize';
Js.name = 'J';
Js.val  = [ linspace(1,N,N) ]';
Js.type = 'set';
Js.load = 'initialize';
As.name = 'A';
As.val  = A;
As.type = 'parameter';
As.form = 'sparse';
bs.name = 'b';
bs.val  = b;
bs.form = 'full';
bs.dim  = 1;
cs.name = 'c';
cs.val  = c;
cs.form = 'full';
cs.dim  = 1;
[s,z,x] = gams('tr3_mod',Is,Js,As,bs,cs);

uels = s.uels{1};
s1   = s.val(:,1);

r = find(strcmp(uels,'modelstat'));
if norm(size(r)-[ 1 1 ]) > 0
  error('could not find modelstat');
end
i = find(s1==r);
modelstat = s.val(i,2);

r = find(strcmp(uels,'solvestat'));
if norm(size(r)-[ 1 1 ]) > 0
  error('could not find solvestat');
end
i = find(s1==r);
solvestat = s.val(i,2);

if norm(modelstat - 1) > 0
  error('bad modelstat, expected 1');
end
if norm(solvestat - 1) > 0
  error('bad solvestat, expected 1');
end

if(z.dim ~= 0)
  error('expected z to have dimension 0');
end
if norm(z.val - z_) > 1e-6
  error('bad z.val');
end
z.val;

fprintf('modelstat = %d (1 = optimal solution)\n', modelstat);
fprintf('solvestat = %d (1 = normal completion)\n\n', solvestat);
fprintf('optimal value = %g\n\n', z.val);
x.val
