function program_cr; 
%  program cr implements algorithm CR described in
%  P. Favati, G. Lotti, O. Menchi
%  "Recursive algorithms for unbalanced banded Toeplitz systems"
%  it solves a system whose matrix and right hand-side are 
%  constructed by function problem_data (an example is given in
%  a separate file). At the end the computed solution is found in x.
%
global m n lnn 
global row col 
global rhs sol x
global h0 v0 hcap t r v1
global fv0 fr ft fv1
global frti fv0ti dv1ti fhti
%
%  function problem_data provides the first row and column of a 
%  Toeplitz matrix with m upper diagonals, n lower diagonals and
%  2^lnn blocks. In a simulation context (as in the example) the 
%  right hand-side is constructed  from a given solution sol.  
%
problem_data;
m, n, lnn
%
nn=2^lnn; dim= n*nn;
%
transformation_cr; 
substitution_cr;
%
%  In a simulation context the computed solution x is compared
%  with the exact solution sol. Drop the following line if sol
%  is not available.
%
error=max(abs(sol-reshape(x(1:nn,1:n,1)',1,dim)))
%

%
% ------------------------------------------------------------------
%

%
function transformation_cr;
%
%    computes the d-vectors and the corresponding FFTs
%    at all levels
%
global m n lnn 
global row col 
global h0 v0 hcap t r v1
global fv0 fr ft fv1
global frti fv0ti dv1ti fhti
%
%       constants
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
an=col(n+1);
lnn1=lnn+1;
%
%        utility vectors
%
h0=zeros(lnn1,n);  v0=zeros(lnn1,n);  hcap=zeros(lnn1,n); 
t=zeros(lnn1,n);  r=zeros(lnn1,m);  v1=zeros(lnn1,m);
%
v0ti=zeros(lnn1,n); hti=zeros(lnn1,n); rti=zeros(lnn1,n);
v1ti=zeros(lnn1,n);
%
fv0=zeros(lnn1,n2); ft=zeros(lnn1,n2); fv1=zeros(lnn1,m2);
fr=zeros(lnn1,m2);
%
frti=zeros(lnn1,n2); fv0ti=zeros(lnn1,n2); dv1ti=zeros(lnn1,n2);
fhti=zeros(lnn1,n2);
%
%                 computes the d-vectors at first level 
%
h0(1,:)=col(1:n);      
v0(1,:)=eye(1,n);
hcap(1,:)=col(1:n);
t(1,:)=col(1:n);
r(1,:)=fliplr(row(2:m+1));
v1(1,:)=zeros(1,m);
%
for i=1:lnn
    %
    %    computes the d-vectors of H0^-1 at current level
    %
    z=inverh0(i);                  
    rti(i,:)=z(:,1)';
    v0ti(i,:)=z(:,2)';
    v1ti(i,:)=z(:,3)';
    hti(i,:)=z(:,4)';
    %
    %                    computes the corresponding FFTs
    %
    fv0(i,:)=ifft(v0(i,:),n2)*sn2;    
    fr(i,:)=ifft(r(i,:),m2)*sm2;
    ft(i,:)=ifft(t(i,:),n2)*sn2;
    fv1(i,:)=ifft(v1(i,:),m2)*sm2;
    frti(i,:)=ifft(rti(i,:),n2)*sn2;
    fv0ti(i,:)=ifft(v0ti(i,:),n2)*sn2;
      te=[zeros(1,n+1),v1ti(i,2:n)];
    dv1ti(i,:)=ifft(te)*sn2;
    fhti(i,:)=ifft(hti(i,:),n2)*sn2;
    %
    %           computes the d-vectors for the next level
    %
    parr={fr(i,:),fv1(i,:)};
    parh1={ft(i,:),fv0(i,:)};
     %
      vetnt=v0ti(i,:);
      vet=fv0ti(i,:);
    v0(i+1,:)=-molth1(an,parh1,vetnt,vet);
      vetnt=rti(i,:);
      vet=frti(i,:);
      z1=molth1(an,parh1,vetnt,vet);
      vetnt=v0ti(i,1:m);
      vet=ifft(vetnt,m2)*sm2;
      z2=moltr(parr,vetnt,vet);
      vetnt=rti(i,1:m);
      vet=ifft(vetnt,m2)*sm2;
    r(i+1,:)=-moltr(parr,vetnt,vet);
      te=[zeros(1,n-m),z2];
    hcap(i+1,:)=hcap(i,:)-an*te;
    t(i+1,:)=t(i,:)-z1;
    v1(i+1,:)=v1(i,:)-z2;
      te=[zeros(1,n-m),v1(i+1,:)];
    h0(i+1,:)=t(i+1,:)+an*te;
end
%

%
% ------------------------------------------------------------------
%

%
function substitution_cr;
%
%                implements forward and backward phases
%
global m n lnn 
global col 
global rhs x
global fv0 fr ft fv1
global frti fv0ti dv1ti fhti
%
%                constants
%
nn=2^lnn;
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
an=col(n+1);
%
%       forward  phase
%
b=reshape(rhs,n,nn)';
w=zeros(nn/2,n);
%
parinv={frti(1,:),fv0ti(1,:),fhti(1,:),dv1ti(1,:)};
for j=1:nn/2
    fy=ifft(b(2*j,:),n2)*sn2;
    w(j,:)=moltinv(parinv,fy);
end
for i=1:lnn
     parinv={frti(i+1,:),fv0ti(i+1,:),fhti(i+1,:),dv1ti(i+1,:)};
     parr={fr(i,:),fv1(i,:)};
     parh1={ft(i,:),fv0(i,:)};
     for j=1:2^(lnn-i)
         z1=w(j,1:m);
         fz1=ifft(z1,m2)*sm2;
         z=moltr(parr,z1,fz1);
         j1=2^i*(j-1)+1;
         b(j1,n-m+1:n)=b(j1,n-m+1:n)-z;
         if j>=2 
             z1=w(j-1,:);
             fz1=ifft(z1,n2)*sn2;
             z=molth1(an,parh1,z1,fz1);
             b(j1,:)=b(j1,:)-z;
         end
     end
     for j=1:2^(lnn-i-1)
         j1=1+2^i*(2*j-1);
         fy=ifft(b(j1,:),n2)*sn2;
         w(j,:)=moltinv(parinv,fy);
     end
end
%
%      backward phase
%
x=zeros(nn,n);
y=b(1,:);
z=inverhcap(y);
x(1,:)=z(:,1)';
%
for i=lnn:-1:1
     parinv={frti(i,:),fv0ti(i,:),fhti(i,:),dv1ti(i,:)};
     parr={fr(i,:),fv1(i,:)};
     parh1={ft(i,:),fv0(i,:)};
     for j=2^(lnn-i):-1:1
         z1=x(j,:);
         fz1=ifft(z1,n2)*sn2;
         z=molth1(an,parh1,z1,fz1);
         j1=1+2^(i-1)*(2*j-1);
         vet=b(j1,:)-z;
         if j<2^(lnn-i)
             y=x(j+1,1:m);
             fy=ifft(y,m2)*sm2;
             vet(n-m+1:n)=vet(n-m+1:n)-moltr(parr,y,fy);
         end
         fvet=ifft(vet,n2)*sn2;
         x(2*j,:)=moltinv(parinv,fvet);
         x(2*j-1,:)=x(j,:);
     end
end
%

%
% ------------------------------------------------------------------
%

%
function z=inverh0(i)
%
%       computes the d-vectors of the H_0^-1 at level i
%
global m n
global h0 v0 r v1
%
rlu=[zeros(1,n-m),r(i,:)];
v1lu=[zeros(1,n-m),v1(i,:)];
%
cv(:,1)=-rlu';
cv(:,2)=h0(i,:)';
cv(:,3)=v1lu';
cv(:,4)=v0(i,:)';
%
dv(:,1)=[0,fliplr(v0(i,2:n))]';
dv(:,2)=[1,-fliplr(v1lu(2:n))]';
dv(:,3)=[0,fliplr(h0(i,2:n))]';
dv(:,4)=[0,fliplr(rlu(2:n))]';
%
vett(:,1)=rlu';
vett(:,2)=v0(i,:)';
vett(:,3)=v1lu';
vett(:,4)=eye(1,n)';
%
alpha=4;
z=levinson(n,alpha,cv,dv,4,vett);
%

%
% ------------------------------------------------------------------
%

%
function z=inverhcap(y)
%
%       solves the system cap H z= y at final level
%
global m n lnn 
global v0 hcap t r v1
%
lnn1=lnn+1;
%
rlu=[zeros(1,n-m),r(lnn1,:)];
rlu1=[zeros(1,n-m),r(1,:)];
v1lu=[zeros(1,n-m),v1(lnn1,:)];
%
cv(:,1)=hcap(lnn1,:)';
cv(:,2)=eye(1,n)';
cv(:,3)=-rlu';
cv(:,4)=v1lu';
%
dv(:,1)=eye(1,n)';
dv(:,2)=[0,fliplr(rlu1(2:n))]';
dv(:,3)=[0,fliplr(v0(lnn1,2:n))]';
dv(:,4)=[0,fliplr(t(lnn1,2:n))]';
%
vett(:,1)=y';
%
alpha=4;
z=levinson(n,alpha,cv,dv,1,vett);
%

%
% ------------------------------------------------------------------
%

%
function fou=molth1(an,par,vetnt,vet)
%      computes the product of H1 by vector vetnt of size n
%      whose transformed vector of size 2n is vet 
%      par contains the transformed d-vectors of H1
global n
%
n2=2*n;
sn2=sqrt(n2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(n+1:n2),n2)*sn2;
fp1=par{1}.*te;
%
fp=par{1}.*vet;
te=fft(fp);
te=ifft(an*vetnt+te(n+1:n2),n2)*sn2;
fp2=par{2}.*te;
%
te=fft(fp2-fp1);
fou=te(1:n);
%

%
% ------------------------------------------------------------------
%

%
function fou=moltr(par,vetnt,vet)
%      computes the product of R by vector vetnt of size m
%      whose transformed vector of size 2m is vet 
%      par contains the transformed d-vectors of R
global m
%
m2=2*m;
sm2=sqrt(m2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(-vetnt+te(m+1:m2),m2)*sm2;
fp1=par{1}.*te;
%
fp=par{1}.*vet;
te=fft(fp);
te=ifft(te(m+1:m2),m2)*sm2;
fp2=par{2}.*te;
%
te=fft(fp2-fp1);
fou=te(1:m);
%

%
% ------------------------------------------------------------------
%

%
function fou=moltinv(par,vet)
%      computes the product of H0^-1 by a vector of size n
%      whose transformed vector of size 2n is vet 
%      par contains the transformed d-vectors of H0^-1
global n 
%
n2=2*n;
sn2=sqrt(n2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(n+1:n2),n2)*sn2;
fp1=par{1}.*te;
%
fp=par{1}.*vet;
te=fft(fp);
te=ifft(te(n+1:n2),n2)*sn2;
fp2=par{2}.*te;
%
fp3=par{3}.*vet;
%
fp4=par{4}.*vet;
%
te=fft(fp1-fp2+fp3+fp4);
fou=te(1:n);
%

%
% ------------------------------------------------------------------
%

%
function z=levinson(nz,alpha,cv,dv,hz,vett)
%  Levinson algorithm for solving a Toeplitz-like system
%  with hz right hand-sides of size nz. 
%  vett contains columnwise the right hand-sides.
%  z contains the solution.
%  the d-vectors of the Toeplitz-like matrix are contained
%  columnwise in the nz times alpha matrices cv and dv

ee=dot(cv(1,1:alpha)',dv(1,1:alpha)'); 
%
%  ee is the Schur complement and should be nonzero
%
z(1,1:hz)=vett(1,1:hz)/ee;              
y=-dot(cv(1,1:alpha)',dv(2,1:alpha)')/ee;
s=dot(cv(2,1:alpha)',dv(1,1:alpha)');
f=cv(1,1:alpha)/ee;
v=-dv(2,1:alpha);
%
for k=2:nz-1
    t=cv(k,1:alpha)-s*f;
    ee=ee-dot(t',v);
    w=(vett(k,1:hz)-s*z)/ee;
    yt=y;
    yt(k)=1;
    z(k,1:hz)=zeros(1,hz);
    z=z+yt'*w;
    s(2:k)=s;
    s(1)=0;
    s=s+cv(k+1,1:alpha)*dv(1:k,1:alpha)';
    f(k,1:alpha)=zeros(1,alpha);
    f=f+yt'*t/ee;
    v=-yt*dv(2:k+1,1:alpha);
    y(2:k)=y;
    y(1)=0;
    y=y+v*f';
end
t=cv(nz,1:alpha)-s*f;
ee=ee-dot(t',v);
w=(vett(nz,1:hz)-s*z)/ee;
yt=y;
yt(nz)=1;
z(nz,1:hz)=zeros(1,hz);
z=z+yt'*w;