   function [optval,L,P,iter]=dualcomp(A,H,toler)
%  function [optval,L,P,iter]=dualcomp(A,H,toler)
%
% **Dual step** algorithm for the approximate matrix completion problem with
%  weight matrix H. A dual-dual interior-point method is used.
%  This program uses the dual step to evaluate the primal step.
%  It is more efficient if H has few 'sufficiently large' elements
%  relative to 0 elements, i.e. H is sparse.
%    Otherwise use the primal step algorithm
%
% User INPUT:  A -  the symmetric matrix A 
%                 (the sparse matlab package is assumed)
%              H -  weight matrix
%                  (diag(H) not 0 is ensured!!)
%                 i.e. 0 corresponds to free elements in A;
%                  >= 1000 corresponds to definitely fixed elements in A.
%                  The program fails if a definitely 
%                       fixed principal minor is not >0.
%          toler -  tolerance for duality gap 
%
%  OUTPUT:  optval - optimal value
%           L      - optimal value of dual variable Lambda 
%                      (and it equals the gradient of obj fn.)
%           P      - optimal primal value
%           iter   - iteration count
%
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
% DEVELOPMENTS:
% 30 Jul 2001: 
%              Attempted to use adaptive sigma-mu but did not succeed
%    in improving speed or accuracy. Perhaps predictor-corrector would
%    be better approach.
%    Extra accuracy is expensive (more iterations).
%%%%%%%%%%%%%%%%%%%%
% 29 Jul 2001: These are the first changes made after six years.
%              Using profile in matlab to try and improve efficiency.
%              Several 'sparse' lines commented out, i.e. time wasted
%               in converting a sparse variable to sparse.
%              Several 'full' statements were removed - definite speedup.
%              Several statements improved since using matlab 6.0, e.g.
%                    H=abs(H) rather than H=abs(full(H));
%              Changed from mineig1 to eigs matlab function for 
%                      finding the minimum eigenvalue.
%
%

% Initialization:
disp(['   ']);
disp([' ----------------------------- ']);
disp(['Start of new approximate completion problem using dual-first-step']);
n=size(A,1);
n2=n^2;
if ~issparse(H),
    disp(['H was not sparse; so made sparse -  WARNING']);
    H=sparse(H);
end
if min(min(H))<0,
  disp(['   ']);
  disp(['H had a negative element - WARNING']);
  disp(['Using abs(H) instead']);
  H=abs(H);
end
while sum(diag(H) <=0) > 0,
  disp(['   ']);
  disp(['diagonal of H had nonpositive elements - WARNING']);
  disp(['solutions of linear systems may become ill-conditioned']);
  disp(['diagonal is perturbed to be nonzero']);
  H=H+diag( (1000*eps)*rand(n,1) );
end

epsilon=toler;  % duality gap stopping tolerance
if ~issparse(A),
   A=sparse(A);
end
indzs=find(H==0);
indnzs=find(H);
[izs,jzs]=find(H);
disp(['The dimension n is: ',num2str(n)]);
disp(['The number of nonzero elements in H is: ',num2str(length(izs))]);
nnzsH= nnz(H); % number of nonzeros
H2=H.*H;    % Hadamard product
tH2=sparse(2*H2);   
%%%%%%%%%%%%%%%%%%%%%%
%%%%% heuristic for starting point
%[ icol, ad, irow, ar] = setsp( A);  % set up data for Lanczos call
%ee=ones(n,1); % !!!Need better initial estimate for SMALLEST eig
%[ae, xnew] = mineig1( A, icol, ad, irow, ar, ee); % Lanczos version
%ae=eigs(A,1);    % get smallest eig of A
%ae=-min(0,ae); 
P=A+(norm(A,1)+1)*eye(n);  
P=sparse(P);
HA=H.*A;
tH2A=tH2.*A;
tH2P=tH2.*P;
L=tH2P-tH2A;   % initial dual variable
L=sparse(L);
[R,lmin]=chol(L);   % check pos def - dual feas.
i=0;step=1.1;
while lmin > 0,
  i=i+1;
  tH2Pn=(step^i)*tH2P;
  L=tH2Pn-tH2A; 
  L=sparse(L);
  [R,lmin]=chol(L);   % check pos def - dual feas.
end
P=sparse((step^i)*P);
gap=trace (P*L);
mu=(1/(2*n))*gap;
disp(['  ']);
disp(['START of ALGORITHM: ']);
fP=(H.*P)-HA;
optval=(norm( fP,'fro' ))^2;  
disp(['initial value of objective function is:  ',num2str(optval)]);
disp(['initial relative duality gap value:  ',num2str(full(gap/(optval+1)))]);
disp(['  ']);
iter=0;
LO=zeros(length(izs));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%MAIN WHILE LOOP%%%%%%%%%%%%%
while min((gap/(optval+1)),optval) > epsilon, 
  iter=iter+1;    %  iteration count
  disp(['   -  ']);
  disp(['iteration number:  ',num2str(iter)]);
%  solve for dual direction q
tLO=cputime;
etLO=clock;
  Li=L\speye(size(L));
%  Li=sparse(Li);   % no need - already sparse after inversion
  Li=(Li+Li)/2;
% mLi=sparse(mu*Li );     % this line used without redoing sparse!
  mLi=mu*Li;
   for s=1:length(izs),
   if s==2,
      tLO2=cputime;
   end
   LO(s,:)=(tH2(jzs(s),izs(s))*((mLi(jzs(s),izs)).*(Li(izs(s),jzs))));
   if s==2,
      disp(['cpu  time for s=2 for operator LO:  ',num2str(cputime-tLO2)]);
   end
%LO(s,:)=full(tH2(jzs(s),izs(s))*(full(mLi(jzs(s),izs)).*full(Li(izs(s),jzs))));
     LO(s,s)=LO(s,s)+1;
   end
  LO=sparse(LO);
disp(['cpu  time for operator LO:  ',num2str(cputime-tLO)]);
disp([' elapsed time for operator LO:  ',num2str(etime(clock,etLO))]);
  LOcond=condest(LO);
disp(['condition number for operator LO:  ',num2str(LOcond)]);
  % check dual feasibility
  tH2P=tH2.*P;
  errdual=sparse( (tH2P-tH2A)-L  );
  dnorm=normest(errdual);
  disp(['dual feasibility before updates:  ',num2str(dnorm)]);
  tqp= sparse(tH2.*sparse(mLi-P) );
  tqp=tqp+errdual;  % Always include dual feasibility fix
  tqp=tqp(indnzs);
tsolve=cputime;
etsolve=clock;
  if LOcond>=inf,   % singular matrix
     q=pinv(full(LO))*tqp;
     disp(['using pseudoinverse for singular operator LO ']);
  else
     q=LO\ tqp ;
  end
disp(['cpu time for system solve:  ',num2str(cputime-tsolve)]);
disp([' elapsed time for system solve: ',num2str(etime(clock,etsolve))]);
  errorsolve=normest(tqp-LO*q);
disp([' norm error in system solve: ',num2str(errorsolve)]);
  q2=zeros(n);
  q2(indnzs)=q;
  q=reshape(q2,n,n);
  q=sparse(q);
  q=( (q+q')/2 );   % ensure symmetry
  h=( mLi-mLi*q*Li-P );  % step for primal P
  h=(h+h')/2;

%%%%%%%%%%%%%%%%%%
%% check on Newton equation - testing only
%  normest( (tH2.*(mLi*q*Li) +q) - (tH2.*(mLi-A)-L)  )
%  normest( (tH2.*(h) -q) + (tH2.*(P-A)-L)  )
%   keyboard
%  end testing lines 
%%%%%%%%%%%%%%%%%%
%  safeguard positive definiteness - primal and dual
%  we use the Cholesky factor. to test pos. def.
    [R,lmind]=chol((L+q));   % check dual feasibility
    [R,lminp]=chol((P+h));   % check primal feasibility
    step=.5;
    i=0;
    qo=q;ho=h;
    while max(lmind,lminp) > 0,
        i=i+1;
        q=sparse( step^i*qo );
        h=sparse ( step^i*ho );
        [R,lminp]=chol(P+h);   % check primal feasibility
        [R,lmind]=chol(L+q);   % check dual feasibility
    end
    if i>0,   % Use the Newton step if possible
       q=sparse( (.97*step^i)*qo );
       h=sparse( (.97*step^i)*ho );
       disp(['using steplength:  ',num2str(.97*step^i)]);
    else
     disp(['Using a Newton step length 1 ']);
       q=qo;
       h=ho;
    end
%%%%%%%%
 %  update mu and P and L
    P=sparse(P+h);
    P=(P+P')/2;
    L=sparse(L+q);    
    L=(L+L')/2;
    gap=trace (P*L);
    mu=(1/(2*n))*gap;
% 
    fP=(H.*P)-HA;
    optval=(norm( fP,'fro' ))^2;  
    disp(['value of objective function after updates:  ',num2str(optval)]);
disp(['rel. dual. gap and gap values  after updates: '...
       ,num2str(full(gap/(optval+1))),'  ',num2str(gap)]);
    disp(['  ']);
end   
%%%%%%%%%%%%%%%%
% end of main while  loop 
% duality gap and dual feas must hold
disp(['  ']);
disp(['optimal P and optimal dual L found']);
optval=(norm( (H.*(P-A)),'fro' ))^2;  
disp(['value of objective function is:  ',num2str(optval)]);
disp(['  ']);
