% ch0_w2.m --------- Chap 10  GMRES
%    Illustrate W and W' as in A1=W*A*W' and |A M - W'|_F
%Usage: ch0_w2(iprob,m,lev)
%e.g.   ch0_w2(1,4,3)  % for prob=1, Daub 4 and 3-levels (implicit)

function ch0_w2(iprob,order,lev)
if nargin<1, L=10, n=2^L; A=cz(n,2); else
  eval(['load matrix' num2str(iprob)]),
end
  n=size(A,1);
if nargin<2, lev= floor(log(n)/log(2))-7; end
warning off MATLAB:flops:UnavailableFunction

rand('state',4210),x=zeros(n,1);
b = rand(n,1); if iprob==3;b = b+i*rand(n,1); end

W=filters(n,order,lev); tol=0.1; method=1;  %% setup 4 ways to make patt
%%patt=fwts(spdiags(b,0,n,n), order, 0, 1, lev); % For pettern only  (a)
%%patt=spdiags([b b b],[-1 0 1],n,n)*W'; %% Make up a pattern        (b)

restart=50; As=A; % Sparse?
if iprob==3
 TOL=mean(mean(abs(A)));As=(abs(A)>TOL*1.8) .* A;  % threshold for dense A
end
if iprob==0
 restart=500;
end
  nzmax=ceil(nnz(As)/n)*15; if iprob>=1, nzmax=ceil(nnz(As)/n)*2; end

  patt=fwts(As, order, 0, 1, lev); % For pettern only          (c)
% patt=W';  % Aim to fill in band-width of nnzmax  (d)

Adpa=0
M = spai2(As,nzmax,tol,method,W',patt,Adpa); %patt=W' not yet the half finger
figure;spy(patt),title('patt matrix')
figure;spy(W'),title('W^T')
figure;spy(M),title('M')
figure
drawnow

 TOL=1.0E-6; maxit=ceil(n/restart);

  [x,res]=gmrest_k(A,b,restart,W,M,TOL,maxit,b);  %%for  P2*A*P1*y = P2*b
% [x,res]=gmrest_k(A,b,restart,[],[],TOL,maxit,b);  %%for  P2*A*P1*y = P2*b

res=res/res(1); semilogy(res,'-ro'); hold on;git=length(res);
Len=max(length(res),40); semilogy([1 Len], [1 1.0E-6],'g:'); hold off
s=sprintf('%d Iterations of a %d-level method for N=%d (%s)',...
  git,lev,n, Matrix);       axis([0 Len 1.0E-6 1]); hold off
s=sprintf('Residuals for Prob %d',iprob);  ylabel(s)
  xlabel(['Test matrix' num2str(iprob) ': Cgce steps (M' ...
          num2str(method) ', levs=' num2str(lev) ')']); 
s=text(20,0.5,Matrix);set(s,'Fontsize',12,'Color','b')
%__________________________________________________________________subprogram
%%filters.m --- for Daubechies' order m filters' matrix
%%
%%Usage: W = filters(N,m,lev)      for matrix A of N x N (order m) in lev levels
%     or W = filters(N,m)          for matrix A of N x N (order m) in max levels
%
%%e.g.  A2=fwts(A,4,1,1,7); W=filters(128,4,7); % Daub 4 in 7 levels
%       A1=W*A*W';   Differ = norm(A1-A2,1)
%       spy(W)                     % ___________________________________________
function W = filters(N,m,levs)
  if nargin<2, help filters; return; end
  W=speye(N); n=N; L_max=floor(log(N)/log(2));
  if nargin<3,levs=L_max;end
  if levs>L_max,error(['Input levs = ' num2str(levs) ' too large!']); end
for lev=1:levs
    Gap=2^(lev-1); nH=floor(n/2);
    row=1:nH; col=1:n; W_L=speye(N);  K=w_matrix(m,n);
 W_L(row,col)   =K(1:2:n,:);  % Mr Sums (c_j)
 W_L(row+nH,col)=K(2:2:n,:);  % Mr Difs (d_j)
 W=W_L*W;                     % Form products
 n = nH;                      % Next level dimension
end
return

function W = w_matrix(m,n)   % used in fwts.m
if nargin<1, help filters; return; end
   p = find([2 4 6 8 16] == m);
if m<2, W=speye(n); return; end
 %  Note these DIAG repeating patterns:   (0,c_{m-1}), (c_{m-2} d_{m-1})
 %   (c_{m-4} d_{m-3}) (c_{m-6} d_{m-5}) ... (c_{0} d_1) and (c_{m-1} 0)
if m==2
c1=sqrt(2)/2; c0=sqrt(2)/2; 
d1=zeros(n,1); d1(1:2:n)=-c0*ones(n/2,1);
d3=c0*ones(n,1);
d5=zeros(n,1); d5(2:2:n)= c0*ones(n/2,1);
W=spdiags([d1 d3 d5],-1:1,n,n);
elseif m==4      %====================================================
c3=sqrt(2)*(1+sqrt(3))/8; c2=sqrt(2)*(3+sqrt(3))/8;
c1=sqrt(2)*(3-sqrt(3))/8; c0=sqrt(2)*(1-sqrt(3))/8;

d1=zeros(n,1); d1(1:2:n)=c0*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c3*ones(n/2,1); d2(2:2:n)=-c1*ones(n/2,1);
d3=c2*ones(n,1);
d4=zeros(n,1); d4(1:2:n)=c1*ones(n/2,1); d4(2:2:n)=-c3*ones(n/2,1);
d5=zeros(n,1); d5(2:2:n)=c0*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5],-1:3,n,n);
W(n-1:n,1:2)=[c1 c0; c2 -c3];
elseif m==6      %====================================================
c0= 3.326705529500826d-01; c1= 8.068915093110927d-01;
c2= 4.598775021184915d-01; c3=-1.350110200102546d-01;
c4=-8.544127388202666d-02; c5= 3.522629188570955d-02;

d1=zeros(n,1); d1(1:2:n)=c5*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c0*ones(n/2,1); d2(2:2:n)=-c4*ones(n/2,1);
d3=zeros(n,1); d3(1:2:n)=c3*ones(n/2,1); d3(2:2:n)= c1*ones(n/2,1);
d4=zeros(n,1); d4(1:2:n)=c2*ones(n/2,1); d4(2:2:n)=-c2*ones(n/2,1);
d5=zeros(n,1); d5(1:2:n)=c1*ones(n/2,1); d5(2:2:n)= c3*ones(n/2,1);
d6=zeros(n,1); d6(1:2:n)=c4*ones(n/2,1); d6(2:2:n)=-c0*ones(n/2,1);
d7=zeros(n,1);                           d7(2:2:n)= c5*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5 d6 d7],-1:5,n,n);
W(n-3:n-2,1:2)=[c4 c5; c1 -c0];
W(n-1:n,3:4)=[c4 c5; c1 -c0];
W(n-1:n,1:2)=[c2 c3; c3 -c2];
elseif m==8      %==================================================== 
c0=-7.5765714789357d-2; c1=-2.9635527645960d-2;
c2= 4.9761866763256d-1; c3= 8.0373875180539d-1;
c4= 2.9785779560561d-1; c5=-9.9219543576956d-2;
c6=-1.2603967262264d-2; c7= 3.2223100604078d-2;

d1=zeros(n,1); d1(1:2:n)=c7*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c0*ones(n/2,1); d2(2:2:n)=-c6*ones(n/2,1);
d3=zeros(n,1); d3(1:2:n)=c5*ones(n/2,1); d3(2:2:n)= c1*ones(n/2,1);
d4=zeros(n,1); d4(1:2:n)=c2*ones(n/2,1); d4(2:2:n)=-c4*ones(n/2,1);
d5=c3*ones(n,1); 
d6=zeros(n,1); d6(1:2:n)=c4*ones(n/2,1); d6(2:2:n)=-c2*ones(n/2,1);
d7=zeros(n,1); d7(1:2:n)=c1*ones(n/2,1); d7(2:2:n)= c5*ones(n/2,1);
d8=zeros(n,1); d8(1:2:n)=c6*ones(n/2,1); d8(2:2:n)=-c0*ones(n/2,1);
d9=zeros(n,1);                           d9(2:2:n)= c7*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5 d6 d7 d8 d9],-1:7,n,n);
W(n-5:n-4,1:2)=[c6 c7; c1 -c0];
W(n-3:n-2,3:4)=[c6 c7; c1 -c0];
W(n-1:n  ,5:6)=[c6 c7; c1 -c0];
W(n-3:n-2,1:2)=[c4 c5; c3 -c2];
W(n-1:n,3:4)=[c4 c5; c3 -c2];
W(n-1:n,1:2)=[c2 c3; c5 -c4];
elseif m==8      %==================================================== 
    c0  = 1.8899503329009d-03; c1  =-3.0292051455166d-04;
    c2  =-1.4952258336794d-02; c3  = 3.8087520140604d-03;
    c4  = 4.9137179673481d-02; c5  =-2.7219029916816d-02;
    c6  =-5.1945838107879d-02; c7  = 3.6444189483599d-01;
    c8  = 7.7718575169981d-01; c9  = 4.8135965125924d-01;
    c10 =-6.1273359067914d-02; c11 =-1.4329423835107d-01;
    c12 = 7.6074873252854d-03; c13 = 3.1695087810348d-02;
    c14 =-5.4213233163559d-04; c15 =-3.3824159513597d-03;
d1=zeros(n,1); d1(1:2:n)=c15*ones(n/2,1);
d2=zeros(n,1); d2(1:2:n)=c0 *ones(n/2,1); d2(2:2:n)=-c14*ones(n/2,1);
d3=zeros(n,1); d3(1:2:n)=c13*ones(n/2,1); d3(2:2:n)= c1 *ones(n/2,1);
d4=zeros(n,1); d4(1:2:n)=c2 *ones(n/2,1); d4(2:2:n)=-c12*ones(n/2,1);
d5=zeros(n,1); d5(1:2:n)=c11*ones(n/2,1); d5(2:2:n)=c3  *ones(n/2,1);
d6=zeros(n,1); d6(1:2:n)=c4 *ones(n/2,1); d6(2:2:n)=-c10*ones(n/2,1);
d7=zeros(n,1); d7(1:2:n)=c9 *ones(n/2,1); d7(2:2:n)= c5 *ones(n/2,1);
d8=zeros(n,1); d8(1:2:n)=c6 *ones(n/2,1); d8(2:2:n)=-c8 *ones(n/2,1);
d9=c7*ones(n,1); 
d10=zeros(n,1); d10(1:2:n)=c8 *ones(n/2,1); d10(2:2:n)=-c6*ones(n/2,1);
d11=zeros(n,1); d11(1:2:n)=c5 *ones(n/2,1); d11(2:2:n)=c9 *ones(n/2,1);
d12=zeros(n,1); d12(1:2:n)=c10*ones(n/2,1); d12(2:2:n)=-c4*ones(n/2,1);
d13=zeros(n,1); d13(1:2:n)=c3 *ones(n/2,1); d13(2:2:n)=c11*ones(n/2,1);
d14=zeros(n,1); d14(1:2:n)=c12*ones(n/2,1); d14(2:2:n)=-c2*ones(n/2,1);
d15=zeros(n,1); d15(1:2:n)=c1 *ones(n/2,1); d15(2:2:n)=c13*ones(n/2,1);
d16=zeros(n,1); d16(1:2:n)=c14*ones(n/2,1); d16(2:2:n)=-c0*ones(n/2,1);
d17=zeros(n,1);                             d17(2:2:n)=c15*ones(n/2,1);

W=spdiags([d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13 d14 d15 d16 d17],-1:15,n,n);
W=W+spdiags([d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13 d14 d15 d16 d17],-n+1:-n+15,n,n);
end
