function demo6
  close all
  format long
  format compact
  NA=16;
  [IP,IPW,TX]=gaussPAX(NA);
  set(0,'DefaultAxesFontSize',16)
  set(0,'DefaultAxesLineWidth',1) 
  for nsub=1:100
    nsub
    [IPP,IPPW]=Poper2(NA,IP,IPW,nsub);
    Pnsub=zeros((2*nsub+4)*NA,4*NA);
    Pnsub(            1  :           NA,     1:  NA)=eye(NA);
    Pnsub(         NA+1  :  (nsub+2)*NA,  NA+1:2*NA)=IPP;
    Pnsub((nsub+2)*NA+1  :(2*nsub+3)*NA,2*NA+1:3*NA)=flipud(fliplr(IPP));
    Pnsub((2*nsub+3)*NA+1:(2*nsub+4)*NA,3*NA+1:4*NA)=eye(NA);   
    figure(1)
    spy(Pnsub)
    title('P_{nsub}')
    drawnow
%
    IQQ=Qoper(NA,TX,nsub);    
    Qnsub=zeros(4*NA,(2*nsub+4)*NA);
    Qnsub(     1:  NA,            1  :           NA)=eye(NA);
    Qnsub(  NA+1:2*NA,         NA+1  :  (nsub+2)*NA)=IQQ;
    Qnsub(2*NA+1:3*NA,(nsub+2)*NA+1  :(2*nsub+3)*NA)=flipud(fliplr(IQQ));
    Qnsub(3*NA+1:4*NA,(2*nsub+3)*NA+1:(2*nsub+4)*NA)=eye(NA);   
    figure(2)
    spy(Qnsub)
    title('Q_{nsub}')
    drawnow
    err_QP=norm(Qnsub*Pnsub-eye(4*NA))
%
    PWnsub=zeros((2*nsub+4)*NA,4*NA);
    PWnsub(            1  :           NA,     1:  NA)=eye(NA);
    PWnsub(         NA+1  :  (nsub+2)*NA,  NA+1:2*NA)=IPPW;
    PWnsub((nsub+2)*NA+1  :(2*nsub+3)*NA,2*NA+1:3*NA)=flipud(fliplr(IPPW));
    PWnsub((2*nsub+3)*NA+1:(2*nsub+4)*NA,3*NA+1:4*NA)=eye(NA);   
    err_PWTP=norm(PWnsub'*Pnsub-eye(4*NA))
%
    figure(3)
    semilogy(nsub,err_PWTP,'*',nsub,err_QP,'o')
    axis([0 100 1e-16 1e0])
    title(['Verfification of {\bfP}_W^T{\bfP}={\bfI} and ', ... 
	   '{\bfQ}{\bfP}={\bfI}'],'FontSize',16)
    hl=legend('${\bf P}_{Wn_{\rm sub}}^T{\bf P}_{n_{\rm sub}}$', ...
	      '${\bf Q}_{n_{\rm sub}}{\bf P}_{n_{\rm sub}}$');
    set(hl,'Interpreter','latex')    
    xlabel('Number of subdivisions {\itn}_{sub}','FontSize',16)
    ylabel('Estimated relative error','FontSize',16)
    axis square
    grid on  
    hold on
  end
  
  function [IPP,IPPW]=Poper(N,ndiv)
  digits(36)
  k=vpa(1:N-1);
  B=k./sqrt(4*k.*k-1);
  [V,D]=eig(diag(B,1)+diag(B,-1));
  DD=diag(D);
  [~,myind0]=sort(double(DD));
  T=DD(myind0);
  W=2*(V(1,myind0).*V(1,myind0))';
  A=vpa(ones(N));
  AA=vpa(ones(2*N,N));
  T2=[T-1;T+1]/2;
  W2=[W;W]/2;
  for k=2:N
    A(:,k)=A(:,k-1).*T;
    AA(:,k)=AA(:,k-1).*T2;   
  end
  IP=AA/A;
  IPP=vpa(eye(N));
  for k=1:ndiv
    NPB=vpa([eye(N*(k-1)) zeros(N*(k-1),N);zeros(2*N,N*(k-1)) IP]);
    IPP=NPB*IPP;
  end
  IPP=double(IPP);
  IPW=IP.*(W2*(1./W)');
  IPPW=vpa(eye(N));
  for k=1:ndiv
    NPB=vpa([eye(N*(k-1)) zeros(N*(k-1),N);zeros(2*N,N*(k-1)) IPW]);
    IPPW=NPB*IPPW;
  end
  IPPW=double(IPPW);

  function IQQ=Qoper(N,T,ndiv)
  digits(36)
  if ndiv==0
    IQQ=vpa(eye(N));
  else
    T1=(T+1)/2;
    IQQ=vpa(zeros(N,N*(ndiv+1)));
    mysum=0;
    T1d=double(T1);
    for idiv=1:ndiv
      Tnow=T1(T1d>(1-1/2^(idiv-1)) & T1d<=(1-1/2^idiv));
      n=length(Tnow);
      if n>0
	TShort=(Tnow-(1-1/2^(idiv-1)))*2^(idiv+1)-1;
	A=vpa(ones(N));
	AS=vpa(ones(n,N));
	for k=1:N-1
	  A(:,k+1)=A(:,k).*T;
	  AS(:,k+1)=AS(:,k).*TShort;
	end
	IQQ(mysum+1:mysum+n,(idiv-1)*N+1:idiv*N)=AS/A;
	mysum=mysum+n;
      end
    end
    Tnow=T1(T1d>(1-1/2^ndiv));
    n=length(Tnow);
    if n>0
      TShort=(Tnow-(1-1/2^ndiv))*2^(ndiv+1)-1;
      A=vpa(ones(N));
      AS=vpa(ones(n,N));
      for k=1:N-1
	A(:,k+1)=A(:,k).*T;
	AS(:,k+1)=AS(:,k).*TShort;
      end
      IQQ(mysum+1:mysum+n,ndiv*N+1:(ndiv+1)*N)=AS/A;
    end
  end
  IQQ=double(IQQ);

  function [IPP,IPPW]=Poper2(N,IP,IPW,ndiv);
  IPP=eye(N);
  for k=1:ndiv
    NPB=[eye(N*(k-1)) zeros(N*(k-1),N);zeros(2*N,N*(k-1)) IP];
    IPP=NPB*IPP;
  end
  IPP=double(IPP);
  IPPW=eye(N);
  for k=1:ndiv
    NPB=[eye(N*(k-1)) zeros(N*(k-1),N);zeros(2*N,N*(k-1)) IPW];
    IPPW=NPB*IPPW;
  end
  IPPW=double(IPPW);
  
  function [IP,IPW,TA]=gaussPAX(NA)
  digits(36)
  k=vpa(1:NA-1);
  B=k./sqrt(4*k.*k-1);
  [V,D]=eig(diag(B,1)+diag(B,-1));
  DD=diag(D);
  [~,myind0]=sort(double(DD));
  TA=DD(myind0);
  WA=2*(V(1,myind0).*V(1,myind0))';
%
  AMAT=vpa(ones(NA));
  AAMAT=vpa(ones(2*NA,NA));
  TAA=[(TA-1)/2;(TA+1)/2];
  WAA=[WA;WA]/2;
  for k=2:NA
    AMAT(:,k)=AMAT(:,k-1).*TA;
    AAMAT(:,k)=AAMAT(:,k-1).*TAA;         
  end
%
  IP=AAMAT/AMAT;
  IPW=IP.*(WAA*(1./WA)');
  IP=double(IP);
  IPW=double(IPW)