
  function [d,xc,io,iter] = ep21(l,u,A,itbound)
%
% This is a Ellipsoid-Scaling Algorithm (Algorithm 2.1) for:   
% problem:    finding x s.t.  l <= Ax <= u
%Input variables
%       d: initial working point
%       l: left-hand-side
%       u: right-hand-side        
%       A: constraint matrix
%       itbound: bound on total # of iterations.itbound=0 implies itbound=35
%
%Output variables
%       d: the feasible weights
%       xc: the center of the Ellipsoid corresponding to dopt
%       io: index for the feasibility of xc, io=1 means xc is feasible
%       
%Written by Aiping Liao, Cornell University, Aug., 1993
%
[m,n] = size(A);
if itbound == 0
   itbound = 35;
end
iter = 0; d = ones(m,1);
%
% Major loop
while ( iter < itbound)
 iter = iter + 1;
 r = 0.5*(u+l); D = diag(d); ADA = A'*D*A;
 xc = (ADA)\(A'*(D*r));
 f = xc'*(ADA)*xc - l'*D*u;
%
% Check feasibility.
%
 io = 1;
 for i=1:m
     if (A(i,:)*xc > u(i)) | (A(i,:)*xc < l(i))
	io = 0;
	j = i;
	break;
     end
 end
 if io == 1
    'feasible solution found'
    break;
 end
% 
% Calculate betta and gamma.
%
 sj = 0.5*(u(j)-l(i));  sj2 = sj*sj;
 gamma = A(j,:)*(ADA\A(j,:)'); betta = sj2/(f*gamma);
%
% Update
%
 if betta <= 0.25
    lambda = 1.0/(2*n*betta*gamma);
    d(j) = d(j)+ lambda;
 end
 if betta > 0.25
    lambda = 2/(n*gamma);
    d(j) = d(j) + lambda;
    d0 = lambda; ajxc = A(j,:)*xc;
    if ajxc < l(j)
       l0 = l(j); u0 = l(j) + sqrt(gamma);
    end
    if ajxc > u(j)
       l0 = u(j) - sqrt(gamma); u0 = u(j);
    end
    if ajxc > u(j)
       l(j) = (d(j)*l(j) + d0*l0)/(d(j) + d0);
    end
    if ajxc < l(j)
       u(j) = (d(j)*u(j) + d0*u0)/(d(j) + d0);
    end
 end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% ep41.m %%%%%%%%%%%%%%%%%%%%%%%%
  function [d,xc,io,iter] = ep41(l,u,A,itbound)
%
% This is a Ellipsoid-Scaling Algorithm (Algorithm 4.1) for:   
% problem:    finding x s.t.  l <= Ax <= u
%Input variables
%       d: initial working point
%       l: left-hand-side
%       u: right-hand-side        
%       A: constraint matrix
%       itbound: bound on total # of iterations.itbound=0 implies itbound=35
%
%Output variables
%       d: the feasible weights
%       xc: the center of the Ellipsoid corresponding to dopt
%       io: index for the feasibility of xc, io=1 means xc is feasible
%
%Written by Aiping Liao, Cornell University, Aug., 1993
%
[m,n] = size(A);
if itbound == 0
   itbound = 35;
end
iter = 0; d = ones(m,1);
%
% Major loop
while ( iter < itbound)
 iter = iter + 1;
 r = 0.5*(u+l); D = diag(d); ADA = A'*D*A;
 xc = (ADA)\(A'*(D*r));
 f = xc'*(ADA)*xc - l'*D*u;
%
% Check feasibility.
%
 io = 1;
 for i=1:m
     if (A(i,:)*xc > u(i)) | (A(i,:)*xc < l(i))
	io = 0;
	j = i;
	break;
     end
 end
 if io == 1
    'feasible solution found'
    break;
 end
% 
% Calculate betta and gamma.
%
 ajxc = A(j,:)*xc;
 sj = 0.5*(u(j)-l(i));  sj2 = sj*sj;
 gamma = A(j,:)*(ADA\A(j,:)'); betta = sj2/(f*gamma);
%
% Update
%
 if betta <= 0.25
    theta0 = (ajxc - l(j))*(ajxc - u(j))/sj2;
    rho = 2*betta - 1/n + (betta*theta0)/n;
    theta = rho*rho + 4*betta*(1-1/n)*(betta*theta0+1/n);
    theta = (sqrt(theta) - rho)/(2*betta*(1-1/n));
    lambda = theta/gamma;
    d(j) = d(j)+ lambda;
 end
 if betta > 0.25
    lambda = 2/(n*gamma);
    d(j) = d(j) + lambda;
    d0 = lambda; 
    if ajxc < l(j)
       l0 = l(j); u0 = l(j) + sqrt(gamma);
    end
    if ajxc > u(j)
       l0 = u(j) - sqrt(gamma); u0 = u(j);
    end
    if ajxc > u(j)
       l(j) = (d(j)*l(j) + d0*l0)/(d(j) + d0);
    end
    if ajxc < l(j)
       u(j) = (d(j)*u(j) + d0*u0)/(d(j) + d0);
    end
 end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% ep42.m %%%%%%%%%%%%%%%%%%%%%
  function [d,xc,io,iter] = ep42(l,u,A,itbound)
%
% This is a Ellipsoid-Scaling Algorithm (Algorithm 4.2) for:   
% problem:    finding x s.t.  l <= Ax <= u
%Input variables
%       d: initial working point
%       l: left-hand-side
%       u: right-hand-side        
%       A: constraint matrix
%       itbound: bound on total # of iterations.itbound=0 implies itbound=35
%
%Output variables
%       d: the feasible weights
%       xc: the center of the Ellipsoid corresponding to dopt
%       io: index for the feasibility of xc, io=1 means xc is feasible
%
%Written by Aiping Liao, Cornell University, Aug., 1993
%
[m,n] = size(A);
if itbound == 0
   itbound = 35;
end
iter = 0; d = ones(m,1);
%
% Major loop
while ( iter < itbound)
 iter = iter + 1;
 r = 0.5*(u+l); D = diag(d); ADA = A'*D*A;
 xc = (ADA)\(A'*(D*r));
 f = xc'*(ADA)*xc - l'*D*u;
%
% Check feasibility.
%
 io = 1;
 for i=1:m
     if (A(i,:)*xc > u(i)) | (A(i,:)*xc < l(i))
	io = 0;
	j = i;
	break;
     end
 end
 if io == 1
    'feasible solution found'
    break;
 end
% 
% Calculate betta and gamma.
%
 ajxc = A(j,:)*xc;
 sj = 0.5*(u(j)-l(i));  sj2 = sj*sj;
 gamma = A(j,:)*(ADA\A(j,:)'); betta = sj2/(f*gamma);
%
% Update
%
 if betta <= 0.25
    theta0 = (ajxc - l(j))*(ajxc - u(j))/sj2;
    rho = 2*betta - 1/n + (betta*theta0)/n;
    theta = rho*rho + 4*betta*(1-1/n)*(betta*theta0+1/n);
    theta = (sqrt(theta) - rho)/(2*betta*(1-1/n));
    lambda = theta/gamma;
    d(j) = d(j)+ lambda;
 end
 if betta > 0.25
    dbar = [d(1:j-1); d(j+1:m)]; Dbar = diag(dbar);
    rbar = [r(1:j-1); r(j+1:m)];
    Abar = [A(1:j-1,:); A(j+1:m,:)];
    ADAbar = Abar'*Dbar*Abar;
    xcbar = (ADAbar)\(Abar'*(Dbar*rbar));
    gammabar = A(j,:)*(ADAbar\A(j,:)');
    zl = xcbar - sqrt(gammabar) *(ADAbar\A(j,:)');
    zu = xcbar + sqrt(gammabar) *(ADAbar\A(j,:)');
    if ajxc < l(j)
       lbar = l(j); u(j) = A(j,:)*zu;
    end
    if ajxc > u(j)
       l(j) = A(j,:)*zl; ubar = u(j);
    end
    theta0 = (ajxc - l(j))*(ajxc - u(j))/sj2;
    rho = 2*betta - 1/n + (betta*theta0)/n;
    theta = rho*rho + 4*betta*(1-1/n)*(betta*theta0+1/n);
    theta = (sqrt(theta) - rho)/(2*betta*(1-1/n));
    lambda = theta/gamma;
    d(j) = d(j)+ lambda;
 end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% ep43.m %%%%%%%%%%%%%%%%%%%%%%
  function [d,xc,io,iter] = ep43(l,u,A,itbound)
%
% This is a Ellipsoid-Scaling Algorithm (Algorithm 4.3) for:   
% problem:    finding x s.t.  l <= Ax <= u
%Input variables
%       d: initial working point
%       l: left-hand-side
%       u: right-hand-side        
%       A: constraint matrix
%       itbound: bound on total # of iterations.itbound=0 implies itbound=35
%
%Output variables
%       d: the feasible weights
%       xc: the center of the Ellipsoid corresponding to dopt
%       io: index for the feasibility of xc, io=1 means xc is feasible
%
%Written by Aiping Liao, Cornell University, Aug., 1993
%
[m,n] = size(A);
if itbound == 0
   itbound = 35;
end
iter = 0; d = ones(m,1);
ll = l; uu = u;
%
% Major loop
while ( iter < itbound)
 iter = iter + 1;
 r = 0.5*(u+l); rr = 0.5*(ll+uu);  D = diag(d); ADA = A'*D*A;
 xc = (ADA)\(A'*(D*rr));
 f = xc'*(ADA)*xc - ll'*D*uu;
 d = d/f; D = diag(d); ADA = ADA/f;
%
% Check feasibility.
%
 io = 1;
 for i=1:m
     if (A(i,:)*xc > u(i)) | (A(i,:)*xc < l(i))
	io = 0;
	j = i;
	break;
     end
 end
 if io == 1
    'feasible solution found'
    break;
 end
% 
% Calculate gamma.
%
 gamma = A(j,:)*(ADA\A(j,:)'); 
%
% Update
%
    lambda = 1/(2*n*gamma);
    d(j) = d(j) + lambda;
    d0 = lambda; ajxc = A(j,:)*xc;
    if ajxc < l(j)
       l0 = l(j); u0 = A(j,:)*xc + sqrt(gamma);
    end
    if ajxc > u(j)
       l0 = A(j,:)*xc - sqrt(gamma); u0 = u(j);
    end
    alpha = d0/(d0+d(j)); delta = 1 - alpha;
    r0 = 0.5*(l0+u0); rj = 0.5*(ll(j)+uu(j)); t0 = l0*u0; tj = ll(j)*uu(j);
    temp1 = alpha*r0 + delta*rj; temp2 = alpha*t0 + delta*tj;
    temp = (temp1)^2 - temp2;
    temp = sqrt(temp);
    ll(j) = temp1 - temp;   uu(j) = temp1 + temp;
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% ep44.m %%%%%%%%%%%%%%%%%%%%%%%
  function [d,xc,io,iter] = ep44(l,u,A,itbound)
%
% This is a Ellipsoid-Scaling Algorithm (Algorithm 4.4) for:   
% problem:    finding x s.t.  l <= Ax <= u
%Input variables
%       d: initial working point
%       l: left-hand-side
%       u: right-hand-side        
%       A: constraint matrix
%       itbound: bound on total # of iterations.itbound=0 implies itbound=35
%
%Output variables
%       d: the feasible weights
%       xc: the center of the Ellipsoid corresponding to dopt
%       io: index for the feasibility of xc, io=1 means xc is feasible
%
%Written by Aiping Liao, Cornell University, Aug., 1993
%
[m,n] = size(A);
if itbound == 0
   itbound = 35;
end
iter = 0; d = ones(m,1);
ll = l; uu = u;
%
% Major loop
while ( iter < itbound)
 iter = iter + 1;
 r = 0.5*(u+l); rr = 0.5*(ll+uu);  D = diag(d); ADA = A'*D*A;
 xc = (ADA)\(A'*(D*rr));
 f = xc'*(ADA)*xc - ll'*D*uu;
 d = d/f; D = diag(d); ADA = ADA/f;
%
% Check feasibility.
%
 io = 1;
 for i=1:m
     if (A(i,:)*xc > u(i)) | (A(i,:)*xc < l(i))
	io = 0;
	j = i;
	break;
     end
 end
 if io == 1
    'feasible solution found'
    break;
 end
% 
% Calculate gamma.
%
 gamma = A(j,:)*(ADA\A(j,:)'); 
%
% Update
%
    lambda = 1/(2*n*gamma);
    d(j) = d(j) + lambda;
    d0 = lambda; ajxc = A(j,:)*xc;
    if ajxc < l(j)
       l0 = ajxc; u0 = A(j,:)*xc + sqrt(gamma);
    end
    if ajxc > u(j)
       l0 = A(j,:)*xc - sqrt(gamma); u0 = ajxc;
    end
    alpha = d0/(d0+d(j)); delta = 1 - alpha;
    r0 = 0.5*(l0+u0); rj = 0.5*(ll(j)+uu(j)); t0 = l0*u0; tj = ll(j)*uu(j);
    temp1 = alpha*r0 + delta*rj; temp2 = alpha*t0 + delta*tj;
    temp = (temp1)^2 - temp2;
    temp = sqrt(temp);
    ll(j) = temp1 - temp;   uu(j) = temp1 + temp;
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
