function [x,info,res] = pardiso_MT_new(A,b,OMP)
% Solve for x in linear system A * x = b.
% --------------------------------------
% Initialize the PARDISO internal data structures. We've told PARDISO to
% handle complex matrices using the sparse direct solver.

tp = tic;
fprintf(['Degree of freedom： ' num2str(size(A,1)) '\n']);
fprintf ('Solution to Ax=b via PARDISO\n') ;
trilA = tril(A);
% eps=1e-8;
% A = A + eps * speye(size(A));

% info = pardisoinit(MTYPE,SOLVER);
%SOLVER：0 sparse direct solver; 1 multi-recursive iterative solver
info = pardisoinit(6,0);

% OMP_NUM_THREADS is honored
info.iparm(3) = OMP;
info.iparm(28) =1;
% Analyze the matrix and compute a symbolic factorization.
verbose = false;%ture
info = pardisoreorder(trilA,info,verbose);
% fprintf('The factors have %d nonzero entries.\n',info.iparm(18));

% Compute the numeric factorization.
t = tic;
info = pardisofactor(trilA,info,verbose);
time_fac = toc(t);
% fprintf('The time of numeric factorization %.2f s.\n',time_fac);

% Compute the solution x using the symbolic factorization.
info.iparm(25)=0 ; %0 indicates that a sequential forward and backward solve is used
t = tic;
[x,info] = pardisosolve(trilA,b,info,verbose);

% x = zeros(size(b));
% for i = 1:size(b,2)
%   [x(:,i),info] = pardisosolve(tril(A),b(:,i),info,verbose);
% end

time_sol = toc(t);
% fprintf('The time of solution %.2f s.\n',time_sol);
% fprintf('Memory factorization and solution %.2f GB.\n', single(info.iparm(17))/1024/1024);
fprintf('Total peak memory consumption is peakmem = %.2f GB.\n', ...
    single( max(info.iparm(15), info.iparm(16)+info.iparm(17)) )/1024/1024);

% Compute the residuals.
res=zeros(2,1);
for i = 1:2%size(b,2)
    res(i) = norm(A*x(:,i)-b(:,i))/norm(b(:,i));
    fprintf('%s %d\n',['Normalized residual for rhs' num2str(i) ' is'], res(i));
end

% Free the PARDISO data structures.
if nargout<2
    pardisofree(info);
    clear info
end

Solve_time = toc(tp);
fprintf('The time of solve Ax=b %.2f s.\n',Solve_time);
end