From octave-sources-request at bevo dot che dot wisc dot edu  Mon Aug 14 03:58:02 2000
Subject: Toeplitz and Hankel solvers and inverses
From: Dirk Laurie <dirk at calvyn dot puk dot ac dot za>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Mon, 14 Aug 2000 10:58:04 +0200


--mYCpIKhGyMATD0i+
Content-Type: text/plain; charset=us-ascii

I have coded the standard Durbin and Trench algorithms for
solving Toeplitz systems and inverting Toeplitz matrices
respectively, but generalized to the unsymmetric case.

The calling protocol provides for the symmetric case,
but the code does not actually exploit it (the savings
would be only 33%, and this routine is not going to be
really fast unless someone makes an .oct version of it.)

The corresponding problems involving Hankel matrices are
handled by the obvious row and column reversal.

Dirk Laurie


--mYCpIKhGyMATD0i+
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="solvtoep.m"

function [x,y,w]=solvtoep(c,r,b)
% solvtoep(c,r,b) == toeplitz(c,r)\b
% solvtoep(c,b) == toeplitz(c)\b
% [x,y,w]=solvtoep(...): 
%    x is solution as above
%    y is last column of inv(T) 
%    w' is last row of inv(T)
% If only y and w needed:
% [y,w]=solvtoep(c,r,[])
% y=solvtoep(c,[]) 

% Dirk Laurie 2000-08-11
% See Golub and Van Loan 3rd edition Section 4.7
% This code does not exploit symmetry if present.  You can
% make it 33% faster in the symmetric case by identifying w with y.
if nargin<3, b=r; r=c; end
if r(1)!=c(1), 
  disp('warning: solvtoep: column wins diagonal conflict');
end
n=length(c); dox=length(b)==n;
t=c(1); c(1)=[]; r(1)=[]; c=c(:); r=r(:);  
den=t; if dox, b=b(:); x=b(1)/den; end
y=zeros(n-1,1); w=y;
if n<2, return, end
w(1)=-c(1)/t; y(1)=-r(1)/t;
for k=1:n-1,
  j=1:k; l=k:-1:1;
  den=(1-w(k)*y(k))*den;
  if dox, 
    x(k+1)=(b(k+1)+w(j)'*b(l))/den;  
    x(j)=x(j)+x(k+1)*y(l);
  end
  if k<n-1,
    wl=w(l);
    w(k+1)=(-c(k+1)-w(j)'*c(l))/den;
    w(j)=w(j)+w(k+1)*y(l);
    y(k+1)=(-r(k+1)-y(j)'*r(l))/den;
    y(j)=y(j)+y(k+1)*wl;
  end
end 
if nargout>1 || ~dox,  
  y=[y(l)/den; 1/den]; w=[w(l)/den; 1/den]; end
if ~dox, x=y; y=w; clear w; end

--mYCpIKhGyMATD0i+
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="invtoep.m"

function A=invtoep(c,r)
% invtoep(c,r) == inv(toeplitz(c,r))

% Dirk Laurie 2000-08-11
% See notes to solvtoep.

if nargin==1, r=c; end
[y,w]=solvtoep(c,r,[]);
n=length(c); A=zeros(n);
A(:,n)=y; A(n,:)=w'; d=y(n);
i=1; ni=n-i;
while i<=ni,
  j=i:ni;
  A(i,j)=A(n+1-j,n+1-i)';
  A(j,i)=A(n+1-i,n+1-j)';
  if j==i, return; end
  j=i:ni-1;
  A(n-j,n-i)=A(i,j)'-(w(j)*y(i)-w(n-i)*y(n-j))/d;
  A(n-i,n-j)=A(j,i)'-(y(j)*w(i)-y(n-i)*w(n-j))'/d;
  i=i+1; ni=ni-1;
end


--mYCpIKhGyMATD0i+
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="solvhank.m"

function x=solvhank(c,r,b)
% solvhank(c,r,b) == hankel(c,r)\b
% solvhank(c,b) == hankel(c)\b

% Dirk Laurie 2000-08-11
n=size(c); k=n:-1:1;
if nargin<3, b=r; r=c(k); end
x=solvtoep(c(k),r,b(k));
 

--mYCpIKhGyMATD0i+
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="invhank.m"

function A=invtoep(c,r)
% invhank(c,r) == inv(hankel(c,r))
 
% Dirk Laurie 2000-08-11 
n=size(c); k=n:-1:1;
if nargin<2, r=c(k); end
A=invtoep(c(k),r); A(:,k)=A;

--mYCpIKhGyMATD0i+--



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------