sirt.m

这是在网上下的一个东东

%
% 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.
%      maxiter  Stop after maxiter iterations.
%
% Outputs:
%      x        Solution.
%
function x=sirt(A,b,tolx,maxiter)
%
%
%
alpha=1.0;
%
% First, get the size of A.
%
[m,n]=size(A);
%
%
%
alpha=1.0;
%
% In the A1 array, we convert all nonzero entries in A to +1.  
%
A1=(A>0);
%
% 
%
AP=A';
A1P=A1';
%
% 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);
NRAYS=zeros(n,1);
for i=1:m
  N(i)=sum(A1P(:,i));
  L(i)=sum(AP(:,i));
end
for i=1:n
  NRAYS(i)=sum(A1(:,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=A1P(:,i)'*newx;
%
% We use the following more accurate formula for delta.
%
    delta=b(i)/L(i)-q/N(i);
%
% This formula is less accurate and doesn't work nearly as well.
%
%    delta=(b(i)-q)/N(i);
%
%
%  Perform updates for those cells touched by ray i.
%
    deltax=deltax+delta*A1P(:,i);
  end
%
%  Now, add the average of the updates to the old model.
%  Note the "./" here.  This means that the averaging is done with
%  respect to the number of rays that pass through a particular cell!
%
  newx=newx+alpha*deltax./NRAYS;
%
% Check for convergence
%
  if (norm(newx-x)/(1+norm(x)) < tolx)
    return;
  end
%
% No convergence, so setup for the next major iteration.
%
  x=newx;
end