%============================================================= 2D 
%ch13.m adapted from
%     cir0.m --- illustration of Toeplitz and Circulant Matrix
%============================================================= BTTB
warning off MATLAB:flops:UnavailableFunction
format compact
if ~exist('n'),n=16, end

if ~exist('iout'),iout=0; end
if n>10, iout=0, end  %% Suppress large output!

%%%%%%%=========================================2D BTTB======================
rand('state',8)
R = round(rand(2*n-1)*9+0.5);    % Toeplitz root in columns for col 1
n2=2*n;  

rand('state',6)
x = round(rand(n)*12+0.5);    % v for w=Kv

%_____Direct Toplitz BTTB__________|____No symmetry_____ Method 1

disp('1/ DIRECT : BTTB direct method (the orginal problem) --------')

t1=cputime;
  T=BTTB(R,0); x1 = T*reshape(x,n*n,1);  x1 = reshape(x1,n,n); 
t1=cputime-t1;

if iout==1
  disp('Here is the root R (col 1''s)')
  mat_pr0(R,1,1,0)
  disp('Here is the T=BTTB:')
  mat_pr0(T,n,1,0)
end
  %clear T

%_____Prepare for BCCB extension__________________________Method 2
% Root extension for Toeplitz BTTB to be a larger BCCB
  R_ext = [zeros(1,n2-1); R]; R_ext0 = [zeros(n2,1)  R_ext];
% Added 1 row and 1 col of zeros to make n2 x n2

disp('2/ DIRECT : BCCB direct method (already `twice'' as large) -------')

  R_ext = fftshift(R_ext0);
    C_L = bccb(R_ext); 
  x_ext = zeros(n2,n2); x_ext(1:n,1:n)=x;  % Extend input vector 
 x2L = reshape(x_ext,n2*n2,1);

t2=cputime;
  x2 = C_L*x2L; x20 = reshape(x2,n2,n2); x2 = x20(1:n,1:n);
t2=cputime-t2;

if iout==1
 disp('Root R_ext for larger BCCB .....  ......  ......  ...... ')
 mat_pr0(R_ext,n2,1,0)
 disp('Matrix C_L .... the resulting BCCB (1st 2 block cols)... ')
 mat_pr0(C_L(:,1:max(n,24)),n2,1,0)  
end
  %clear C_L

%_____The power of BCCB done by FTT2 _____________________Method 3

disp('3/ FFT2 :  for BCCB extended method (already `twice'' as large) -------')
t3=cputime;
x30 = real( ifft2( fft2(R_ext,n2,n2) .* fft2(x_ext,n2,n2) ));
t3=cputime-t3;

x3 = x30(1:n,1:n);

%_____Direct Toplitz BTTB______________|____Symmetry_____ Method 4
disp('4/ DIRECT : BTTB direct method (a new symmetric problem) --------')
R=R(1:n,1:n);  %% Reset for a symmetric BTTB (smaller root)
t4=cputime;
T=BTTB(R,1);   x4 = T*reshape(x,n*n,1);  x4 = reshape(x4,n,n); 
t4=cputime-t4;

%_____The power of BCCB done by FTT2 _____________________Method 5

disp('5/ FFT2 :  for BCCB extended method (already `twice'' as large) -------')
r = [1 n:-1:2]; F_ext = [R R(:,r)]; % Horizontal enlargement
       F_ext = [F_ext; F_ext(r,:)]; % Vertical   enlargement
t5=cputime;
x50 = real( ifft2( fft2(F_ext,n2,n2) .* fft2(x_ext,n2,n2) ));
t5=cputime-t5;

x5 = x50(1:n,1:n);

if iout==1
 disp('x (data)  x1 (BTTB)  x2 (BCCB)  x3 (FFT2 BCCB)  x4 (BTTB symm) x5 (FFT2 symm) -- Multiplications:')
 mat_prt(x(:),[x1(:) x2(:) x3(:) x4(:) x5(:)],' |','x ','4 Products: x1 x2 x3 x4 x5',0)
end

diff_12 = norm( x1 - x2, 1); diff_23 = norm( x2 - x3, 1);
                             diff_45 = norm( x5 - x4, 1);
fprintf('\n\tComparisons:  error(x1-x2)_1 = %e\n',diff_12)
fprintf('\t              error(x2-x3)_1 = %e\n',diff_23)
fprintf('\t              error(x4-x5)_1 = %e\n',diff_45)
fprintf('\tCPU for Methods 1-5 = \n\t %f, %f, %f, %f, %f [n=%d]\n',...
                         t1,t2,t3,t4,t5,n)
disp('===============================================================')
