43．bicg函数

在MATLAB中，提供了bicg函数用于实现双共轭梯度法。函数的语法格式为：

x＝bicg（A，b）：针对x对线性方程组A∗x＝b求解。n×n系数矩阵A必须是方阵，并且应为大型稀疏矩阵。列向量b的长度必须为n。A可以是函数句柄afun，这样，afun（x，'notransp'）返回A∗x，afun（x，'transp'）返回A'∗x。

如果bicg收敛，则会显示一条关于该结果的消息。如果bicg无法在达到最大迭代次数后收敛或出于任何原因暂停，则会输出一条包含相对残差norm（b－A∗x）/norm（b）以及该方法停止或失败时所达到的迭代数的警告消息。

bicg（A，b，tol）：指定该方法的容差。如果tol为[]，bicg使用默认值1e－6。

bicg（A，b，tol，maxit）：指定最大迭代次数。如果maxit为[]，bicg使用默认值min（n，20）。

bicg（A，b，tol，maxit，M）和bicg（A，b，tol，maxit，M1，M2）：使用预设子条件M或M＝M1∗M2，并高效求解关于x的方程组inv（M）∗A∗x＝inv（M）∗b。如果M为[]，bicg不会应用预设子条件。M可以是函数句柄mfun，这样，mfun（x，'notransp'）返回M\x，mfun（x，'transp'）返回M'\x。

bicg（A，b，tol，maxit，M1，M2，x0）：指定初始估计值。如果x0为[]，bicg使用默认值（即全部为零的向量）。

[x，flag]＝bicg（A，b，…）：也返回一个收敛标志flag，其取值见表1-1。

提示：如果flag不为0，返回的解x具有在所有迭代中最小的范数残差。如果指定flag输出，则不会显示消息。

[x，flag，relres]＝bicg（A，b，…）：如果flag为0，relres≤tol，还返回相对残差norm（b－A∗x）/norm（b）。

[x，flag，relres，iter]＝bicg（A，b，…）：其中0≤iter≤maxit，还返回计算x时所达到的迭代数。

[x，flag，relres，iter，resvec]＝bicg（A，b，…）：还返回每次迭代中的残差范数的向量（包括norm（b－A∗x0））。

【例1-45】实例演示如何使用预设子条件。

     %加载A = west0479,它是一个非对称的479×479实稀疏矩阵
     load west0479;
     A = west0479;
     %定义b以使实际解是全为1的向量
     b = full(sum(A,2));
     %设置容差和最大迭代次数
     tol = 1e-12;
     maxit = 20;
     %使用bicg根据请求的容差和迭代次数求解
     [x0,fl0,rr0,it0,rv0] = bicg(A,b,tol,maxit);

bicg未在请求的20次迭代内收敛至请求的容差1e－12，因此fl0为1。实际上，bicg的行为太差，因此初始估计值（x0＝zeros（size（A，2），1））是最佳解，并会返回最佳解（如it0＝0所示）。MATLAB将残差历史记录存储在rv0中。

     %绘制bicg的行为
     semilogy(0:maxit,rv0/norm(b),'-o');
     xlabel('迭代数');
     ylabel('相对残差');

运行程序，效果如图1-18所示。

图1-18表明该解未收敛，可以使用预设子条件改进结果，将在例1-47中讲解。