% Subroutine of rcompletedPPA.m

function [x,lambda,iter,t] = rcompletedPPAinit(q,bx,m,n,i,j,mu,tol,maxiter)

% Initialization

tstart = tic;

fprintf(' Calculate the initial point by PADM2 \n')
meps = 1e-8; % machine epsilon
beta = min(i,j)*mu/(2*norm(bx,2)); 
taubeta = 1.61*beta;

% Initialization
freq = 10;

x = zeros(m,n*(i+j-1));
X = zeros(m*i,n*j);
lambda = zeros(m*i,n*j);
qx = q.*bx;
tq = 1 - q;

Hankeltime = 0;
HankelAdjtime = 0;
svd_time = 0;

Hankeltime0 = clock;
Hx = Hkx(x,m,n,i,j);
Hankeltime = Hankeltime + etime(clock,Hankeltime0);

fbest = inf;
xbest = x;

% Loop

iter = 0;

while iter < maxiter

  xold = x;
  Xold = X;
  lambdaold = lambda;
  
  G = lambda/beta - Hx;
  svd_time0 = clock;
  [R,S,T] = svd(G,'econ');
  svd_time = svd_time + etime(clock,svd_time0);
  
  X = R*max(S - mu/beta,0)*T'; % matrix shrinkage
  
  lX = lambda - beta*X;
  HankelAdjtime0 = clock;
  hxL = Hkxad(lX,m,n,i,j);
  HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);
  
  x = invmap(hxL + qx,beta,q,m,n,i,j);
  
  Hankeltime0 = clock;
  Hx = Hkx(x,m,n,i,j);
  Hankeltime = Hankeltime + etime(clock,Hankeltime0);
  
  lambda = lambda - taubeta*(Hx + X); 
  
  iter = iter + 1;
  
  if max([norm(x-xold,'fro'),norm(lambda-lambdaold,'fro'),norm(X-Xold,'fro')]) < meps % early termination
  
    svd_time0 = clock;
    [R,S,T] = svd(lambda,'econ');
    svd_time = svd_time + etime(clock,svd_time0);
    
    lambda2 = R*min(S,mu)*T';
    
%     HankelAdjtime0 = clock;
%     hLL = Hkxad(lambda2,m,n,i,j);
%     HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);
    xbest = x;
    
    fprintf(' early termination \n')
    break
    
  end 
  
  if mod(iter,freq) == 0
    
    fval = 1/2*norm(q.*x-qx,'fro')^2 + mu*sum(svd(Hx));
    
    svd_time0 = clock;
    [R,S,T] = svd(lambda,'econ');
    svd_time = svd_time + etime(clock,svd_time0);
    
    lambda2 = R*min(S,mu)*T';
    
    HankelAdjtime0 = clock;
    hLL = Hkxad(lambda2,m,n,i,j);
    HankelAdjtime = HankelAdjtime + etime(clock,HankelAdjtime0);
    
    dfeas = norm(tq.*hLL,'fro')/max([1,norm(hLL,'fro')]);
    dval = -0.5*norm(q.*hLL,'fro')^2 - mytrace(qx,hLL);
    
    if fbest > fval
      fbest = fval;
      xbest = x;
    end
    
%    fprintf(' iter, fval, dval, dfea, time = %g, %g, %g, %g, %g\n', iter, fval, dval, dfeas, toc)
    gap = abs(fbest-dval)/max(abs(dval),1);
    if max(gap,dfeas) < tol
      break
    end
    
  end
end

x = xbest;
lambda = lambda2;

t = toc(tstart);
fprintf('#iter=%g, time=%g, \n ', iter,t)
