function program_mcr; 
%  program mcr implements algorithm MCR 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 rti v0ti v1ti hti
global fv0 fr ft fv1
global frti fv0ti fv1ti fhti
global dv0 dt dr dv1 dv1ti
%
%  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_mcr; 
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_mcr;
%
%    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 rti v0ti v1ti hti
global fv0 fr ft fv1
global frti fv0ti fv1ti fhti
global dv0 dt dr dv1 dv1ti
%
%       constants
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
an=col(n+1);
%
%        utility vectors
%
h0=zeros(lnn+1,n); v0=zeros(lnn+1,n); hcap=zeros(lnn+1,n);
t=zeros(lnn+1,n); r=zeros(lnn+1,m); v1=zeros(lnn+1,m);
%
v0ti=zeros(lnn+1,n); hti=zeros(lnn+1,n); rti=zeros(lnn+1,n);
v1ti=zeros(lnn+1,n);
%
fv0=zeros(lnn+1,n2); ft=zeros(lnn+1,n2); fr=zeros(lnn+1,m2);
fv1=zeros(lnn+1,m2);
%
fv0ti=zeros(lnn+1,n2); fhti=zeros(lnn+1,n2); frti=zeros(lnn+1,n2);
fv1ti=zeros(lnn+1,n2);
%
dv0=zeros(lnn+1,n2); dt=zeros(lnn+1,n2); dr=zeros(lnn+1,m2);
dv1=zeros(lnn+1,m2); dv1ti=zeros(lnn+1,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);
%
%                           FFT of the d-vectors at first level
%
fv0(1,:)=ifft(v0(1,:),n2)*sn2;    
ft(1,:)=ifft(t(1,:),n2)*sn2;
fr(1,:)=ifft(r(1,:),m2)*sm2;
fv1(1,:)=ifft(v1(1,:),m2)*sm2;
%
 te=[zeros(1,n+1),v0(1,2:n)];
dv0(1,:)=ifft(te)*sn2;
 te=[zeros(1,n+1),t(1,2:n)];
dt(1,:)=ifft(te)*sn2;
 te=[zeros(1,m+1),r(1,2:m)];
dr(1,:)=ifft(te)*sm2;
 te=[zeros(1,m+1),v1(1,2:m)];
dv1(1,:)=ifft(te)*sm2;
%
%                    d-vectors at current level are already  known                              
for i=1:lnn
    %
    %    computes the d-vectors of H0^-1 at current level
    %
    parh1={ft(i,:),fv0(i,:)};
    parr={fr(i,:),fv1(i,:)};
    % 
    if i==1
        z=inverh0(i);
    else
        z=inverh0mcr(i);
    end 
    rti(i,:)=z(:,1)';
    v0ti(i,:)=z(:,2)';
    v1ti(i,:)=z(:,3)';
    hti(i,:)=z(:,4)';
    %
    %                    computes the corresponding FFTs
    %
    frti(i,:)=ifft(rti(i,:),n2)*sn2;
    fv0ti(i,:)=ifft(v0ti(i,:),n2)*sn2;
    fv1ti(i,:)=ifft(v1ti(i,:),n2)*sn2;
    fhti(i,:)=ifft(hti(i,:),n2)*sn2;
       te=[zeros(1,n+1),v1ti(i,2:n)];
    dv1ti(i,:)=ifft(te)*sn2;
    %
    %           computes the d-vectors for the next level
    %                   
       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;
    %
    %                    computes the corresponding FFTs
    %
    fv0(i+1,:)=ifft(v0(i+1,:),n2)*sn2;    
    ft(i+1,:)=ifft(t(i+1,:),n2)*sn2;
    fr(i+1,:)=ifft(r(i+1,:),m2)*sm2;
    fv1(i+1,:)=ifft(v1(i+1,:),m2)*sm2;
       te=[zeros(1,n+1),v0(i+1,2:n)];
    dv0(i+1,:)=ifft(te)*sn2;
       te=[zeros(1,n+1),t(i+1,2:n)];
    dt(i+1,:)=ifft(te)*sn2;
       te=[zeros(1,m+1),r(i+1,2:m)];
    dr(i+1,:)=ifft(te)*sm2;
       te=[zeros(1,m+1),v1(i+1,2:m)];
    dv1(i+1,:)=ifft(te)*sm2;
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)
%
% compute the d-vectors of the inverse of matrix H_0 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 tou=inverh0mcr(i)
%
% compute the d-vectors of the inverse of matrix H_0 at  level i
% by using SMW formula
%
global m n
global col
global h0 v0 hcap t r v1
global rti v0ti v1ti hti
global fv0 fr ft fv1
global frti fv0ti fv1ti fhti
global dv0 dt dr dv1 dv1ti
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
fon=ones(1,n2)/sn2;
on1=zeros(1,n);
on1(n-m+1)=1;
fon1=ifft(on1,n2)*sn2;
an=col(n+1);
%
%      set parameters for calling multiplication routines
%
parinv={frti(i-1,:),fv0ti(i-1,:),fhti(i-1,:),dv1ti(i-1,:)};
parv={fv1(i,:),dt(i,:),fr(i,:),dv0(i,:),fr(i-1,:),dv0(i-1,:), ...
        fv1(i-1,:),dt(i-1,:)};
paru={ft(i,:),dv1(i,:),fv0(i,:),dr(i,:),fv0(i-1,:),dr(i-1,:), ...
        ft(i-1,:),dv1(i-1,:)};
%
%             computes right hand-sides vectors
%
ou=zeros(4,n);
vett=zeros(m2,4);
%
ou(4,:)=moltinv(parinv,fon);
fy=ifft(ou(4,:),n2)*sn2;
vett(1:m,4)=moltv(an,parv,ou(4,:),fy)';
vett(m+1:m2,4)=ou(4,1:m)';
%
te=[zeros(1,n-m),r(i,:)];
fy=ifft(te,n2)*sn2;
ou(1,:)=moltinv(parinv,fy);
fy=ifft(ou(1,:),n2)*sn2;
vett(1:m,1)=moltv(an,parv,ou(1,:),fy)';
vett(m+1:m2,1)=ou(1,1:m)';
%
te=[zeros(1,n-m),v1(i,:)];
fy=ifft(te,n2)*sn2;
ou(3,:)=moltinv(parinv,fy);
fy=ifft(ou(3,:),n2)*sn2;
vett(1:m,3)=moltv(an,parv,ou(3,:),fy)';
vett(m+1:m2,3)=ou(3,1:m)';
%
ou(2,:)=moltinv(parinv,fv0(i,:));
fy=ifft(ou(2,:),n2)*sn2;
vett(1:m,2)=moltv(an,parv,ou(2,:),fy)';
vett(m+1:m2,2)=ou(2,1:m)';
%            
%            computes  d-vectors of P(i)
%
cv=zeros(m2,9);
dv=zeros(m2,9);
%
y=moltinv(parinv,fon1);
fy=ifft(y,n2)*sn2;
cv(1:m,1)=-moltv(an,parv,y,fy)';
cv(m+1:m2,1)=-y(1:m)';
cv(1,1)=cv(1,1)+1;
%
y=moltinv(parinv,fv0(i,:));
fy=ifft(y,n2)*sn2;
cv(1:m,2)=-moltv(an,parv,y,fy)';
cv(m+1:m2,2)=-y(1:m)';
%
cv(1:m,3)=-r(i,:)';
%
cv(1:m,4)=-moltv(an,parv,v0ti(i-1,:),fv0ti(i-1,:))';
cv(m+1:m2,4)=-v0ti(i-1,1:m)';
%
cv(1:m,5)=-moltv(an,parv,rti(i-1,:),frti(i-1,:))';
cv(m+1:m2,5)=-rti(i-1,1:m)';
cv(1:m,5)=cv(1:m,5)+r(i-1,1:m)';
%
fz1=-ft(i,:)+ft(i-1,:);
y=moltinv(parinv,fz1);
fy=ifft(y,n2)*sn2;
cv(1:m,6)=-moltv(an,parv,y,fy)';
cv(m+1:m2,6)=-y(1:m)';
if i==2, cv(m+1,6)=cv(m+1,6)+1; end
%
cv(1:m,7)=-v1(i,:)';
cv(m+1,7)=1;
%
y=moltv(an,parv,v1ti(i-1,:),fv1ti(i-1,:));
cv(1:m,8)=v1(i,1:m)'-v1(i-1,1:m)'+y';
cv(m+1:m2,8)=v1ti(i-1,1:m)';
%
cv(m+1,9)=1;
%
dv(1,1)=1;
dv(2:m2,2)=flipud(cv(2:m2,3));
dv(2:m2,3)=-flipud(cv(2:m2,2));
dv(2:m2,4)=flipud(cv(2:m2,5));
dv(2:m2,5)=-flipud(cv(2:m2,4));
dv(2:m2,6)=flipud(cv(2:m2,7));
dv(2:m2,7)=-flipud(cv(2:m2,6));
dv(2:m2,8)=flipud(cv(2:m2,9));
dv(2:m2,9)=-flipud(cv(2:m2,8));
%
if i==2, rango=7; else rango=9; end
%
z=levinson(m2,rango,cv,dv,4,vett);
%
for j=1:4 
    y=z(m+1:m2,j)';
    fy=ifft(y,m2)*sm2;
    yy1=moltu(paru,y,fy);
    yy1(n-m+1:n)=yy1(n-m+1:n)+z(1:m,j)';
    fy=ifft(yy1,n2)*sn2;
    ou(j,:)=ou(j,:)+moltinv(parinv,fy);
end
tou=ou';
%

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

%
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=moltu(par,vetnt,vet)
%      computes the product of U1 by vector vetnt of size m
%      whose transformed vector of size 2m is vet 
%      par contains the transformed d-vectors of U1
global m n
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(1:m)-vetnt,n2)*sn2;
fp1=par{1}.*te;
%
fp=par{4}.*vet;
te=fft(fp);
te=ifft(te(1:m),n2)*sn2;
fp2=par{3}.*te;
%
fp=par{6}.*vet;
te=fft(fp);
te=ifft(te(1:m),n2)*sn2;
fp3=par{5}.*te;
%
fp=par{8}.*vet;
te=fft(fp);
te=ifft(te(1:m)-vetnt,n2)*sn2;
fp4=par{7}.*te;
%
te=fft(fp1-fp2+fp3-fp4);
fou=te(1:n);
%

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

%
function fou=moltv(an,par,vetnt,vet)
%      computes the product of V1 by vector vetnt of size n
%      whose transformed vector of size 2n is vet 
%      par contains the transformed d-vectors of V1
global m n
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(1:m)+an*vetnt(1:m),m2)*sm2;
fp1=par{1}.*te;
%
fp=par{4}.*vet;
te=fft(fp);
te=ifft(te(1:m),m2)*sm2;
fp2=par{3}.*te;
%
fp=par{6}.*vet;
te=fft(fp);
te=ifft(te(1:m),m2)*sm2;
fp3=par{5}.*te;
%
fp=par{8}.*vet;
te=fft(fp);
te=ifft(te(1:m)+an*vetnt(1:m),m2)*sm2;
fp4=par{7}.*te;
%
te=fft(-fp1+fp2-fp3+fp4);
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;