稳定的双共轭梯度法matlab程序

以下是一个稳定的双共轭梯度法的Matlab程序，其中使用了预条件：

function [x, flag, relres, iter, resvec] = stable_bicg(A, b, tol, maxit, M)

% Inputs:
% A - 系数矩阵
% b - 右手边向量
% tol - 相对残差容限
% maxit - 迭代次数上限
% M - 预条件矩阵

% Outputs:
% x - 解向量
% flag - 表示迭代是否成功收敛
% relres - 相对残差
% iter - 实际迭代次数
% resvec - 残差向量

% 确保 A 和 b 是列向量
b = b(:);
n = length(b);

% 初始化变量
x = zeros(n, 1);
r = b;
rtilde = r;
p = zeros(n, 1);
ptilde = zeros(n, 1);
resvec = zeros(maxit, 1);
flag = 0;

% 计算初始残差和范数
normb = norm(b);
normr = norm(r);
if (normr <= tol*normb)
    flag = 0;
    relres = normr/normb;
    iter = 0;
    resvec(1) = relres;
    return;
end

for iter = 1:maxit
    rho = rtilde'*r;
    if (rho == 0)
        flag = 4;
        break;
    end
    
    if (iter == 1)
        p = r;
        ptilde = rtilde;
    else
        beta = rho/rhoold;
        p = r + beta*p;
        ptilde = rtilde + beta*ptilde;
    end
    
    z = M\p;
    ztilde = M'\ptilde;
    q = A*z;
    qtilde = A'*ztilde;
    
    alpha = rho/(qtilde'*p);
    x = x + alpha*z;
    r = r - alpha*q;
    rtilde = rtilde - alpha*qtilde;
    
    normr = norm(r);
    resvec(iter) = normr/normb;
    if (normr <= tol*normb)
        flag = 0;
        break;
    end
    
    rhoold = rho;
end

if (flag == 0)
    relres = normr/normb;
else
    relres = NaN;
end
resvec(iter+1:maxit) = [];