function demo10c
% *** Inner product preserving + compensated summation: composed operators ***
  close all
  format long
  format compact
%
% *** user specified quantities **********
  theta=pi/2     % corner opening angle
  npan=11;       % number of coarse panels
  evec=1;        % a unit vector 
% ****************************************
  T16=Tinit16;
  W16=Winit16;
%
  for nsub=1:100
% *** panel breakpoints and interval lengths in parameter ***
    npanfin=npan+2*nsub;
    sinterfin=zeros(npanfin+1,1);
    sinterfin(1:npan+1)=linspace(0,1,npan+1);
    sinterfindiff=ones(npan+2*nsub,1)/npan;
    for k=1:nsub
      sinterfin(3:end)=sinterfin(2:end-1);
      sinterfin(2)=(sinterfin(1)+sinterfin(2))/2;   
      sinterfindiff(3:end)=sinterfindiff(2:end-1);
      sinterfindiff(2)=sinterfindiff(1)/2;
      sinterfindiff(1)=sinterfindiff(1)/2;   
    end
    sinterfin(end-nsub:end)=1-flipud(sinterfin(1:nsub+1));
    sinterfindiff(end-nsub-1:end)=flipud(sinterfindiff(1:nsub+2));
%
% *** discretization points and weights ***  
    npfin=16*npanfin
    sfin=zeros(npfin,1);
    wfin=zeros(npfin,1);
    for k=1:npanfin
      myind=k*16-15:k*16;
      sdif=sinterfindiff(k)/2;
      sfin(myind)=(sinterfin(k)+sinterfin(k+1))/2+sdif*T16;
      wfin(myind)=W16*sdif;
    end
    zfin=zfunc(sfin,theta) ;
    zpfin=zpfunc(sfin,theta);
    zppfin=zppfunc(sfin,theta);
% *** some extra presicion gained from symmetry ***
    zfin(npfin/2+1:npfin)=conj(flipud(zfin(1:npfin/2)));
    zpfin(npfin/2+1:npfin)=-conj(flipud(zpfin(1:npfin/2)));
    zppfin(npfin/2+1:npfin)=conj(flipud(zppfin(1:npfin/2)));
% *************************************************
    nzfin=-1i*zpfin./abs(zpfin);
    wzpfin=wfin.*zpfin;
    sqrtwfin=sqrt(wfin);
    wzpfinPres=sqrtwfin.*zpfin;
%
% *** The K_fin and M_fin matrices are set up ***
    Kfin=M1Ainit(zfin,zpfin,zppfin,nzfin,sqrtwfin,wzpfinPres,npfin);
    Mfin=M1RinitPres(zfin,zpfin,zppfin,wfin,wzpfin,sqrtwfin,npfin);
%    
% *** Solving main linear system ***
    g=2*sqrtwfin.*real(conj(evec)*nzfin);     
    [rhofin,it]=myGMRES(Mfin,Kfin,g,npfin,100,eps);
    GMRES_iter=it
%    
% *** Post processing ***
    qfin=myinner(rhofin,real(conj(evec)*zfin.*abs(wzpfinPres)))
    qref=1.95243329423584
    myerr=abs(qref-qfin)/abs(qref);
    myerr(myerr<eps)=eps;
%
    set(0,'DefaultAxesFontSize',16)
    set(0,'DefaultAxesLineWidth',1) 
%
    figure(1)
    semilogy(nsub,myerr,'*')
    axis([0 100 1e-16 1e0])
    title('Convergence of {\itq} with mesh refinement','FontSize',16)
    xlabel('Number of subdivisions {\itn}_{sub}','FontSize',16)
    ylabel('Estimated relative error in {\itq}','FontSize',16)  
    axis square
    grid on  
    hold on
%  
    figure(2)
    semilogy(nsub,it,'*')
    axis([0 100 1 100])
    title('GMRES iterations for full convergence','FontSize',16)
    xlabel('Number of subdivisions {\itn}_{sub}','FontSize',16)
    ylabel('Number of iterations ','FontSize',16)  
    axis square
    grid on  
    hold on  
  
  end
  
  function M1=M1Ainit(z,zp,zpp,nz,w,wzp,N)
% *** adjoint of double layer potential ***   
  M1=zeros(N);
  for m=1:N
    M1(:,m)=abs(wzp(m))*real(nz./(z(m)-z));
    M1(m,m)=-w(m)*imag(zpp(m)/zp(m))/2;      
  end
  M1=M1/pi;

  function M1=M1RinitPres(z,zp,zpp,w,wzp,sqrtw,N)
% *** double layer potential ***   
  M1=zeros(N);
  for m=1:N
    M1(:,m)=sqrtw.*imag(wzp(m)./(z(m)-z));
    M1(m,m)=sqrtw(m)*w(m)*imag(zpp(m)/zp(m))/2;
  end
  M1=M1/pi;
  
  function zout=zfunc(s,theta)
  zout=sin(pi*s).*exp(1i*theta*(s-0.5));

  function zpout=zpfunc(s,theta)
  zpout=(pi*cos(pi*s)+1i*theta*sin(pi*s)).*exp(1i*theta*(s-0.5));
  
  function zppout=zppfunc(s,theta)
  zppout=(2i*pi*theta*cos(pi*s)-(theta^2+pi^2)*sin(pi*s)).* ...
	 exp(1i*theta*(s-0.5));  

  function [x,it]=myGMRES(A,B,b,n,m,tol)
% *** GMRES with low-threshold stagnation control ***
  V=zeros(n,m+1);
  H=zeros(m);
  cs=zeros(m,1);
  sn=zeros(m,1);
  bnrm2=mynorm(b);
  V(:,1)=b/bnrm2;
  s=bnrm2*eye(m+1,1);
  for it = 1:m                                  
    it1=it+1;                                   
    w=mymatvec(A,mymatvec(B,V(:,it)));
    for k = 1:it
      H(k,it)=myinner(V(:,k),w);
      w=w-H(k,it)*V(:,k);
    end
    H(it,it)=H(it,it)+1;
    wnrm2=mynorm(w);
    V(:,it1)=w/wnrm2;
    for k=1:it-1                                
      temp     = cs(k)*H(k,it)+sn(k)*H(k+1,it);
      H(k+1,it)=-sn(k)*H(k,it)+cs(k)*H(k+1,it);
      H(k,it)  = temp;
    end
    [cs(it),sn(it)]=rotmat(H(it,it),wnrm2);     
    H(it,it)= cs(it)*H(it,it)+sn(it)*wnrm2;
    s(it1) =-sn(it)*s(it);                      
    s(it)  = cs(it)*s(it);                         
    myerr=abs(s(it1))/bnrm2;
    if (myerr<=tol)|(it==m)                     
      predicted_residual=myerr
%      y=triu(H(1:it,1:it))\s(1:it);            
      y=mybacksubst(triu(H(1:it,1:it)),s(1:it));             
      x=fliplr(V(:,1:it))*flipud(y);
      true_residual=mynorm(x+mymatvec(A,mymatvec(B,x))-b)/bnrm2
      break;
    end
  end

  function [c,s]=rotmat(a,b)
  if  b==0
    c=1;
    s=0;
  elseif abs(b)>abs(a)
    temp=a/b;
    s=1/sqrt(1+temp^2);
    c=temp*s;
  else
    temp=b/a;
    c=1/sqrt(1+temp^2);
    s=temp*c;
  end

  function b=mymatvec(A,x)
  [m,n]=size(A);
  b=zeros(m,1);
  for k=1:m
    b(k)=myinner(A(k,:),x.');
  end
  
  function c=mynorm(a)
  c=sqrt(myinner(a,a));

  function c=myinner(a,b)
  c=compsum(a.*b);

  function y=mybacksubst(U,y)
  n=length(y);
  for i=n:-1:1
    y(i)=(y(i)-compsum(U(i,i+1:n)'.*y(i+1:n)))/U(i,i);
  end
  
  function sumout=compsum(x)
  n=length(x);
  sumout=0;
  err=0;
  for k=1:n
    temp=sumout;
    y=x(k)+err;
    sumout=temp+y;
    err=(temp-sumout)+y;
  end
  sumout=sumout+err;
  
  function T=Tinit16
% *** 16-point Gauss-Legendre nodes ***  
  T=zeros(16,1);
  T( 1)=-0.989400934991649932596154173450332627;
  T( 2)=-0.944575023073232576077988415534608345;
  T( 3)=-0.865631202387831743880467897712393132;
  T( 4)=-0.755404408355003033895101194847442268;
  T( 5)=-0.617876244402643748446671764048791019;
  T( 6)=-0.458016777657227386342419442983577574;
  T( 7)=-0.281603550779258913230460501460496106;
  T( 8)=-0.095012509837637440185319335424958063;
  T( 9)= 0.095012509837637440185319335424958063;
  T(10)= 0.281603550779258913230460501460496106;
  T(11)= 0.458016777657227386342419442983577574;
  T(12)= 0.617876244402643748446671764048791019;
  T(13)= 0.755404408355003033895101194847442268;
  T(14)= 0.865631202387831743880467897712393132;
  T(15)= 0.944575023073232576077988415534608345;
  T(16)= 0.989400934991649932596154173450332627;

  function W=Winit16
% *** 16-point Gauss-Legendre weights ***  
  W=zeros(16,1); 
  W( 1)= 0.027152459411754094851780572456018104;
  W( 2)= 0.062253523938647892862843836994377694;
  W( 3)= 0.095158511682492784809925107602246226;
  W( 4)= 0.124628971255533872052476282192016420;
  W( 5)= 0.149595988816576732081501730547478549;
  W( 6)= 0.169156519395002538189312079030359962;
  W( 7)= 0.182603415044923588866763667969219939;
  W( 8)= 0.189450610455068496285396723208283105;
  W( 9)= 0.189450610455068496285396723208283105;
  W(10)= 0.182603415044923588866763667969219939;
  W(11)= 0.169156519395002538189312079030359962;
  W(12)= 0.149595988816576732081501730547478549;
  W(13)= 0.124628971255533872052476282192016420;
  W(14)= 0.095158511682492784809925107602246226;
  W(15)= 0.062253523938647892862843836994377694;
  W(16)= 0.027152459411754094851780572456018104;