%此函数为主函数，需要输入三个参数，维数n，初始点x0，构造H的辅助向量z，注意x0和z都是列向量；输出x和F。
%函数数据从文件data_g.m，data_y.m，data_z.m获取
%[x F] = ncp(2,[5 5]',[1 1]')在主程序窗口输入这样形式的数据即可。
function [x,F] = ncp(n,x0,z)
X = sym(zeros(1,n)); %为符号向量分配空间
max_iter = 5000; % 最大迭代步数
EPS = 1e-6; % 判断算法停止的误差界
%下面是参数
rho = 10; % rho > 0
p= 3; % p > 2
theta = 0.1; % 0 < theta < 1/2
k = 0; xk = x0; % 初始化
for i = 1:1:n
X(i) = sym(['x',num2str(i)]); %赋符号向量
end
[symF,symgrad_F] = Function(n); %调用函数 Function 得到符号函数symF及其梯度symgrad_F
while k < max_iter
F = symF; %由于每次迭代xk都会改变，这里需要每次都调用符号函数symF和symgead_F，然后再重新将两函数中的符号X用xk代替
grad_F = symgrad_F;
k = k + 1;
for j = 1:1:n %用数值代替符号得到函数F和grad_F
F = subs(F,X(j),xk(j));
grad_F = subs(grad_F,X(j),xk(j));
end
F = eval(F); %将F和grad_F恢复正常精度
grad_F = eval(grad_F);
[F_Phi,grad_Phi] = Phi(n,xk,F,grad_F); %调用Phi函数，得到函数F_Phi及其梯度gead_Phi
[F_Psi,grad_Psi] = Psi(F_Phi,grad_Phi); %调用Psi函数，得到函数F_Psi及其梯度gead_Psi
H = F_H( xk,n,F,grad_F,grad_Phi,z); %调用F_H函数，得到H
if norm( grad_Psi ) < EPS %算法停止条件判断
break;
else
if det( H ) == 0 %H奇异，取负梯度方向
dk = -grad_Psi;
else
dk = -inv(H)* F_Phi; %H非奇异，像这样取
if grad_Psi' * dk > -rho * (norm(dk))^p %下降条件判断
dk = -grad_Psi;
end
end
ik = get_k( n,xk,dk,symF,F_Psi,grad_Psi,theta); %调用函数get_k，得到最小ik
xk = xk + 2^(-ik) * dk; %更新迭代点
end
fprintf('This is %dth iteration\n', k);
xk
end
x = xk; %将最后结果赋给x输出
return;
