function demo9
  close all
  format long
  format compact
%
% *** user specified quantities **********
  lambda=0.999   % parameter lambda
  theta=pi/2     % wedge opening angle
% ****************************************
  T16=Tinit16;
  W16=Winit16;
  [IP,IPW]=IPWinit(T16,W16);
%
  Pbc=zeros(96,64);
  Pbc( 1:16, 1:16)=eye(16);
  Pbc(17:48,17:32)=IP;
  Pbc(49:96,33:64)=flipud(fliplr(Pbc(1:48,1:32)));
%
  PWbc=zeros(96,64);
  PWbc( 1:16, 1:16)=eye(16);
  PWbc(17:48,17:32)=IPW;
  PWbc(49:96,33:64)=flipud(fliplr(PWbc(1:48,1:32)));
%
% *** Initializer for the recursion ***
  [R0,err0]=initializeR(W16,T16,lambda,theta,PWbc,Pbc);
  [R00,err00]=initializeN(W16,T16,lambda,theta,PWbc,Pbc);
  agreement=norm(R00-R0)/norm(R0)
%
  set(0,'DefaultAxesFontSize',16)
  set(0,'DefaultAxesLineWidth',1) 
  figure(1)
  semilogy(1:length(err0),err0,'ro',1:length(err00),err00,'b*')
  axis([0 100 1e-16 1e0])
  title('Fixed-point iteration versus Newton''s method','FontSize',16)
  xlabel('Number of iteration steps','FontSize',16)
  ylabel('Estimated relative error in {\bfR}_*','FontSize',16)  
  legend('Fixed-point','Newton')
  grid on  
  
  function [R,errvec]=initializeR(W16,T16,lambda,theta,PWbc,Pbc)
  myerr=1;
  sloc=[T16/4+0.25;T16/4+0.75;T16/2+1.5];
  wloc=[W16/4;W16/4;W16/2];
  zloc=sloc*exp(-1i*theta/2);
  zploc=ones(48,1)*exp(-1i*theta/2);
  zpploc=zeros(48,1);
  zloc=[conj(flipud(zloc));zloc];
  zploc=[-conj(flipud(zploc));zploc];
  zpploc=[conj(flipud(zpploc));zpploc];
  wloc=[flipud(wloc);wloc];
  nzloc=-1i*zploc./abs(zploc);
  wzploc=wloc.*zploc;
  Kloc=M1Ainit(zloc,zploc,zpploc,nzloc,wloc,wzploc,96);
  MAT=eye(96)+lambda*Kloc;
  R=inv(MAT(17:80,17:80));
  iter=0;
  while myerr>eps
    Rold=R;
    R=SchurBana(PWbc,Pbc,MAT,R);
    myerr=norm(R-Rold,'fro')/norm(R,'fro');
    iter=iter+1;
    errvec(iter)=myerr;
  end
  Init_iter_fixed=iter
  
  function [R,errvec]=initializeN(W16,T16,lambda,theta,PWbc,Pbc)
  myerr=1;
  sloc=[T16/4+0.25;T16/4+0.75;T16/2+1.5];
  wloc=[W16/4;W16/4;W16/2];
  zloc=sloc*exp(-1i*theta/2);
  zploc=ones(48,1)*exp(-1i*theta/2);
  zpploc=zeros(48,1);
  zloc=[conj(flipud(zloc));zloc];
  zploc=[-conj(flipud(zploc));zploc];
  zpploc=[conj(flipud(zpploc));zpploc];
  wloc=[flipud(wloc);wloc];
  nzloc=-1i*zploc./abs(zploc);
  wzploc=wloc.*zploc;
  Kloc=M1Ainit(zloc,zploc,zpploc,nzloc,wloc,wzploc,96);
  MAT=eye(96)+lambda*Kloc;
  R=inv(MAT(17:80,17:80));
  iter=0;
  while myerr>eps
    [A,B,C]=SchurBanaL2(PWbc,Pbc,MAT,R);
    myerr=norm(C,'fro')/norm(R,'fro');
    R=R+GMRESylv(A,B,C,64,40,eps);
    iter=iter+1;
    errvec(iter)=myerr;  
  end
  Init_iter_Newton=iter

  function A=SchurBana(PW,P,K,A)
  starL=17:80;
  circL=[1:16 81:96];
  starS=17:48;
  circS=[1:16 49:64];
  VA=K(circL,starL)*A;
  PTA=PW(starL,starS)'*A;
  PTAU=PTA*K(starL,circL);
  DVAUI=inv(K(circL,circL)-VA*K(starL,circL));
  DVAUIVAP=DVAUI*(VA*P(starL,starS));
  A(starS,starS)=PTA*P(starL,starS)+PTAU*DVAUIVAP;
  A(circS,circS)=DVAUI;
  A(circS,starS)=-DVAUIVAP;
  A(starS,circS)=-PTAU*DVAUI;
      
  function [K1,K2,K3]=SchurBanaL2(PWbc,Pbc,K,A,NA)
  starL=17:80;
  circL=[1:16 81:96];
  starS=17:48;
  circS=[1:16 49:64];
  VA=K(circL,starL)*A;
  VAP=VA*Pbc(starL,starS);
  PTA=PWbc(starL,starS).'*A;
  PTAU=PTA*K(starL,circL);
  DVAUI=inv(K(circL,circL)-VA*K(starL,circL));
  DVAUIVAP=DVAUI*VAP;
  UDVAUI=K(starL,circL)*DVAUI;
  DVAUIV=DVAUI*K(circL,starL);
%
  K3=zeros(64);
  K3(starS,starS)=PTA*Pbc(starL,starS)+PTAU*DVAUIVAP;
  K3(circS,circS)=DVAUI;
  K3(circS,starS)=-DVAUIVAP;
  K3(starS,circS)=-PTAU*DVAUI;
  K3=K3-A;
%
  K1=zeros(64);
  K1(starS,:)=PWbc(starL,starS).'+PTAU*DVAUIV;
  K1(circS,:)=-DVAUIV;
%   
  K2=zeros(64);
  K2(:,starS)=Pbc(starL,starS)+UDVAUI*VAP;
  K2(:,circS)=-UDVAUI;

  function X=GMRESylv(C,D,E,n,m,tol)
  n2=n*n;
  V=zeros(n2,m);
  H=zeros(m);
  cs=zeros(m,1);
  sn=zeros(m,1);
  bnrm2=norm(E(:));
  V(:,1)=E(:)/bnrm2;
  s=bnrm2*eye(m+1,1);
  for it = 1:m        
    it1=it+1;                                   
    Vmat=shapeup(V(:,it),n);
    wmat=-C*Vmat*D; 
    w=wmat(:);
    for k = 1:it
      H(k,it)=V(:,k)'*w;
      w=w-H(k,it)*V(:,k);
    end
    H(it,it)=H(it,it)+1;
    wnrm2=norm(w);
    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);                         
    error=abs(s(it1))/bnrm2;
    if (error<=tol)||(it==m)                   
      y=triu(H(1:it,1:it))\s(1:it);
      X=shapeup(fliplr(V(:,1:it))*flipud(y),n);
      GMRESylv_iter=it;
      break;
    end
    V(:,it1)=w/wnrm2;   
  end

  function B=shapeup(b,n)
  B=zeros(n);
  B(:)=b;
  
  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 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 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 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;

  function [IP,IPW]=IPWinit(T16,W16)
  A=ones(16);
  AA=ones(32,16);
  T32=[(T16-1)/2;(T16+1)/2];
  W32=[W16;W16]/2;
  for k=2:16
    A(:,k)=A(:,k-1).*T16;
    AA(:,k)=AA(:,k-1).*T32;         
  end
  IP=AA/A;
  IPW=IP.*(W32*(1./W16)');