function program_st;
%  program st implements algorithm ST 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 fb1cap fb2cap
global xup vup vdo lup 
global gdo pup pdo qup qdo 
global fxup fxdo flup fldo
global fp1ti fp2ti fp3ti fp4ti
global fme1 fme3
global dudo dvup dvdo dgi dhi dwup
global dkup dkdo
global dq1ti dq2ti dq3ti dme2
%
%  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_st; 
substitution_st; 
%
%  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_st;
%    creates the d-vectors and the coresponding FFTs
%    at all levels
global m n lnn 
global row col 
global fb1cap fb2cap
global xup vup vdo lup 
global gdo pup pdo qup qdo 
global fxup fxdo flup fldo
global fp1ti fp2ti fp3ti fp4ti
global fme1 fme3
global dudo dvup dvdo dgi dhi dwup
global dkup dkdo
global dq1ti dq2ti dq3ti dme2
%
%       constants
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
lnn1=lnn+1;
an=col(n+1);
  te=[zeros(1,n+1),col(2:n)];
dvb0=ifft(te)*sn2;
%
%    utility vectors
%
xdo=zeros(lnn1,n); uup=zeros(lnn1,n); udo=zeros(lnn1,n);
ldo=zeros(lnn1,n); gi=zeros(lnn1,n); wup=zeros(lnn1,n);
kup=zeros(lnn1,n);
%
xup=zeros(lnn1,m); vup=zeros(lnn1,m); vdo=zeros(lnn1,m);
lup=zeros(lnn1,m); hi=zeros(lnn1,m); gdo=zeros(lnn1,m);
pup=zeros(lnn1,m); pdo=zeros(lnn1,m); qup=zeros(lnn1,m);
qdo=zeros(lnn1,m); kdo=zeros(lnn1,m);
%
me1=zeros(lnn1,n); me3=zeros(lnn1,n);
%
p1ti=zeros(lnn1,m); p2ti=zeros(lnn1,m); p3ti=zeros(lnn1,m);
p4ti=zeros(lnn1,m); q1ti=zeros(lnn1,m); q2ti=zeros(lnn1,m);
q3ti=zeros(lnn1,m); me2=zeros(lnn1,m);
%
fxdo=zeros(lnn1,n2); fldo=zeros(lnn1,n2);
%
fxup=zeros(lnn1,m2); flup=zeros(lnn1,m2);
%
fme1=zeros(lnn1,n2); fme3=zeros(lnn1,n2);
%
fp1ti=zeros(lnn1,m2); fp2ti=zeros(lnn1,m2); fp3ti=zeros(lnn1,m2);
fp4ti=zeros(lnn1,m2); 
%
dudo=zeros(lnn1,n2); dgi=zeros(lnn1,n2); dwup=zeros(lnn1,n2); 
dkup=zeros(lnn1,n2);
%
dvup=zeros(lnn1,m2); dvdo=zeros(lnn1,m2); dhi=zeros(lnn1,m2);
dkdo=zeros(lnn1,m2);
%
dq1ti=zeros(lnn1,m2); dq2ti=zeros(lnn1,m2); dq3ti=zeros(lnn1,m2);
dme2=zeros(lnn1,m2);
%
%    computes the d-vectors at first level 
%
z=inverb0;
%
b1cap=z(:,1)';
b2cap=z(:,2)';
%
fb1cap=ifft(b1cap,n2)*sn2;
fb2cap=ifft(b2cap,n2)*sn2;
%
%    computes the d-vectors of the corner blocks at first level
%
xup(1,:)=b1cap(1:m);
fxup(1,:)=ifft(xup(1,:),m2)*sm2;
xdo(1,:)=b1cap;
fxdo(1,:)=fb1cap;
%
rb1=fliplr(col(2:n+1));
frb1=ifft(rb1,n2)*sn2;
  te=ifft(fliplr(b1cap),n2)*sn2;
  te=fft(frb1.*te,n2);
uup(1,:)=te(1:n);
  te=ifft(fliplr(b2cap),n2)*sn2;
  te=fft(frb1.*te);
udo(1,:)=te(1:n);
  te=[zeros(1,n+1),fliplr(udo(1,1:n-1))];
dudo(1,:)=ifft(te)*sn2;
  te=[row(m+1),zeros(1,m),row(2:m)];
dvrr1=ifft(te,m2)*sm2;
  te=fliplr(b1cap);
  te=ifft(te(n-m+1:n),m2)*sm2;
  te=fft(dvrr1.*te);
vup(1,:)=te(1:m);
  te=[zeros(1,m+1),fliplr(vup(1,1:m-1))];
dvup(1,:)=ifft(te)*sm2;
  te=fliplr(b2cap);
  te=ifft(te(n-m+1:n),m2)*sm2;
  te=fft(dvrr1.*te);
vdo(1,:)=te(1:m);
  te=[1,zeros(1,m),fliplr(vdo(1,1:m-1))];
dvdo(1,:)=ifft(te)*sm2;
wup(1,:)=fliplr(col(1:n))-uup(1,1:n);
  te=[zeros(1,n+1),fliplr(uup(1,1:n-1))];
  te=ifft(te)*sn2;
dwup(1,:)=dvb0-te+an/sn2;
ldo(1,:)=b2cap;
fldo(1,:)=fb2cap;
lup(1,:)=b2cap(1:m);
flup(1,:)=ifft(b2cap(1:m),m2)*sm2;
gi(1,:)=-fliplr(eye(1,n));
  te=[zeros(1,n+1),fliplr(gi(1,1:n-1))];
dgi(1,1:n2)=ifft(te)*sn2;
hi(1,:)=zeros(1,m);
dhi(1,:)=zeros(1,m2);
% 
for i=1:lnn
    %
    %                    computes the d-vectors at current level
    %
    parhdo={fxdo(i,:),dvdo(i,:),fldo(i,:),dvup(i,:),dhi(i,:)};
    pargup={fxup(i,:),dudo(i,:),flup(i,:),dwup(i,:)};
    parghtra={fxup(i,:),dudo(i,:),dvdo(i,:),flup(i,:),dwup(i,:), ...
            dvup(i,:)};
    pargdo={fxdo(i,:),dudo(i,:),fldo(i,:),dwup(i,:),dgi(i,:)};
    parhh={fxup(i,:),fxdo(i,:),dvdo(i,:),flup(i,:),fldo(i,:), ...
            dvup(i,:),dhi(i,:)};
    %
    %                    computes the d-vectors of P(i)    
    %
      vet=fldo(i,:);
    pup(i,:)=moltgup(pargup,vet);
      vet=fxdo(i,:);
    pdo(i,:)=eye(1,m)-moltgup(pargup,vet);
      vet=ifft(uup(i,:),n2)*sn2;
      te=uup(i,1:m);
      vet1=ifft(te,m2)*sm2;
    qup(i,:)=row(2:m+1)-molthdotra(parhdo,vet,vet1);
    %
    
      vet=ifft(udo(i,:),n2)*sn2;
      te=udo(i,1:m);
      vet1=ifft(te,m2)*sm2;
    qdo(i,:)=molthdotra(parhdo,vet,vet1);
    %
    %                  returns the d-vectors of P(i)^-1
    %
    z=inverp(i);    
    %                  
    p1ti(i,1:m)=z(:,1)';
    p2ti(i,1:m)=z(:,2)';
    p3ti(i,1:m)=z(:,3)';
    p4ti(i,1:m)=z(:,4)';
    %
      vet=ifft(gi(i,1:n),n2)*sn2;
      te=gi(i,1:m);
      vet1=ifft(te,m2)*sm2;
    gdo(i,:)=molthdotra(parhdo,vet,vet1);
    %        
    z=inverptra(i);  
    %                  
    q1ti(i,1:m)=z(:,1)';
    q2ti(i,1:m)=z(:,2)';
    q3ti(i,1:m)=z(:,3)';
    me2(i,1:m)=z(:,4)';
    %
    %    computes the required FFTs of the d-vectors at current level
    %
    fp1ti(i,:)=ifft(p1ti(i,:),m2)*sm2;   
    fp2ti(i,:)=ifft(p2ti(i,:),m2)*sm2;   
    fp3ti(i,:)=ifft(p3ti(i,:),m2)*sm2; 
    fp4ti(i,:)=ifft(p4ti(i,:),m2)*sm2; 
      te=[zeros(1,m+1),fliplr(q1ti(i,1:m-1))];
    dq1ti(i,:)=ifft(te)*sm2;
      te=[zeros(1,m+1),fliplr(q2ti(i,1:m-1))];
    dq2ti(i,:)=ifft(te)*sm2;
      te=[zeros(1,m+1),fliplr(q3ti(i,1:m-1))];
    dq3ti(i,:)=ifft(te)*sm2;
      te=[zeros(1,m+1),fliplr(me2(i,1:m-1))];
    dme2(i,:)=ifft(te)*sm2;
      te=ifft(q3ti(i,:),m2)*sm2; 
    [foua,foub]=moltghtra(parghtra,te);
    kdo(i,:)=foub;
      te=ifft(foua,n2)*sn2;
    kup(i,:)=moltgdotra(pargdo,te);
      te=[zeros(1,m+1),fliplr(kdo(i,1:m-1))];
    dkdo(i,:)=ifft(te)*sm2;
      te=[zeros(1,n+1),fliplr(kup(i,1:n-1))];
    dkup(i,:)=ifft(te)*sn2;
      te=molthdo(parhdo,fp4ti(i,:));
      te=ifft(te,n2)*sn2;
    me3(i,:)=moltgdo(pargdo,te);
    fme3(i,:)=ifft(me3(i,:),n2)*sn2;
      te=molthdo(parhdo,fp2ti(i,:));
      te=ifft(te,n2)*sn2;
      te=moltgdo(pargdo,te);
    me1(i,:)=ldo(i,:)+te;
    fme1(i,:)=ifft(me1(i,:),n2)*sn2;
    %
    %     computes the d-vectors of the corner blocks for the
    %     next level and corresponding FFTs
    %
    [foua,foub]=molthh(parhh,fp1ti(i,:));
    xup(i+1,:)=-foua;
    fxup(i+1,:)=ifft(xup(i+1,:),m2)*sm2;
      te=ifft(foub,n2)*sn2;
      te=moltgdo(pargdo,te);
    xdo(i+1,:)=xdo(i,:)+te;
    fxdo(i+1,:)=ifft(xdo(i+1,:),n2)*sn2;
    %
    [foua,foub]=molthh(parhh,fp3ti(i,:));
    lup(i+1,:)=lup(i,:)+foua;
    flup(i+1,:)=ifft(lup(i+1,:),m2)*sm2;
      te=ifft(foub+ldo(i,:),n2)*sn2;
    ldo(i+1,:)=-moltgdo(pargdo,te);
    fldo(i+1,:)=ifft(ldo(i+1,:),n2)*sn2;
    %
      te=ifft(q1ti(i,:),m2)*sm2;   
    [foua,foub]=moltghtra(parghtra,te);
      te=ifft(foua+udo(i,:),n2)*sn2;
    udo(i+1,:)=-moltgdotra(pargdo,te);
      te=[zeros(1,n+1),fliplr(udo(i+1,1:n-1))];
    dudo(i+1,:)=ifft(te)*sn2;
    vdo(i+1,:)=vdo(i,:)+foub;
    te=[1,zeros(1,m),fliplr(vdo(i+1,1:m-1))];
    dvdo(i+1,:)=ifft(te)*sm2;
    %
      te=ifft(q2ti(i,:),m2)*sm2; 
    [foua,foub]=moltghtra(parghtra,te);
      te=ifft(foua,n2)*sn2;
    uup(i+1,:)=uup(i,:)+moltgdotra(pargdo,te);
    vup(i+1,:)=-foub;
      te=[zeros(1,m+1),fliplr(vup(i+1,1:m-1))];
    dvup(i+1,:)=ifft(te)*sm2;
    %
    wup(i+1,1:n)=fliplr(col(1:n))-uup(i+1,1:n);
       te=[zeros(1,n+1),fliplr(uup(i+1,1:n-1))];
       te=ifft(te)*sn2;
    dwup(i+1,:)=dvb0-te+an/sn2;
      te=ifft(me2(i,:),m2)*sm2; 
    [foua,foub]=moltghtra(parghtra,te);
      te=ifft(foua+gi(i,:),n2)*sn2;
    gi(i+1,:)=-moltgdotra(pargdo,te);
      te=[zeros(1,n+1),fliplr(gi(i+1,1:n-1))];
    dgi(i+1,:)=ifft(te)*sn2;
    hi(i+1,:)=hi(i,:)+foub;
      te=[zeros(1,m+1),fliplr(hi(i+1,1:m-1))];
    dhi(i+1,:)=ifft(te)*sm2;
end
%

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

%
function substitution_st;
%
%    implements forward and backward phases
%
global m n lnn 
global row col 
global rhs x
global fb1cap fb2cap
global fxup fxdo flup fldo
global fp1ti fp2ti fp3ti fp4ti
global fme1 fme3
global dudo dvup dvdo dgi dhi dwup
global dkup dkdo
global dq1ti dq2ti dq3ti dme2
%
%                constants
%
nn=2^lnn; dim= n*nn;
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
an=col(n+1);
%
xsop=zeros(nn,m);
xsot=zeros(nn,n);
fxsop=zeros(nn,m2);
fxsot=zeros(nn,n2);
ysot=zeros(nn,n);
py=zeros(nn,m);
xn1=zeros(1,n);
  te=[col(n+1),zeros(1,n),col(2:n)];
drb1=ifft(te)*sn2;
  te=[zeros(1,m+1),fliplr(row(2:m))];
dvrr=ifft(te)*sm2;
%
%       forward  phase
%
 for i=1:nn
      fxn1=ifft(xn1,n2)*sn2;
      te=fft(drb1.*fxn1);
      tni=rhs((i-1)*n+1:i*n)-te(1:n);
      vet=ifft(tni,n2)*sn2;
      xn1=moltb0inv(fb1cap,fb2cap,tni,vet);
    ysot(i,:)=xn1;
    if mod(i,2)==0
        j=i;
        k=1;
        while mod(j,2)==0
            j=j/2;
            jk=j+nn*(1-2/2^k);
            jk2=2*j+nn*(1-4/2^k);
            if k==1, te=ysot(2*j,1:m); else te=xsop(jk2,:); end
            vet=ifft(te,m2)*sm2;
            vet1=ifft(te,n2)*sn2;
            parpim={fp4ti(k,:),fme3(k,:),fp1ti(k,:),fxdo(k+1,:), ...
                    dq1ti(k,:),fp2ti(k,:),fme1(k,:),dvrr,fp3ti(k,:), ...
                fldo(k+1,:),dq2ti(k,:),dq3ti(k,:),dme2(k,:)};
            [foua,foub]=moltpim(parpim,vet,vet1);
            py(jk,:)=foua;
            if k==1, te=ysot(2*j,:); else te=xsot(jk2,:); end
            xn1=te+foub;
            xsot(jk,:)=xn1;
            vet=ifft(py(jk,:),m2)*sm2;
            parhup={fxup(k,:),dvdo(k,:),flup(k,:),dvup(k,:)};
            fou=molthup(parhup,vet);
            if k==1, te=ysot(2*j-1,1:m); else te=xsop(jk2-1,:); end
            xsop(jk,:)=te-fou;
            k=k+1;
        end
    end
end
 %
 %      backward phase
 %
 for i=lnn:-1:2
     parlk={fp1ti(i,:),dudo(i+1,:),fp3ti(i,:),dwup(i+1,:), ...
             dudo(i,:)+dkup(i,:),dvdo(i+1,:),dvup(i+1,:),dkdo(i,:)};
     parghdo={fxdo(i,:),dudo(i,:),fldo(i,:),dwup(i,:),dgi(i,:), ...
             dvdo(i,:),dvup(i,:),dhi(i,:)};
     for j=1:2^(lnn-i)
         ji=j+nn*(1-2/2^i);
         ji2=2*j+nn*(1-4/2^i);
         xsop(ji2-1,:)=xsop(ji,:);
         fxsop(ji2-1,:)=fxsop(ji,:);
         xsot(ji2,:)=xsot(ji,:);
         if j==1
             vet1=zeros(1,n2);
         else
             k=2*j-2;
             h=0;
             while mod(k,2)==0
                 k=k/2;
                 h=h+1;
             end
             vet1=fxsot(k+nn*(1-4/2^(i+h)),:);
         end
         vet2=fxsop(ji+1,:);
         fou=moltlk(parlk,vet1,vet2);
         xsop(ji2,:)=py(ji,:)+fou;
         fxsop(ji2,:)=ifft(xsop(ji2,:),m2)*sm2;
         vet2=fxsop(ji2,:);
         vet=-moltghdo(parghdo,vet1,vet2);
         xsot(ji2-1,:)=xsot(ji2-1,:)+vet;
         fxsot(ji2-1,:)=ifft(vet,n2)*sn2;
     end
 end
     for j=nn/2:-1:1 
         xsot(2*j,:)=xsot(j,:);
        if j==1
             vet1=zeros(1,n2);
         else
             k=2*j-2;
             h=0;
             while mod(k,2)==0
                 k=k/2;
                 h=h+1;
             end
             vet1=fxsot(k+nn*(1-2/2^h),:);
         end
         vet2=ifft(xsot(j,1:m),m2)*sm2;
         parghdo1={fxdo(1,:),dudo(1,:),fldo(1,:),dwup(1,:),dgi(1,:), ...
                 dvdo(1,:),dvup(1,:),dhi(1,:)};
         vet=-moltghdo(parghdo1,vet1,vet2);
         xsot(2*j-1,:)=ysot(2*j-1,:)+vet;
     end
      x=xsot;   
%

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

%
function z=inverb0
%
%       computes the d-vectors of B_0^-1
%
global m n 
global row col 
%
cv(:,1)=col(1:n)';
cv(:,2)=eye(1,n)';
%
dv(:,1)=eye(1,n)';
dv(:,2)=[0,row(2:m+1),zeros(1,n-m-1)]';
%
vett(:,1)=[zeros(1,n-m),fliplr(row(2:m+1))]';
vett(:,2)=eye(1,n)';
%
z=levinson(n,2,cv,dv,2,vett);
%

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

%
function z=inverp(i)
%
%         computes left d-vectors of P(i)^-1
%
global m
global xup vup vdo lup 
global pup pdo qup qdo 
%
alpha=4;
cv=zeros(m,alpha);
dv=zeros(m,alpha);
%
cv(:,1)=-xup(i,1:m)';
cv(:,2)=-lup(i,1:m)';
cv(:,3)=pup(i,1:m)';
cv(:,4)=pdo(i,1:m)';
%
dv(2:m,1)=qdo(i,1:m-1)';
dv(2:m,2)=qup(i,1:m-1)';
dv(2:m,3)=vup(i,1:m-1)';
dv(2:m,4)=vdo(i,1:m-1)';
dv(1,4)=1;
%
vett(:,1)=xup(i,1:m)';
vett(:,2)=lup(i,1:m)';
vett(:,3)=pup(i,1:m)';
vett(:,4)=eye(1,m)';
%
z=levinson(m,alpha,cv,dv,4,vett);
%

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

%
function z=inverptra(i)
%
%               computes right d-vectors of P(i)^-1
%
global m
global xup vup vdo lup 
global gdo pup pdo qup qdo 
%
alpha=4;
cv=zeros(m,alpha);
dv=zeros(m,alpha);
%
dv(:,1)=-xup(i,1:m)';
dv(:,2)=-lup(i,1:m)';
dv(:,3)=pup(i,1:m)';
dv(:,4)=pdo(i,1:m)';
%
cv(2:m,1)=qdo(i,1:m-1)';
cv(2:m,2)=qup(i,1:m-1)';
cv(2:m,3)=vup(i,1:m-1)';
cv(2:m,4)=vdo(i,1:m-1)';
cv(1,4)=1;
%
vett(:,1)=qdo(i,1:m)';
vett(:,2)=vup(i,1:m)';
vett(:,3)=vdo(i,1:m)';
vett(:,4)=gdo(i,1:m)';
%
z=levinson(m,alpha,cv,dv,4,vett);
%

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

%
function fou=moltb0inv(u,ut,vetnt,vet)
%      computes the product of B0^-1 by vector vetnt of size n
%      whose transformed vector of size 2n is vet 
%      u and ut are the transformed d-vectors of B0^-1
global n
%
n2=2*n;
sn2=sqrt(n2);
%
fp=ut.*vet;
te=fft(fp);
te=ifft(te(n+1:n2),n2)*sn2;
fp1=u.*te;
%
fp=u.*vet;
te=fft(fp);
te=ifft(-vetnt+te(n+1:n2),n2)*sn2;
fp2=ut.*te;
%
te=fft(fp1-fp2);
fou=te(1:n);

%

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

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

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

%
function fou=moltgdotra(par,vet)
%         computes the product of underlined G^T by a vector of size n
%         whose transformed vector of size 2n is vet
%         par contains the transformed d-vectors of underlined G

global n
%
n2=2*n;
sn2=sqrt(n2);
%
fp=conj(par{1}).*vet;
te=fft(fp);
te=ifft(te(1:n),n2)*sn2;
fp1=conj(par{2}).*te;
%
fp=conj(par{3}).*vet;
te=fft(fp);
te=ifft(te(1:n),n2)*sn2;
fp2=conj(par{4}).*te;
%
fp3=conj(par{5}).*vet;
%
te=fft(fp1+fp2+fp3);
fou=te(1:n);
%

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

%
function fou=moltghdo(par,vet,vet1)
%         computes the product of underlined G by a vector of size n 
%         whose transformed vector of size 2n is vet 
%         plus the product of underlined H by a vector of size m
%         whose transformed vector of size 2m is vet1 
%         par contains the transformed d-vectors of underlined G

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:n),n2)*sn2;
fp1=par{1}.*te;
%
fp=par{4}.*vet;
te=fft(fp);
te=ifft(te(1:n),n2)*sn2;
fp2=par{3}.*te;
%
fp3=par{5}.*vet;
%
%
fp=par{6}.*vet1;
te=fft(fp);
te=ifft(te(1:m),n2)*sn2;
fp4=par{1}.*te;
%
fp=par{7}.*vet1;
te=fft(fp);
te=ifft(te(1:m),n2)*sn2;
fp5=par{3}.*te;
%
fp6=par{8}.*vet1;
%
te=fft(fp1+fp2+fp3+fp4-fp5);
fou=te(1:n);
te=fft(fp6);
fou(1:m)=fou(1:m)+te(1:m);

%

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

%
function [foua,foub]=moltghtra(par,vet)
%      computes the product of overlined G^T by a vector of size m 
%      whose transformed vector of size 2m is vet 
%      and the product of overlined H^T by the same vector
%      par contains the transformed d-vectors of overlined G
%      and overlined H
global m n
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
%
fp=conj(par{1}).*vet;
te=fft(fp);
te1=ifft(te(1:m),n2)*sn2;
fp1a=conj(par{2}).*te1;
te1=ifft(te(1:m),m2)*sm2;
fp1b=conj(par{3}).*te1;
%
fp=conj(par{4}).*vet;
te=fft(fp);
te1=ifft(te(1:m),n2)*sn2;
fp2a=conj(par{5}).*te1;
te1=fft(fp1a+fp2a);
foua=te1(1:n);
te1=ifft(te(1:m),m2)*sm2;
fp2b=conj(par{6}).*te1;
%
te1=fft(fp1b-fp2b);
foub=te1(1:m);
%

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

%
function fou=moltgup(par,vet)
%         computes the product of overlined G by a vector of size n 
%         whose transformed vector of size 2n is vet
%         par contains the transformed d-vectors of overlined G
global m
%
m2=2*m;
sm2=sqrt(m2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(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;
%
te=fft(fp1+fp2);
fou=te(1:m);
%

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

%
function fou=molthdo(par,vet)
%         computes the product of underlined H by a vector of size m
%         whose transformed vector of size 2m is vet
%         par contains the transformed d-vectors of underlined H
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),n2)*sn2;
fp1=par{1}.*te;
%
fp=par{4}.*vet;
te=fft(fp);
te=ifft(te(1:m),n2)*sn2;
fp2=par{3}.*te;
%
fp3=par{5}.*vet;
%
te=fft(fp1-fp2);
fou=te(1:n);
te=fft(fp3);
fou(1:m)=fou(1:m)+te(1:m);
%

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

%
function fou=molthdotra(par,vet,vet1)
%         computes the product of underlined H^T by a vector x 
%         (not appearing among the parameters) of size n
%         vet is the transformed vector of size 2n of x
%         vet1 is the transformed vector of size 2m of the 
%         first m components of x
%         par contains the transformed d-vectors of underlined H
global m
%
m2=2*m;
sm2=sqrt(m2);
%
fp=conj(par{1}).*vet;
te=fft(fp);
te=ifft(te(1:m),m2)*sm2;
fp1=conj(par{2}).*te;
%
fp=conj(par{3}).*vet;
te=fft(fp);
te=ifft(te(1:m),m2)*sm2;
fp2=conj(par{4}).*te;
%
fp3=conj(par{5}).*vet1;
%
te=fft(fp1-fp2+fp3);
fou=te(1:m);
%

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

%
function [foua,foub]=molthh(par,vet)
%         computes the product of overlined H by a vector of size m 
%         whose transformed vector of size 2m is vet 
%         and the product of underlined H by the same vector
%         par contains the transformed d-vectors of overlined H
%         and underlined H
global m n 
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
%
fp=par{3}.*vet;
te=fft(fp);
te1=ifft(te(1:m),m2)*sm2;
fp1a=par{1}.*te1;
te1=ifft(te(1:m),n2)*sn2;
fp1b=par{2}.*te1;
%
fp=par{6}.*vet;
te=fft(fp);
te1=ifft(te(1:m),m2)*sm2;
fp2a=par{4}.*te1;
te1=fft(fp1a-fp2a);
foua=te1(1:m);
te1=ifft(te(1:m),n2)*sn2;
fp2b=par{5}.*te1;
%
fp3=par{7}.*vet;
%

te=fft(fp1b-fp2b);
foub=te(1:n);
te=fft(fp3);
foub(1:m)=foub(1:m)+te(1:m);
%

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

%
function fou=molthup(par,vet)
%         computes the product of overlined H by a vector of size m
%         whose transformed vector of size 2m is vet
%         par contains the transformed d-vectors of overlined H

global m
%
m2=2*m;
sm2=sqrt(m2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(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;
%
te=fft(fp1-fp2);
fou=te(1:m);
%

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

%
function fou=moltlk(par,vet,vet1)
%         computes the product of L by a vector of size n 
%         whose transformed vector of size 2n is vet 
%         plus the product of K by a vector of size m
%         whose transformed vector of size 2m is vet1 
%         par contains the transformed d-vectors of L and K
global m
%
m2=2*m;
sm2=sqrt(m2);
%
fp=par{2}.*vet;
te=fft(fp);
te=ifft(te(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;
%
fp3=par{5}.*vet;
%
%
fp=par{6}.*vet1;
te=fft(fp);
te=ifft(te(1:m),m2)*sm2;
fp4=par{1}.*te;
%
fp=par{7}.*vet1;
te=fft(fp);
te=ifft(te(1:m),m2)*sm2;
fp5=par{3}.*te;
%
fp6=par{8}.*vet1;
%
te=fft(-fp1+fp2-fp4-fp5+fp6);
fou=te(1:m);
te=fft(fp3);
fou(1:m)=fou(1:m)-te(1:m);
%

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

%
function [foua,foub]=moltpim(par,vet,vet1)
%         computes the product of P^-1 by a vector x 
%         (not appearing among the parameters) of size m
%         and the product of M by the same vector
%         vet is the transformed vector of size 2m of x
%         vet1 is the transformed vector of size 2n of x
%         par contains the transformed d-vectors of P^-1 and of M
global m n
%
n2=2*n;
sn2=sqrt(n2);
m2=2*m;
sm2=sqrt(m2);
%
fp=par{5}.*vet;
te=fft(fp);
te1=ifft(te(1:m),m2)*sm2;
fp1a=par{3}.*te1;
te1=ifft(te(1:m),n2)*sn2;
fp1b=par{4}.*te1;
%
fp=par{8}.*vet;
te=fft(fp);
te1=ifft(te(1:m),m2)*sm2;
fp2a=par{6}.*te1;
te1=ifft(te(1:m),n2)*sn2;
fp2b=par{7}.*te1;
%
fp=par{11}.*vet;
te=fft(fp);
te1=ifft(te(1:m),m2)*sm2;
fp3a=par{9}.*te1;
te1=ifft(te(1:m),n2)*sn2;
fp3b=par{10}.*te1;
%
fp=par{12}.*vet;
fp1=par{1}.*vet;
%
te=fft(fp1+fp1a+fp2a-fp3a-fp);
foua=te(1:m);
%
fp=par{13}.*vet;
fp1=par{2}.*vet1;
%
te=fft(fp1b+fp2b+fp3b+fp1);
foub=te(1:n);
te=fft(fp);
foub(1:m)=foub(1:m)+te(1:m);
%

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

%
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;

