The version in decomposition does some optional steps of iterative refinement: It uses the same solution you have above based on calling lu (up to minor round-off because things are stored in different format). However, it then checks the resulting residual and does another solve with just that residual as a right-hand side if necessary:
rng(42)
DeltaSize = 500;
M = DeltaSize:DeltaSize:2000; % Reduced length for run-time
ErrorLU = zeros(numel(M),1);
ErrorLUrefined = zeros(numel(M),1);
ErrorLUObject = zeros(numel(M),1);
for cnt = 1:numel(M)
disp(num2str(cnt))
MatrixDimension = M(cnt);
NumberOfNNZ = 0.1 * MatrixDimension^2;
% Build S and B
row = randi([1,MatrixDimension],NumberOfNNZ,1);
col = randi([1,MatrixDimension],NumberOfNNZ,1);
v = rand(NumberOfNNZ,1) + 1i * rand(NumberOfNNZ,1);
S = sparse(row,col,v,MatrixDimension,MatrixDimension) + speye(MatrixDimension);
B = rand(MatrixDimension,1) + 1i * rand(MatrixDimension,1);
% use LU
[l,u,p,q,d] = lu(S);
X = q * (u \ (l \ (p * (d \ B))));
ErrorLU(cnt) = norm(full(B - S * X));
% use LU with iterative refinement
R = B - S * X;
Xerr = q * (u \ (l \ (p * (d \ R))));
X = X + Xerr;
ErrorLUrefined(cnt) = norm(full(B - S * X));
% use decomposition object
dS = decomposition(S,'lu');
X = dS \ B;
ErrorLUObject(cnt) = norm(full(B - S * X));
end
LogErrorLU = log10(ErrorLU);
LogErrorLUObject = log10(ErrorLUObject);
figure(1)
semilogy(M,ErrorLU,'.-','Displayname','lu(S)')
hold on
semilogy(M,ErrorLUrefined,'x-','Displayname','lu(S) iterative refinement')
semilogy(M,ErrorLUObject,'.-','Displayname','decomposition(S,lu)')
hold off
ylim([1e-17 1e-5])
legend show
grid on;
ax = gca;
ax.YMinorGrid = 'off';
This is because decomposition does its best to match exactly what S\B does.
