function [xstar,ubflag,fconv,xconv] = affscale(A,b,c,x0,tol,beta)
%AFFSCALE Solves standard form linear programs using affine scaling alogorithm
%
%  Solves standard form linear programming problems:
%        min c'*x    subject to:   A*x = b 
%         x
%
%  Inputs:
%        A: Constraint coefficient matrix  (m x n)
%        b: Constraint right hand side vector (m x 1)
%        c: Cost vector (n x 1)
%       x0: Initial feasible solution (n x 1)
%      tol: Optimality tolerance
%     beta: Step-length scale factor
%
%  Outputs:
%    xstar: Optimum solution (if not found then [])
%   ubflag: Unboundedness flag (1 if unbounded; 0 otherwise)
%    fconv: Record of cost at each iteration (last value is optimal)
%    xconv: Record of x at each iteration (last value is optimal)
%
%  Nirav B. Shah 10/2002 for 6.251J

% check the inputs for errors

if (size(A,2) ~= size(x0,1))
    error('A matrix size does not agree with x0');
elseif (size(A,1) ~= size(b,1))
    error('A matrix size does not agree with b');
elseif (size(c) ~= size(x0))
    error('c vector size does not match x0');
elseif (size(x0,2) ~= 1)
    error('x0 must be a column vector');
elseif (size(c,2) ~= 1)
    error('c must be a column vector');
elseif (size(b,2) ~= 1)
    error('b must be a column vector');
elseif ~((sum(A*x0 == b)==length(b)) & (sum((x0 >= 0))==length(x0)))
    error('initial solution x0 is not feasible');
elseif (sum((x0 > 0)) < size(x0,1))
    error('initial solution x0 must strictly greater than 0');
end;

% all good lets start the algorithm
xk = x0;
xconv = [xk];
k=0;
fk = c'*xk;
fconv = [fk];
ubflag = 0;

e = ones(size(xk));
done = 0;

while (~done)
    % Computation of dual estimates and reduced costs
    Xk = diag(xk);
    pk = inv(A*Xk*Xk*A')*A*Xk*Xk*c;
    rk = c-A'*pk;

    % Optimality check

    if ((rk > -0.00001) & (e'*Xk*rk < tol))
        % Optimal solution found
        xstar = xk;
        done = 1;
%        disp('Optimal found at iteration:');
%        disp(k);
        continue;
    end;
    
    % Unboundedness check
    if (-Xk*Xk*rk >= 0)
        % Optimal cost is -infinity
        xstar = [];
        ubflag = 1;
        done = 1;
        disp('Unboundedness determined at iteration:');
        disp(k);
        continue;
    end;
    
    % Update the primal
    xk = xk - beta*(Xk*Xk*rk/norm(Xk*rk));
    xconv = [xconv xk];
    fconv = [fconv c'*xk];
    k=k+1;
end;