code

Signed-off-by: 刘明宏 <liuminghong@mail.sdu.edu.cn>

bin/Solver/Pardiso/pardiso_MT_new.m

@@ -0,0 +1,63 @@
+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