% ch3_fmm.m ------- Driver to call "FMM toolbox" from www.madmaxoptics.com
% (1) Create the data for testing
%       P(i) = \sum_{j=1}^N q_j \log | x_i - y_j |
% (2) Compute the product P_d = Aq directly
% (3) Compute the product P_f = Aq by the FMM

function p=use_fmm(N)
   v=version; if str2num(v(1))>5,
      warning off MATLAB:flops:UnavailableFunction, end
   if nargin<1, help ch3_fmm, return, end 
   iout=0; if N<25, iout=1; end %% OK to display
disp('(1) --------- Data creation ------- for (x_j,y_j) and q_j')
   rand('state',25);
   if iout==0
   x = rand(N,1); y = rand(N,1); q = rand(N,1)*N; 
   else % integer for easier display
   x = round(rand(N,1)*N); y = round(rand(N,1)*N);
   q = round(rand(N,1)*N); 
   end
   D = createsrc2d(x,y, q);  % Note there are many other options + models
   figure;
    plot(x,y,'xr',x,y,'go',x,y,'b+');
    t=title(['ch3\_fmm.m -- N =' num2str(N)]);set(t,'color','r')

disp('(2) --------- Direct product ------ for P_d=A*q')
flops(0), t=cputime; A = matrixq(x,y); P_d = A*q; f1=flops;
t1=cputime-t;

disp('(3) --------- The FMM product ----- for P_f=A*q')
flops(0),t=cputime; P_f = fmmcoul2d( D, x,y); t2=cputime-t;
f2=flops;

if iout==1
mat_prt4(A,q,'|',' Matrix A',' RHS q')
mat_prt([x y q], [P_d P_f], ':',' (x y) and q ','P_d direct and P_f FMM')
else  % not meaningful to show cpu for small N
fprintf('Computing p=A*q takes %f seconds by a direct product\n',t1)
fprintf('                  and %f seconds by the FMM (N=%d)\n',t2,N)
end
if f1>0,
 fprintf('Direct flops = %d and the FMM =%d\n',f1,f2) % Matlab Ver 5
end

fprintf('\tDifference in 2-norm = %e\n',  norm(P_d-P_f) )

function A=matrixq(x,y)  %% Note p_j=(x_j,y_j)
  n = length(x); A=[];
for i=1:n
for j=1:n
 dist = sqrt( (x(i)-x(j))^2 + (y(i)-y(j))^2) ;
 if dist < eps, dist=1; end
 A(i,j) = log(dist);
end
end 
