  function  [X,y,Z,muu,itersinfo,meritfn]=mcsdp(toler,crosstoler)
%           [X,y,Z,muu,itersinfo,meritfn]=mcsdp(toler,crosstoler)
%
%
%    MCSDP solves the SDP relaxation of the max-cut problem
%%      max trace QX    s.t.  diag(X)=ones     X psd
%%  The algorithm solves the equivalent perturbed optimality conditions:
%%   ( usMat(x,n) +I)(Diag(y)-Q) - muu I = 0, 
%%                           X=( usMat(x,n) +I); Z=(Diag(y)-Q); muu >0 (all pos def)
%%         and simultaneously (adaptively) drives muu to 0.
%%  URL:
%%
%   This file mcsdp.m can be tested using the file testmcsdp.m
%   This file calls:
%  mcsdp.m ->  solvlss.m -> MSlsqr.m -> MSFprecondfun.m ->
%                            zops.m,xops.m,zop.m,xop.m, usMat.m, usvec.m
%%
%% Input:
%%
%%    The following line, with global variables, is needed in the call pgm!
%
%             global Xu Q q2kron indsu ;
%
%%    Q  -  given real symmetric n by n matrix (global variable!!)
%%            (e.g. obtained from the weights from the max-cut problem)
%%   toler - relative tolerance for stopping criteria of mcsdp pd-ip algorithm
%
%    crosstoler -      (optional input) tolerance for start of 
%                         crossover to sigmaa=0 and no pos. check. 
%                       Default is .1 (see below)
%
%% Output:
%%     X   -  optimum correlation (for max-cut) matrix
%%     y,Z   -  optimum dual vector,matrix
%%     muu  -  duality gap at end  (trace ZX)
%%   itersinfo -  2 dim vector with
%%                 (number of iterations used, crossover iteration no.)
%%   meritfn - if included as desired input (optional), this is a 
%%             vector with the value of ||F||| at each iteration
%%             and a plot at completion.
%%
% Parameters:
%      iterbnd   - upper bound on number of iterations
%
% DEVELOPMENTS:
% August 6, 2001:
%          Start to modify and develope the algorithm.
%              
% August 7, 2001:
%          Add correct operators and adjoints and diagonal preconditioner
%          Just preliminary operator operations, i.e. symmetry 
%            structure should be exploited in multiplications!
%          Next time: Write mfiles for operators and adjoints F',  F'*
% August 9, 2001:
%          Adaptive change in tolerance for lsqr - using sqrt(muu) value,
%          i.e. inexact GN algorithm.
% August 10, 2001:
%          working on preconditioner
% August 11, 2001:
%            Testing on the accuracy of the lsqr solve, i.e.
%            this seems to depend on the residual as well as nnz
% August 11, 2001:
%             Diagonal Preconditioning is working well. 
%             But precond fails when the adaptive atol,btol is used???
% August 16, 2001:
%            Start working on crossover.
% August 21, 2001:
%            improve preconditioning calculations in solvlss.m          
% Sept. 9, 2001:
%            start testing on SDPLIB probs - remove zero cols in Q
% Sept. 13, 2001:
%            improvements on the *op*.m functions and on usMat use.
% Sept. 15, 2001:
%            improvements to reduce time spent in MSFprecondfun, i.e
%             xu and Qy are updated in solvlss
%          

global Xu Q q2kron indsu ;

norig=size(Q,1);n=norig;
Qorig=Q;
if ~issparse(Q),
  disp('to User from mcsdp.m: Inputted matrix Q of weights was not sparse - BE AWARE!! ');
  disp('Q is made sparse for the algorithm.');
  disp('If Q is not sparse, then another version  of the algor. should be used');
  Q=sparse(Q);
end


y=(sum(abs(Q))');    % initial estimate of dual column vector
indsQ0=find(y==0);      % avoid 0 cols of Q 
Qcolszero=0;
if sum(indsQ0)>0,
  Qcolszero=1;
  indsnnz=1:n;indsnnz=setdiff(indsnnz,indsQ0);   % nonzero cols
  y=y(indsnnz);
  Q=Q(:,indsnnz);Q=Q(indsnnz,:);
  disp(['Warning from mcsdp.m: zero columns in Q; User can simplify original  MC problem by using disconnected graphs']);
  disp([num2str(length(indsQ0)),'  zero columns in Q removed during algorithm;']);
  n=size(Q,1);
end

dQ=diag(Q);doptval=sum(dQ);
diagnz=0;
if norm(dQ)>0,
  diagnz=1;
  disp('to User from mcsdp.m beware: ');
  disp('diagonal of input matrix Q was set to 0 for the algorithm, since it is not needed');
     Q=Q-diag(dQ);
end
temp=Q-Q';
if temp(:)'*temp(:)>eps,
     disp('to User from mcsdp.m beware: ');
     disp('Q-Q^T~=0 in function mcsdp');
     disp('Q is symmetrized for the user');
     Q=(Q-Q')/2;
end
clear temp
%  Parameters:  % conservative at the start
sigmaa=1;      % centering parameter - changes adaptively
step=1;        % initial stepsize - changes adaptively
              %  initially was 2 - using 1 for cond. numb tests
tau1=.05;      % initial backtrack value
           % using: step=min(.99997,1-min(tau1,abs(tau1*gapp)   )    )*step;
iterbnd=100;   % bound on number of major iterations

%% Initializations
X=zeros(n);
Xu=zeros(n);          % Xu=usMat(x); updated in solvlss before MSlsqr call
nnzQ=nnz(Q);
Q2=spalloc(n,n,nnzQ);
n2=n^2;tn=n*(n+1)/2;tnm1=n*(n-1)/2;
indsu=find(triu(ones(n),1));    % strict upper triangular indices
indsoff=setdiff(1:n2,1:n+1:n2);  % off diagonal indices
q2=Q(indsu);
                 % for opt. diagonal preconditioning
  q2=q2.*q2;
  Q2(indsu)=q2;
  Q2=Q2+Q2';
  q2=sum(Q2)';       % vector of column norms
  inds=[1:n];
  inds=diag(inds)*ones(n);
  indsprec=zeros(tnm1,2);
  indsprec(:,1)=inds(indsu);
  inds=inds';
  indsprec(:,2)=inds(indsu);
  q2kron=q2(indsprec(:,1))+ q2(indsprec(:,2));
  clear inds Q2 q2

x=zeros(tnm1,1);
X=usMat(x,n)+eye(n);
y=2*y;
%Z=sparse(diag(y)-Q);     % to ensure dual feas 
Z=spdiags(y,0,-Q);       %stays sparse??? - efficient???
%%      better feasible solution???? how??? Rendl approach???
%%% X,y,Z are feasible primal-dual variables; X is represented by x


%% check dual feasibility   - temporary
if  min(eig(Z)) < 0,
   'check initial feas'
   keyboard
end
ZX=Z*X;
wantmeritfn=nargout>=6;
if nargin<3,
   crosstoler=.1;    % default tolerance to start crossover
end
gapp=trace(ZX);
muu=gapp/n;
noiter=0;
optval=trace(Q*X);optval=optval+doptval;   % adjust for diag(Q)
crossover=0;          % indicator for crossover algor.
itersinfo(2)=0;       % in case no crossover takes place
ZXnorm=norm(ZX,'fro');    %% for crossover
if wantmeritfn,meritfn=[];end

      disp('noiter   -log10(relgapp)    step        sigmaa      optval         cputime.dir');
while ( ~crossover*(gapp/(abs(optval)+1) > toler)  & (noiter<iterbnd)  )...
         |(crossover*(ZXnorm)/(abs(optval)+1) > toler);

   noiter=noiter+1;
   % compute RHS at current iterate and current muu estimate
   Fmu=  ZX-sigmaa*muu*eye(n)   ;
   % find the search direction at the current iterate
   tfordirc=cputime;
   % find a least squares soln F'deltaxy=-Fmu
   [deltax,deltay]=solvlss(x,y,-Fmu(:),muu,...
                    noiter,nnzQ,toler,indsprec);
   dirtimec=cputime-tfordirc;
%%%%%%%%%%%%%%%%%%%%%%%
%%% start crossover - needs lots of polishing
  if   (gapp/(abs(optval)+1) < crosstoler)  & ~crossover ,
     crossover=1;
     step=1;      % step stays 1 after crossover starts
     sigmaa=0;     % stays 0 after crossover starts
     itersinfo(2)=noiter;
      disp(['*********CROSSover started at: ',num2str(gapp/(abs(optval)+1))]);
      disp('noiter -log10(relnormZX)     step       sigmaa     optval          cputime.dir,');
  end
  if ~crossover,
%     disp(['starting the back tracking  ']);
     Xp=X+step* usMat(deltax,n);
     dy=step* deltay;
     Zp=Z+spdiags([ dy],0,n,n);
  %%%%%%%%%change to Lanczos? then do not form matrices??? space saving??
     [Rx,choltrueX]=chol(Xp);           % redo eff using Lanczos! ???
     [Rz,choltrueZ]=chol(Zp);           % redo eff using Lanczos!
     while choltrueX+choltrueZ>0,
        step=.8*step;
        Xp=X+step* usMat(deltax,n);
        dy=step* deltay;
        Zp=Z+spdiags([ dy],0,n,n);
        [Rx,choltrueX]=chol(Xp);
        [Rz,choltrueZ]=chol(Zp);
     end
  %  step=.999997*step;     % original backtrack strategy
     step=min(.99997,1-min(tau1,abs(tau1*gapp)   )    )*step;
  %  step=.97*step;
   end
%%% end crossover modif - below for sigmaa as well.
   Xp=X+step* usMat(deltax,n);
   dy=step* deltay;
   Zp=Z+spdiags([ dy],0,n,n);
%  dirtimeb=cputime-dirtimeb;
   X=Xp;
   %Xu(indsoff)=X(indsoff)/sqrt(2);
   %keyboard
   x=usvec(X,n);
   Z=Zp;
   y=y+dy;
   ZX=Z*X;
   ZXnorm=norm(ZX,'fro');    %% for crossover check
   gapp=trace(ZX);
   muu=gapp/(2*n);
   optval=trace(Q*X);optval=optval+doptval;   % adjust for diag(Q)
%  reler is relative error: duality-gapp/n1/(optval+1)
  if crossover,
  muu=ZXnorm/(abs(optval)+1);  % need positive for lsqr tolerance atol,btol
disp([sprintf(' %3g      %6.2e         %4g    %10g     %5.4e      %2.4e  %2.4e %2.4e', ...
       noiter,-log10(ZXnorm/(abs(optval)+1)), ...
        1,0,optval,dirtimec)]);
  else
disp([sprintf(' %3g      %6.2e         %3.1e    %2.1e     %2.4e      %2.4e  %2.4e %2.4e', ...
       noiter,-log10(gapp/(optval+1)), ...
        step,sigmaa,optval,dirtimec)]);
% update muu and sigmaa depending on stepsize
      if step < .001,muu = muu * 2.1;sigmaa=1;step=.9;  %muu increased
      elseif step < .01, muu = muu * 2;sigmaa=1;step=1;  % muu unchanged
      elseif step < .25, muu = muu * 1.5; sigmaa=1;step=1; 
%%   testing for influence of condition number on lsqr
%%% steplength restricted to one when using lsqr - seems to help
%%% the next 3 lines are changed to get steplength <=1
%% This appears to be TRUE in case of crossover soon after,
%%   i.e. stay well centered till crossover.
%     elseif step < .95, muu = muu * 1.2; sigmaa=1-.3*min(step,1);step=1.5;
%     elseif step < 1.05, muu = muu *.5; sigmaa=1-.3*min(step,1); step=2.0;
%     else  sigmaa=1-.3*step;step=2;muu=muu*.8;
      elseif step < .95, muu = muu * 1.2; sigmaa=1-.3*min(step,1);step=1.0;
      elseif step < 1.05, muu = muu *.5; sigmaa=1-.3*min(step,1); step=1.0;
      else  sigmaa=1-.3*step;step=1;muu=muu*.8;
      end
  end
%%% end crossover modif - see above as well

   if wantmeritfn,meritfn=[meritfn ZXnorm];end

end          % end of main while loop
   if wantmeritfn,
     plot(-log10(meritfn/(abs(optval)+1)));
     title('plot of -log10(norm(ZX)/(abs(optval)+1))')
   end
if noiter <=iterbnd,
  disp('to USER from mcsdp.m: Tolerance achieved - algorithm ends')
  disp(['for crossover; final rel. norm(ZX) is: ',num2str(ZXnorm/(abs(optval)+1))]); 
  disp(['for crossover; rel. min of eigs is: ',...
     num2str( min(  min(eig(X)),min(eig(Z)) )/(abs(optval)+1)  ) ]); 
else
  disp(['to USER from mcsdp.m: ERROR;  exceeded iteration bound:  ',...
                num2str(noiter)]);
end
itersinfo(1)=noiter;
  if diagnz,
    disp(['nonzero diagonal of Q is restored']);
    Q=Q+diag(dQ);
  end
  if Qcolszero,
  Q1=Q;y1=y;X1=X;Z1=Z;
  Q=spalloc(norig,norig,nnz(Q));y=zeros(norig,1);X=eye(norig);
  Z=spalloc(norig,norig,nnz(Z));
  Q(indsnnz,indsnnz)=Q1;y(indsnnz)=y1;X(indsnnz,indsnnz)=X1;
  Z(indsnnz,indsnnz)=Z1;
  disp(['zero cols of Q restored; X,y,Z adjusted']);
  n=norig;
  end
