MATLAB实现矩阵雅可比（Jacobi）迭代

包含三个function，下列代码最后一段是主函数，其他都是function。三个function建议从下往上看。

function [x,k,resvec,DD,ID,JD,D,Ab] = jacobis(AA,IA,JA,b,x,tol,kmax)

%This function is an implementation of Jacobi’s
%iterative method for solving Ax=b.
%The input vector x is the initial guess, while tol
%is the tolerance used for stopping the iteration; kmax is the
%maximum number of iterations to be run. The output is the kth
%approximation to the exact solution, while resvec is a vector
%containing the 2-norms of the residual vectors.
n=length(b);
if ~exist(‘x’), x=zeros(n,1);end
if ~exist(‘tol’), tol=1e-6;end
if ~exist(‘kmax’), kmax=1e4;end

r=b-matvecs(AA,IA,JA,x);
res0=norm®;res=norm®;
resvec=res/res0;res=1;k=0;
%extract a vector D containing the main diagonal of matrix A
[DD,ID,JD]=diags(AA,IA,JA);
D=zeros(ID(end),1);
for i =1:length(DD)
D(ID(i))=D(ID(i))+DD(i);
end
while res>tol & k<kmax
k=k+1;
z=r./D;
x=x+z;
r=b-matvecs(AA,IA,JA,x);
res=norm®/res0;
resvec=[resvec res];
fprintf(’%2.9f\t %3d\n’,resvec(end),k)
end
%-----------------------------------------------------------
function w=matvecs(AA,IA,JA,v)

%Sparse matrix-vector product c=A*v: the matrix is provided
%in coordinate format AA, IA, JA, with m=IA(end), n=JA(end).
%The input vector v is assumed to be full.
%The output vector w is assumed to be full.

ma=IA(end);na=JA(end);nb=length(v);
if na~=nb,error(‘Wrong sizes: cannot perform matrix-vector product’);end

w=zeros(ma,1);
for k=1:length(AA)
w(IA(k))=w(IA(k))+AA(k)*v(JA(k));%edit this line to provide the correct expression for w
end
%------------------------------------------------------------
function [DD,ID,JD]=diags(AA,IA,JA)

%Jacobi splitting:
%Input: sparse matrix A in coordinate format (AA, IA, JA)
%Output: sparse diagonal matrix D in coordinate format (DD, ID, JD)
diag_values=[];
diag_row_index=[];
diag_col_index=[];
for k=1:length(AA)
if IA(k)==JA(k)
diag_row_index=[diag_row_index;IA(k)];
diag_col_index=[diag_col_index;JA(k)];
diag_values=[diag_values;AA(k)];
end
end
DD=diag_values;
ID=[diag_row_index;IA(end)];
JD=[diag_col_index;JA(end)];

n=100;
e = ones(n,1);
A = spdiags([-e 4*e -e], -1:1, n, n);
[IA,JA,AA]=find(A);
IA=[IA;n];JA=[JA;n];
[x,k,resvec] = jacobis(AA,IA,JA,e);