Equilibrate a matrix with a large condition number to improve the efficiency and stability of a linear system solution with the iterative solver gmres.

Load the west0479 matrix, which is a real-valued 479-by-479 sparse matrix. Use condest to calculate the estimated condition number of the matrix.

load west0479
A = west0479;
c1 = condest(A)

c1 = 1.4244e+12

Try to solve the linear system $\mathrm{Ax}=\mathit{b}$ using gmres with 450 iterations and a tolerance of 1e-11. Specify five outputs so that gmres returns the residual norms of the solution at each iteration (using ~ to suppress unneeded outputs). Plot the residual norms in a semilog plot. The plot shows that gmres is not able to achieve a reasonable residual norm, and so the calculated solution for $\mathit{x}$ is not reliable.

b = ones(size(A,1),1);
tol = 1e-11;
maxit = 450;
[x,flx,~,~,rvx] = gmres(A,b,[],tol,maxit);
semilogy(rvx)
title('Residual Norm at Each Iteration')

Use equilibrate to permute and rescale A. Create a new matrix B = R*P*A*C, which has a better condition number and diagonal entries of only 1 and -1.

[P,R,C] = equilibrate(A);
B = R*P*A*C;
c2 = condest(B)

c2 = 5.1036e+04

Using the outputs returned by equilibrate, you can reformulate the problem $\mathrm{Ax}=\mathit{b}$ into $\mathrm{By}=\mathit{d}$, where $\mathit{B}=\mathrm{RPAC}$ and $\mathit{d}=\mathrm{RPb}$. In this form you can recover the solution to the original system with $\mathit{x}=\mathrm{Cy}$.

Use gmres to solve $\mathrm{By}=\mathit{d}$ for $\mathit{y}$, and then replot the residual norms at each iteration. The plot shows that equilibrating the matrix improves the stability of the problem, with gmres converging to the desired tolerance of 1e-11 in fewer than 200 iterations.

d = R*P*b;
[y,fly,~,~,rvy] = gmres(B,d,[],tol,maxit);
hold on
semilogy(rvy)
legend('Original', 'Equilibrated', 'Location', 'southeast')
title('Relative Residual Norms (No Preconditioner)')
hold off

semilogy(rvy)
hold on

[L1,U1] = ilu(B,struct('type','ilutp','droptol',1e-1,'thresh',0));
[yp1,flyp1,~,~,rvyp1] = gmres(B,d,[],tol,maxit,L1,U1);
semilogy(rvyp1)

[L2,U2] = ilu(B,struct('type','ilutp','droptol',1e-2,'thresh',0));
[yp2,flyp2,~,~,rvyp2] = gmres(B,d,[],tol,maxit,L2,U2);
semilogy(rvyp2)

legend('No preconditioner', 'ILUTP(1e-1)', 'ILUTP(1e-2)')
title('Relative Residual Norms with ILU Preconditioner (Equilibrated)')
hold off