%% Numerical integration
% This example shows the basics of working with spmd statements, and how they
% provide an interactive means of performing parallel computations.  We do
% this by decomposing the domain for parallel numerical integration. 
%
% Related Documentation:
% 
% * <matlab:doc('spmd') spmd reference page> in the 
% Parallel Computing Toolbox(TM) User's Guide
%
% Related Examples:
%
% * <docid:distcomp_examples.example-ex58432651 Using GOP to Achieve MPI_Allreduce Functionality>
%   Copyright 2007-2013 The MathWorks, Inc.
%  Modified by HA
%% Introduction
% Integrate "labvise" in parallel.
%
% $$ \int_0^1 {4 \over {1 + x^2}} dx$$
%
% We intend to have the parallel pool perform the calculations in
% parallel, and to use the |spmd| keyword to mark the parallel blocks of
% code.  We first look at the size of the parallel pool that is currently
% open. 
p = gcp;
numlabs=p.NumWorkers
%% 
% We divide the work between the workers (labs) by having each worker calculate 
% the integral of the function over a subinterval of [0, 1] as shown in the
% picture. 
%
% <<../paralleldemo_quadpi_mpi_01.png>>

%%
% We define the variables |a| and |b| on all the workers, but let their values
% depend on |labindex| so that the intervals [a, b] correspond to the
% subintervals shown in the figure. We then verify that the intervals are
% correct.  Note that the code in the body of the spmd statement is
% executed in parallel on all the workers in the parallel pool.
spmd
    a = (labindex - 1)/numlabs;
    b = labindex/numlabs;
    fprintf('Subinterval: [%-4g, %-4g]\n', a, b);
end

%% 
% We let all the workers now use a MATLAB quadrature method to approximate
% each integral.  They all operate on the same function, but on the
% different subintervals of [0,1] shown in the figure above.
f=@(x) 4./(1+x.^2)
%%
spmd
    LabIntegral = integral(f, a, b);
    fprintf('Subinterval: [%-4g, %-4g]   LabIntegral: %4g\n', ...
            a, b, LabIntegral);
end
%% Add the Results
% The workers have all calculated their portions of the integral of the
% function, and we add the results together to form the entire integral
% over [0, 1].  We use the |gplus| function to add |LabIntegral| across all
% the workers and return the sum on all the workers.
spmd
    sumIntegrals = gplus(LabIntegral,1);
end
%% Inspect Results in the Client  
% Since the variable |piApprox| was assigned to inside an spmd statement,
% it is accessible on the client as a Composite.  Composite objects
% resemble cell arrays with one element for each worker.  Indexing into a
% Composite brings back the corresponding value from the worker to the client.
format long
Int_on_client=sumIntegrals{1}
%%
delete(gcp)
%%
numlabs=24;
parpool(numlabs)
%%
 clear
%%
numlabs=24;
f=@(x) 4./(1+x.^2)
%% The whole job without extra displays.
spmd
    a = (labindex - 1)/numlabs;
    b = labindex/numlabs;
    LabIntegral = integral(f, a, b);
    sumIntegrals = gplus(LabIntegral,1);        
end
Int_on_Client=sumIntegrals{1}

%% Bring result to client 
% Has to be done outside the spmd-block, before pool is closed.
sumIntegrals{1}
%% Which variables are left on client after:
delete(gcp)
whos
%{
>> a 
a = 
    Invalid Composite (the parallel pool in use has been closed).
%}
%% Develop: Use a lower degree method on small subintervals
%  Sum up to credit from parallel computing.
% trapzd, better still: write the code: simpson.
%% Or: Function with fast changing parts on a long interval.
% Especially piecewise defined function with (almost) singularities.

