% This code computes a matrix in the convex hull of {A_1,...,A_m}
% with minimum condition number, where A_i are positive definite
% n by n matrices. It applies the alternating direction method to solve
%
% min  lambda_max(X)
% s.t. X - I >= 0, y >= 0, X = y_1*A_1 + ... + y_m*A*m.
%% 
% Input
% M       - a matrix whose ith column consists of the vector formed by 
%           stacking the columns of A_i one by one
% tol     - tolerance parameter
% maxiter - maximum number of iterations allowed
%%
% Output
% t       - approximate minimum condition number
% X       - matrix in the convex hull with cond(X) = t
% time    - cpu time
% iter    - number of iterations
%%
% Subroutines called:
% cond_update_X, cond_update_y, simplex_proj
%%
% March 5, 2011
% Written by Zhaosong Lu and Ting Kei Pong

function [t,X,time,iter] = condno_adm(M,tol,maxiter)

tic
% parameters
n = sqrt(size(M,1));
m = size(M,2);
beta = 1;
fprintf(' beta = %g\n', beta)
M = @(x,mode) handleMatrix(M,x,mode);

% initialization
y = zeros(m,1);
My = reshape(M(y,1),n,n);
w = ones(n,1)/n/beta;
Q = -eye(n)/n;

iter = 1;

while iter <= maxiter
  
  [X,w] = cond_update_X(My,n,w,Q,beta,tol/100);
  X = (X + X')/2;
  
  y = cond_update_y(M,My,n,X,y,Q,beta,tol/100);
  My = reshape(M(y,1),n,n);
  My = (My + My')/2;
  
  Q = Q + beta*(X - My);
  Q = (Q + Q')/2;
  
  % compute dual function value
  tmpeig = eig(My);
  pval = max(tmpeig);
  pfeas = (1 - min(tmpeig))/max(1,norm(tmpeig));
    
  q = eig(Q);
  dval = 1 + sum(q);
  dfeas = max([M(reshape(Q,n^2,1),2);sum(max(-q,0)) - 1])/max(1,norm(q));
  
  if abs(pval - dval)/max(1,abs(pval)) < 10*tol && max(dfeas,pfeas) < tol
    fprintf(' Terminate:\n iter = %4d, pval = %4.3f, dval = %4.3f, pfeas = %4.3f, dfeas = %4.3f, time = %4.1f\n', iter, pval, dval, pfeas, dfeas, toc)
    break
  end
  
  iter = iter + 1;
%   if mod(iter,100) == 0
%     fprintf(' Xiter = %4d, yiter = %4d\n', Xiter, yiter)
%     fprintf(' iter = %4d, pval = %4.3f, dval = %4.3f, pfeas = %4.3f, dfeas = %4.3f, time = %4.1f\n', iter, pval, dval, pfeas, dfeas, toc)
%   end
  
end

X = My/sum(y);
t = cond(X);
fprintf(' cond number = %4.1f\n', t)
time = toc;
  