function [X, val_obj] = Correlation_Newton(G,C,m,n,tau,tau0,tol)

%%%%%%%%%%%%%%%%%%%                           %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% This code is for the Local Correlation Stress Testing of the  Lagrangian Dual approach described in
%%%%%%  Houduo Qi and  Defeng Sun, "Correlation Stress Testing for Value-at-Risk: An Unconstrained
%%%%%%  Convex Optimization Approach", Department of Mathematics, National
%%%%%%  University of Singapore, March 2007 
%%%%%%%%%%%%%%%%%%%%%%%                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% 
% for solving 
%%  min  0.5 || X -  C ||^2 
%%  s.t.  X_ii = 1, i=1, ..., n
%%        X_ij = C_ij, i=1, ..., m; j =1, ..., n (Recall X = X^T)
%%        X- tau*I >=0 (symmetric and positive semi-definite)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% %%%%%%%%%%%%%%%%%%%%%%%
% G: the unstresed correlation matrix (n by n)
% C: stressed correlation matrix (n by n)
% C(1:m,1:m) is fixed
% tau0: the threshold for positive definiteness of  C  
% tau: the lower bound for the smallest eigevalue of X 
% tol: stopping crierion for the KKT system: 1.0e-5 ~ 1.0e-6
% X: the calibrated correlation matrix
% val_obj: final objective function value

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Last updated on March 20, 2007 %%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ---Newton method starts--- ')
t0=cputime;


G =(G+G')/2; % make G symmetric
C =(C+C')/2;

if m<1
    tau =tau0;
else
    tau =min(tau, 0.90*tau0);
end

n_2=n-m;

if m>0
    rhs1 =ones(n,1);
    [P,D] = eig(G);
    d = diag(D);
    d =real(d);
    if min(d) < tau0
        fprintf('\n')
        fprintf('The smallest eignevalue of the initial G < %d \n', tau0)
        disp(' ---Start the Pre-processing--- ')
        [G,y]=CorNewton1(G,rhs1,tau0);
        C(1:m,1:n) = G(1:m,1:n); % update G to make sure  G is positive definite
        C((m+1):m,1:m) =  G((m+1):m,1:m);
        disp(' ---Pre-processing finished--- ')
    else
        fprintf('\n')
        fprintf('The smallest eignevalue of the initial matrix G(1:m,1:m) >= %d \n', tau0)
        disp(' ---No pre-processing needed--- ')
    end
end

C =C-tau*eye(n);  % reset C: C>= tau *I;

C_1=zeros(m, m);
C_2=zeros(m, n_2);

C_1=C(1:m, 1:m);
C_2=C(1:m, (m+1):n);

% Calculate Schur Complement
Y=C_2'*inv(C_1)*C_2;
dd=diag(Y);

C_3=zeros(n_2, n_2);
C_3=C((m+1):n, (m+1):n);


b=(1-tau)*ones(n_2,1)-dd;

C_3 = C_3-Y; % C_3 updated 
  
[Y0,y0] = CorNewton1(C_3,b,0); % updated C_3 is calibrated

C_3 = Y0+Y; % C_3 is updated back

X = [C_1,C_2; C_2',C_3]; 
G = X-C;
val_obj = sum(sum(G.*G))/2;
X = X + tau*eye(n); % optimal X


fprintf('Newton: Final function value %d \n', val_obj)

time_used = cputime -t0