1、使用稀疏矩阵定义线性方程组
当您使用分布式函数时，MATLAB会使用默认的集群设置自动启动一个并行池。这个例子使用了来自MATLAB gallery函数的Wathen矩阵。该矩阵是一个稀疏的、对称的、总体维数为
的随机矩阵。
n = 400;
A = distributed(gallery('wathen',n,n));
使用“local”配置文件启动并行池(parpool)...
连接到并行池 (number of workers: 6).
N = 481601
在此例中,b被定义为A的行和, 从Ax=b得到 x的精确解
b = sum(A,2);
因为sum作用于一个分布式数组，所以b也是分布式的，它的数据存储在并行池的worker的内存中。最后，您可以定义精确的解，以便与使用迭代方法得到的解进行比较。
xExact = ones(N,1,'distributed');
2、用共轭梯度法解线性方程组
[xCG_1,flagCG_1,relres_CG1,iterCG_1,resvecCG_1] = pcg(A,b);
errCG_1 = abs(xExact-xCG_1);
figure(1)
hold off
semilogy(errCG_1,'o');
title('System of Linear Equations with Sparse Matrix');
ylabel('Absolute Error');
xlabel('Element in x');
计算结果误差较大。
在此示例中，解决方案未收敛于默认的最大迭代次数内，这会导致较高的误差。
为了增加收敛的可能性，您可以自定义公差和最大迭代步骤数的设置。
tolerance = 1e-12;
maxit = N;
tCG = tic;
[xCG_2,flagCG_2,relresCG_2,iterCG_2,resvecCG_2] = pcg(A,b,tolerance,maxit);
tCG = toc(tCG);
flagCG_2
errCG_2 = abs(xExact-xCG_2);
figure(2)
hold off
semilogy(errCG_1,'o');
hold on
semilogy(errCG_2,'d');
title('Comparison of Absolute Error');
ylabel('Absolute Error');
xlabel('Element in x');
legend('Default tolerance and iterations','Improved tolerance and iterations');
hold off
与以前的解决方案相比，这个解决方案有一个改进的绝对误差。
但显然，默认步骤数不足以实现该系统的良好解决方案。
2、用预处理共轭梯度法求解线性方程组
使用预条件共轭梯度(PCG)方法可以提高解系统的效率。
首先，使用预处理矩阵对线性方程组M进行预处理。接下来，使用CG方法求解你的预处理方程组。
PCG方法比CG方法的迭代次数少得多。
pcg方法也使用MATLAB函数pcg。你可以提供一个合适的预处理矩阵作为额外的输入。理想的预调节器矩阵是其逆M−1近似于系数矩阵的逆a−1，但更容易计算的矩阵。这个例子使用A的对角线来确定线性方程组的前提条件。
M = spdiags(spdiags(A,0),0,N,N);
tPCG = tic;
[xPCG,flagPCG,relresPCG,iterPCG,resvecPCG]=pcg(A,b,tolerance,maxit,M);
tPCG = toc(tPCG);
figure(4)
hold off;
semilogy(resvecCG_2./resvecCG_2(1))
hold on;
semilogy(resvecPCG./resvecPCG(1))
title('Evolution of Relative Residual');
ylabel('Relative Residual');
xlabel('Iteration Step');
legend('Residuals of CG','Residuals of PCG with M \approx diag(A)')
PCG方法与非预处理系统相比，需要更少的步骤来收敛。这个结果也反映在执行时间上。
fprintf([...
'\nTime to solve system with CG: %d s', ...
'\nTime to solve system with PCG: %d s'],tCG,tPCG);
Time to solve system with CG: 1.244593e+01 s
Time to solve system with PCG: 7.657432e-01 s
除了用更少的迭代步骤来解决这个示例系统之外，PCG方法还可以返回更精确的解。
errPCG = abs(xExact-xPCG);
figure(5)
hold off
semilogy(errCG_1,'o');
hold on
semilogy(errCG_2,'d');
semilogy(errPCG,'x');
title('Comparison of absolute error');
ylabel('Absolute error');
xlabel('Element in x');
legend('CG default','CG custom','PCG');
在完成计算之后，可以删除并行池。gcp函数返回当前并行池对象，以便您可以删除当前池。
delete(gcp('nocreate'))
