%=========================================================== circ_toep
%circ_toep.m -- Illustration of Toeplitz and Circulant Matrix
%               (calling formating Mfiles: mat_prt.m and mat_prt4.m)
%Tasks (1) compute f=Cu 
%      (2) compute v=Tg
%      (2) solve   f=C\u
%Input vectors: c and r for Toeplitz
%             : c       for Circulant case
%             : u (RHS) for Circulant case
%             : g (RHS) for Toeplitz 
% or don't set c, r, u, g to use default values (to reproduce Examples)
%
%=========================================================== circ_toep
format compact
help circ_toep
if ~exist('mat_prt');
error(['Download or copy mat_prt.m and mat_prt4.m to PATH before run this Mfile: ' mfilename])
end

disp(' TASK (1) :   f=Cu +++++++++++++++++++++++++++++++++++++++++++++++++ ')

disp('Example 2.1================1D Circulant============  circ_toep.m')
if ~exist('iout');iout=1;end, global n_o sep_on; % output control

if ~exist('c') 
   n = 8; rand('seed',0); c = round(randperm(n)); % column
end
   n =length(c); C8 = gallery('circul',c)';
   mat_prt(C8,[],' ','C8 matrix')

if ~exist('u') 
   rand('seed',311192157)
   u = round(rand(n,1)*10+0.5), % test vector for f = C8*u
end

f = C8*u;     %DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDirect --- the hard way
f_FFT=real( ifft( fft(C8(:,1)) .* fft(u) )); %---------- the FFT way 
fprintf('\nFFT_of_Col_1__________\n          FFT_of_u____________\n')
fprintf('                  Point_product____________')
 mat_prt4( [fft(C8(:,1)) fft(u) fft(C8(:,1)).*fft(u)] )
fprintf('\nRHS u ________________\n          f=Cu via FFT________')
 mat_prt4( [u f_FFT] )

disp(' TASK (2) :  v=Tg +++++++++++++++++++++++++++++++++++++++++++++++++ ')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Example 2.2================1D Toeplitz to Circulant========circ_toep.m')

if ~exist('c') | ~exist('r')
   n = 8; rand('seed',0); r = round(randperm(n));
   c = round(rand(1,n)*9+pi); c(1)=r(1);   %  r=row vector, c=col
end
   n =length(c);

T = toeplitz(c , r); disp(['Matrix T size = ' num2str(size(T))]);
 mat_prt(T,[],' ','T8 matrix')
t = [T(1,n:-1:1) T(2:n,1)'];
ce= [t(n:2*n-1) 0 t(1:n-1)];

C16 = gallery('circul',ce)';
 disp(['Matrix C size = ' num2str(size(C16))]);
 mat_prt(C16,[],' ','C16 matrix')
 n_o=1;

if ~exist('g') 
   rand('seed',774845808)
   g = rand(n,1); % new test vector for RHS g:   v = T*g
end

g_tilde = [g' zeros(1,n)]'; %%% Extension to match C dimension 
v = T*g;     %DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDirect --- the hard way
v_tilde = C16*g_tilde; %EEEEEEEEEEEEEEEEEEEEEExtended --- the hard way (phew)
v_FFT=real( ifft( fft(C16(:,1)) .* fft(g_tilde) )); %%%%% the FFT way 

if iout==1
%mat_prt(T,C16, '  V', 'T Toeplitz','C Circulant - See Chapter two',1)
 mat_prt([g v],[g_tilde v_tilde v_FFT], '  |', 'Left g T*g',...
               'g_til v_tilde v_FFT')
fprintf('\nFFT_of_Col_1__________\n          FFT_of_v_tilde______\n')
fprintf('                  Point_product____________')
 mat_prt4( [fft(C16(:,1)) fft(g_tilde) fft(C16(:,1)).*fft(g_tilde)] )

fprintf('\nRHS g ________________\n          v=Tu via FFT________')
 mat_prt4( [g v_FFT(1:n)] )
end
    v_tilde=v_tilde(1:n); v_FFT=v_FFT(1:n);
fprintf(' ========== 1D Toeplitz via circulant extension ==========\n')
fprintf(' M1:direct T_nxn,  M2:direct C_2nx2n,  M3:FTT 2nx1 vector\n')
diff_12 = norm( v - v_tilde, 1)
diff_23 = norm( v_tilde - v_FFT, 1) 

fprintf('\n TASK (3) :   f=C\\u +++++++++++++++++++++++++++++++++++++++++++ \n')

f_direct = C8 \ u

f_FFT=real( ifft( fft(u) ./ fft(C8(:,1)) )); %---------- the FFT way 

fprintf('\nFFT_of_u______________\n          FFT_of_Col_1 (F_n h)______\n')
fprintf('                   Point Division___________')
  mat_prt4( [fft(u) fft(C8(:,1)) fft(u)./fft(C8(:,1))] )
fprintf('\nRHS u ________________\n          f=Cu via FFT________\n')
fprintf('                   f=C\\u directly___________')
 mat_prt4( [u f_FFT f_direct] )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cheers
