levinson_algorithm.m

列文松法解对称TOEPLITZ矩阵组成的线性方程组

% This program is used to solve the Toeplitz matrix using the Levinson
% algorithm.  Preference:"数值分析与算法",page 118.  4/16/08.
function X=levinson_algorithm(A,B)

% A=[1 2 3 4 5 6;
%    2 1 2 3 4 5;
%    3 2 1 2 3 4;
%    4 3 2 1 2 3;
%    5 4 3 2 1 2;
%    6 5 4 3 2 1];
% h=[12 3 54 67 45 41]';
% B=A*h;

[N,N]=size(A);
x=zeros(N,N);
y=zeros(N,N);
a=zeros(1,N+1);
q=zeros(1,N);
b=zeros(1,N);
c=zeros(1,N);
w=zeros(1,N);

t(1:N)=A(1,:);
a(1)=t(1);
y(1,1)=1;
x(1,1)=B(1)/t(1);              % All computation  Mul   2N^2+6N-8
for k=1:N-1                    %                  Add   2N^2+6N-8         
    p1=0;
    p2=0;
    for j=1:k                         
        p1=p1+t(k-j+2)*x(k,j);      % k Mul   k  Add 
        p2=p2+t(j+1)*y(k,j);        % k Mul   k Add   
    end
    q(k)=p1;
    b(k)=p2;
    c(k)=-b(k)/a(k);                % 1 Mul
    
    y(k+1,1)=c(k)*y(k,k);           % 1 Mul
    if k>=2
        for i=2:k
            y(k+1,i)=y(k,i-1)+c(k)*y(k,k-i+1);   % k-1 Mul  K-1 Add
        end
    end
    y(k+1,k+1)=y(k,k);
    
    a(k+1)=a(k)+c(k)*b(k);           % 1 Mul  1 Add
    w(k)=(B(k+1)-q(k))/a(k+1);         %  1 Mul  1 Add
    for i=1:k
        x(k+1,i)=x(k,i)+w(k)*y(k+1,i);       % k MUl  k  Add
    end
    x(k+1,k+1)=w(k)*y(k+1,k+1);       % 1 Mul
end

% disp('X=');disp(x);

X=x(N,:)';