x3 = Q * (U \ (L \ P*b)) ; t3=toc ; nz3=nnz(L+U);

fprintf (<span class="string">'1: nnz(L+U): %5d time: %8.4f resid: %e\n'</span>, nz1,t1, norm(b-A*x1));
fprintf (<span class="string">'2: nnz(L+U): %5d time: %8.4f resid: %e\n'</span>, nz2,t2, norm(b-A*x2));
fprintf (<span class="string">'3: nnz(L+U): %5d time: %8.4f resid: %e\n'</span>, nz3,t3, norm(b-A*x3));
</pre><pre class="codeoutput">1: nnz(L+U): 16151 time:   0.0070 resid: 1.426954e-10
2: nnz(L+U): 15826 time:   0.0147 resid: 9.176936e-11
3: nnz(L+U)