pcg.m

％ SSOR预处理的共轭梯度法求解方程Ax=b ％ 输入参数说明 ％ A 正定矩阵[n*n] ％ b 右边向量 ％ omega SSOR预处理参数(0--2) ％ Times 迭代次数 ％

function [NewX,avgerr]=pcg(A,b,omiga,Times,errtol)
% SSOR预处理的共轭梯度法求解方程Ax=b
% 输入参数说明
%  A     正定矩阵[n*n]
%  b     右边向量
%  omega SSOR预处理参数(0--2)
%  Times 迭代次数
%  errtol 给定误差终止条件
%
%输出参数
%  NewX  方程Ax=b的x近似解
%  avgerr 求解的当前平均绝对误差

%omiga=1;
%tic;   %计时器开始，与计算无关
[cols,rows]=size(b);
if cols > rows
   x=ones(cols,1);% x的初始值
else 
   x=ones(rows,1);% x的初始值
   b=b';
end;
[D,L,U]=dlu(A);% 对A进行分解子函数，见下面
M=(D-omiga*L)*inv(D)*(D-omiga*U)/(omiga*(2-omiga));%计算ssor预处理矩阵M

%下面为pcg方法的计算过程，完全翻译，所以不做标注
r=b-A*x;
invm=inv(M);
z=invm*r;
p=z;
x0=x-1;
k=1;
Error = x-x0; 
avgerr = 99999;
while (k<=Times)&&(max(abs(Error))>errtol)
    zr=z'*r;
    alpha=zr/(p'*(A*p));
    x0=x;
    x=x+alpha*p;
    r=r-alpha*A*p;
    z=invm*r;
    beta=(z'*r)/zr;
    p=z+beta*p;
    Error = x-x0; 
    ErrR=b-A*x;
    MSEError=sum(abs(ErrR))/length(x); % 求平均绝对误差
    if MSEError < avgerr
        avgerr = MSEError;
        NewX=x;
    end;
    k=k+1;
end
%maxerr=max(abs(x-x0))
%minerr=min(abs(x-x0))
%avgerr=sum(abs(OldError))/length(NewX);
%disp(['迭代次数k=',num2str(k)]);%先是迭代次数k和耗时
%disp(['pcg method 耗时(s) :',num2str(toc)]);

 

function [D,L,U]=dlu(A) %dlu分解子函数
D=zeros(size(A));
L=D;
U=D;
for i=1:size(A,1)
    D(i,i)=A(i,i);
end
for i=1:size(A,1)
    for j=1:size(A,2)
        L(i,j)=-(i<j)*A(i,j);
        U(i,j)=-(i>j)*A(i,j);
    end
end;