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
