function [s,x]=sirt(A,b,tolx,maxiter)
% Syntax: x=sirt(A,b,tolx,maxiter)
%
% Inputs:
% A = Constraint matrix.约束矩阵
% b = right hand side.
% tolx = Terminate when the relative diff between
% two iterations is less than tolx.
% 当两个迭代之间的相对差异小于tolx时终止。
% maxiter = Stop after maxiter iterations.
% 麦克斯特后停止迭代
% Outputs:
% x = Solution.结果
% cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
% First, get the size of A.
[m,n]=size(A);
%
% In the A1 array, we convert all nonzero entries in A to +1.
A1=A>0;
%
% Start with the zero solution.
x=zeros(n,1);
%
% Start the iteration count at 0.
iter=0;
%
% Precompute N(i) and L(i) factors.
N=zeros(m,1);
L=zeros(m,1);
for i=1:m
N(i)=sum(A1(i,:));
L(i)=sum(A(i,:));
end;
%
% Now, the main loop.
while (1==1)
%
% Check to make sure that we haven?t exceeded maxiter.
iter=iter+1;
if (iter > maxiter)
disp('Max Iterations Exceeded.')
return;
end;
%
% Start the next round of updates with the current solution.
newx=x;
%
% Now, compute the updates for all of the rays and
% all cells, and put them in a vector called deltax.
%
deltax=zeros(n,1);
for i=1:m
%
% Compute the approximate travel time for ray i.
q=A1(i,:)*newx;
%
% Perform updates for those cells touched by ray i.
for j=1:n
if (A(i,j)>0)
%
% We can use either of our two approximate formulas for the update.
% deltax(j)=deltax(j)+(b(i)-q)/N(i);
deltax(j)=deltax(j)+b(i)/L(i)-q/N(i);
end;
end;
end;
%
% Now, add the average of the updates to the old model.
%
newx=newx+deltax/m;
%
% Check for convergence
%
if (norm(newx-x)/(1+norm(x)) < tolx)
return;
end;
%
% No convergence, so setup for the next major iteration.
%
x=newx;
s=post2grid(x);
%
end;
% cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function g=post2grid(p)
g=zeros(4,4);
% Post to Grid
for i=1:4
for j=1:4
ndx=4*(i-1)+j;
g(i,j)=p(ndx);
end;
end;
评论0