art.m

来自「这是在网上下的一个东东」· M 代码 · 共 110 行

%
% x=art(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=art(A,b,tolx,maxiter)
%
% Alpha is a damping factor.  If alpha<1, then we won't take full steps
% in the ART direction.  Using a smaller value of alpha (say alpha=.75)
% can help with convergence on some problems.  
%
alpha=1.0;
%
% 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);
%
% Get transposed copies of the arrays for faster access.
%
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);
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.');
    x=newx;
    return;
  end
%
% Start the next round of updates with the current solution.
%
  newx=x;
%
% Now, update each of the m rays.
%
  for i=1:m
%
%  Compute the approximate travel time for ray i.
%
    q=A1P(:,i)'*newx;
%
% We use the more accurate formula for delta.
% 
    delta=b(i)/L(i)-q/N(i);
%
% This alternative formula is less accurate and doesn't work nearly as well.
%
%    delta=(b(i)-q)/N(i);
%
%
% Now do the update.
%
    newx=newx+alpha*delta*A1P(:,i);
  end
%
% Check for convergence
%
  if (norm(newx-x)/(1+norm(x)) < tolx)
    x=newx;
    return;
  end
%
% No convergence, so setup for the next major iteration.
%
  x=newx;
end