%% Agenda
% Multi threading
% Parallel pools
% cluster profiles (intro)
% Parallel Server cluster profiles
% parfor construct
% spmd construct
% distributed construct
% GPU computing
% OPTIONAL case study exercises
%
%%
%% Multithreading
%
% Matlab can execute a large number of operators and functions in parallel through multi treading
% Does anyone know what multi threading is and/or how Matlab does it?
%
% Examples of Matlab functions that are suitable for multi-threading.
% element wise operators (e.g. every dot function)
% standard operators (e.g. * , + ,  -. / )
% This is all done completely transparant to the user.
% By default Matlab will use as many threads as are available a node

% to check the number of computational thread
num_threads=feature('NumThreads')

% Is the answer what you expected?
%

% to set the number of computational threads
num_threads=feature('NumThreads', 2)

% Why is it important when working on ACES to control the number of threads?
%
%   Average Desktop has 4 to 8 cores. ACES has 96 cores per node (64 on Ice Lake Nodes)
%
%% Example: simple matrix multiplication
%
% Let's do a simmple matrix-matrix multiplication. Matlab will use an optimized multi threaded library function to do this
%

% create random matrix
A=rand(15000);
feature('NumThreads',1)
for count = 1:3
    tic;
    C=A*A;
    mytime = toc;
    fprintf('Computation time is %8.7f seconds.\n',mytime);
end

%%
%% EXERCISE: In the code above
% Vary the number  of threads, and see how run time changes
% change the size of the matrix (try with very small matrix dimensions )

%set num threads back to 8, will see later why
feature('NumThreads',8);

%
%%
%% PARALLEL POOLS
%
%% Matlab Workers and pools
%
% The Main concept in parallel processing is a Matlab worker. A worker is basically an indepedent Matlab process. A pool is a set of Matlab workers who can communicate with each other. The idea is to utilize the workers to split up the "problem" in smaller parts and  the workers work independently on their part. Hopefully speeding up the computation.
%%
%% Creating a parallel pool


% start a parallel pool
N=2;
mypool = parpool(N)

% How to check if they were actually started?
%
%
% For documentation about parpool and all the options:

%check documentation of parpool
doc parpool

%
%
%% Retrieving the parallel pool object
%
% Ideally, when you create the pool you assign it to a variable (like we did above), but if you didn't, no worry, Matlab provides an easy way to retrieve way to retrieve the pool object. Handy when the pool object was cleared. There is no real need to assign the parpool object to a variable

mypool = gcp('nocreate') % gcp --> Get Current Pool

%
%
%% Deleting a pool
%
% When you are done with the workers, don't forget to close the pool.  This will free up resources

delete(mypool) % or mypool.delete()

%
%% Let's take a look at the parallel preferences (might explain the 'nocreate' param above)
%
%
%%
%% CLUSTER PROFILES
%
% In the previous section we introduced parallel pools. Now we will get a bit deeper into this. So far we created a pool and all the workers were created on the same node. Matlab actually provides a lot of features to decide where and how workers should run. Some things a cluster profile does:
% How and where to start workers
% keep track of all "jobs"
% So far, we only used 'local' profile (named  Process in R2023a)
% all workers must run on the same node
% Part of the Matlab Paralel Toolbox
% limit up to to 512 in theory ( in practice limit is number of cores on the node)
%
%
%
% Let's see how to see all the available cluster profiles

%to see all available cluster profiles
parallel.clusterProfiles

%to see the default cluster profile
parallel.defaultClusterProfile()

%
% Now, let's take a look how to find all the cluster profiles in the GUI
%
%
% The next functions show how to get cluster objects, and how to change them

%to see all available cluster profiles
cluster_object = parcluster('local')

% What do we see?
%
% Let's change the number of workers and threads in the cluster object


% let's make some changes to this cluster profile
cluster_object.NumThreads=4;
cluster_object.NumWorkers=2;
cluster_object

% to save the cluster profile
cluster_object.saveAsProfile('mylocal');

%
%
% We can also change/check the default cluster profile

%to change the default cluster profile
parallel.defaultClusterProfile('mylocal');
parallel.defaultClusterProfile()

%let's justset default back to local
parallel.defaultClusterProfile('local');

% Let's check the cluster profiles in the GUI again?
%
%
%
%%  Parallel Server: Cluster Profiles
%
% workers can be distributed over multiple nodes
% number of workers limited only by MDCS license
% integrates with batch schedulers (SLURM on Grace)
%% Importing cluster profile
% Matlab comes standard with the 'local' profile. If you want to use addtional profiles you will have to create them or import them. HPRC already created a cluster profile for you, so you only need to import it

% to import a cluster profile using matlabbuilt-in function
% parallel.importProfile("location of profile"); 

% tamu convenience function
profile = tamuprofile.importProfile()

% Let's check the newly created cluster profile  in the GUI again
%

% What is the command (should say commands) to show contents of the profile
% above. 


%% EXERCISE: In the code above
% Investigate the AdditionalProperties field and find out what it contains. Does it look familiar.
%
%%
%% TAMUClusterProperties
%
% The TAMUClusterProperties class is a Matlab class developed by HPRC, to easily set requirement and properties (e.g. how many workers to use). The parallel cluster object will use this information when starting a parpool or batch job.
% NOTE, in the current Matlab version, everything can be done by directly manipulating the AdditionalProperties. In older versions that was not the case. This is a convenience funcitonality
%

% see the help page
help TAMUClusterProperties

%

% initialize a TAMUClusterProperties object
tp = TAMUClusterProperties()

%
%
% NOTE: some of the properties that can be set are needed to run on our cluster (e.g. wall time and memory)

% set requirements; 
tp.workers(8) % how many workers to start
tp.workers_per_node(8) % how many workers to start per node
tp.walltime('03:00')  % how long to use the workers

%
%
% Retrieve cluster object. Will use that later to start a parpool or offload work using batch

% to get a fully populated cluster object
% created by HPRC
profile = tamuprofile.parcluster(tp)

%
%
%%
%% Parpool revisited
% let's open a parpool again. This time using our cluster profile. We will see it's avery similar to starting a pool using the 'local' cluster profile

% open a parallel pool using the cluster profile
delete(gcp)
% is there anything that we need to do before creating pool
pool = parpool(profile) 

% Let's check where the parpool is running
%
%%
%% HPRC convenience functions
% HPRC developed convenience functions that use the TAMYClusterProperties object directly, so there is no need to get a cluster objet
%

% delete the current pool
delete(gcp);

%start the pool again using the convience function
tamuprofile.parpool(tp)

% Use the Slurm command to check the job
%
%%
%% Cool things you can do with cluster profiles
% We can do some cool things when we start workers on multiple nodes
% every worker can use multi threading for aditional speedup
% distribute the workers over different nodes (one worker per node) and each worker utilize a different GPU
% distribute the workers over different nodes (one worker per node) and evey worker can utilize all the memory on that node (useful for calculations that need much more memory than available on a single node)
% Anything else you can think of?
%%
% NOTE, In traditional clusters, every node has a fixed (( in most cases a limited) number of GPUs per node,so the only way to utilize more GPUs is to use multiple nodes. ACES is a composable cluster, so you can request a large number of GPUs and batch scheduler will give you access to those GPUs on a single node.
%%
%
%%
%% PARFOR
%
% In the previous section we started a parallel with Matlab workers. now we will actually start utilizing them. The first concept we will look at is parfor
% parfor is probably the most common and straightforward parallel construct in Matlab
%% Example 1

%feature('NumThreads',1);

% let's turn this regular for loop into a parfor loop
%pool=parpool(2);

tic
n = 800; A = 300; a = zeros(n);
for count = 1:n
    a(count) = max(abs(eig(rand(A))));
end
toc

%
%
%
%
%% You can adjust the number of workers in parfor
%

clearvars;
n = 800; A = 300; a = zeros(n);
tic;
% use the format below to specify the number of workers to use in parfor
parfor (count = 1:n, 2)
    a(count) = max(abs(eig(rand(A))));
end
toc

%% EXERCISE: In the code above
% Do a performance study. Change the code to vary the the number of workers and see how it affects the  runtime.
%
%% Limitations with parfor
% Loop index must be consecutive increasing integers ( e.g. i= N : 0 is not allowed)
% there is always some overhead in using parfor (both starting and iterating)
% Not every for loop can be changed into a parfor loop.

%classifying variables, what happens if we make this a parfor
N=1000;
a=rand(N);
for i = 2 : N
    a(i) = a(i) + a(i-1);
end
a

%% EXERCISE: In the code above
% Change the above code in a parfor and try to run it. What do you see?
%
% Checkout how Matlb classifies variables to determine if loop can run in parallel:
% https://www.mathworks.com/help/coder/ug/classification-of-variables-in-parfor-loops.html
%
%
%% DEMO: processing number of input files concurently
% (time permitted)
% Suppose we have a large number of input files (e.g. input1.mat, input2.mat ,....) that we want to process, normally will run them one by one, as follows:
% Read input file
% Run matlab on this input file
% Write results back to disk
% To process all, we just put a for loop around it and run it (or just run multiple matlab instances). How can we process the input files in parallel?
%


N=4

% create dummy input files for our demo
for i=1:N
    myinvar = ones(10,1) * i;
    save(['input' num2str(i) '.mat'],'myinvar');
end


% let's process the input matrices one by one
parfor count=1:N
    loadvar=load(['input' num2str(count) '.mat'])
    myvar= loadvar.myinvar
    myvar
    % normally this would be where all processing is done
    myoutvar = myvar*10;
    myoutvar
    % save output to file
  
    mysave(myoutvar,count)
    %save(['output',num2str(count),'.mat'],'myoutvar');
end
dummy=1

% create a parallel version of this below

% The above demo is just one example how to concurrently apply a function on multiple inpus. If the inputs were variables/values instead of files the approach would be similar (even a bit more straight forward).
%
%
%% SPMD
%
% The next parallel construct we will discuss is the SPMD construct. In an spmd block all the workers will execute the block concurrently. SPMD blocks are actually extensively used by other parallel constructs.
%

clearvars
N=2
spmd
    worker_id=spmdIndex; % in older Matlab versions it's called labindex
    num_workers=spmdSize; % in older Matlab versions it's called numlabs
    A=zeros(N);
    A=A+spmdIndex
end

% How to check the value of A and what do you think it is

%
%%
%% Data persistence
%
% One of the features of spmd blocks is data persistence. Even after the spmd block is finished, ALL variables created on the workers are still available

% worker_id, num_workers, and A are still accessible on the client
whos

worker_id

% to access individual elments use the {} operator
A{2} % will return the value of array A on worker 2

% can even change the value
A{2}=magic(N);

A{2}



%
% variables created inside a spmd block will stay on the workers after the spmd block finished.

spmd
    B=A*spmdSize
end
B

% what do you thing the value of B will be on the different workers?

%%
%
%% Communication between workers
%
% Matlab provided a number of functions for Matlab workers to send data to eachother
% labSend(var,id)         \to send variable var to worker with id id
% var = labReceive(id) \to receive variable from worker with id id and assign to var
% vf =  labSendReceive(it,if,vt) \to \;combine lab send and receive
%


spmd
    if (spmdIndex == spmdSize)
        sent=rand(4)*spmdIndex;
        % enter the command to send the value of "sent" here

    elseif (spmdIndex == 1)
        % enter the command to receive the sent value here
        received = 0; 

        fprintf('worker 1 received %f from worker %i',received,spmdSize);
    end
end

% what is the type of variable received?
%received

%what happens if we try to access received{2}
%received{2}

%% EXERCISE: In the code above
% Change the above code such that the worker with worker id sends a message and the worker with the lowest id ( that will be worker 1) receives the message and prints it
% What will be the type for variable "received" outside of the SPMD block
%
%
%
%% Distributed Arrays
%
%% Distributed array are conceptually the same as regular arrays/matrices. However, Under the hood the data is actually distributed among all the workers.
% most regular matrix operators will work on distributed array (overloaded functions)
% matrix elements can be of any type
%% How to create distributed arrays?
% create distributed array from regular array
% use the Matlab constructor functions (similar to ones(), zeros(),... )
% create directly on individual workers (ADVANCED, will not discuss here)

clearvars;
N=4;
M=8;

% create distributed array from regular array
A=ones(N,M) .* (1:M);
Ad=distributed(A);

whos

%

% create distributed array directly on workers
Bd=distributed.rand(N);

whos

%


% to check if a variable is distributed
a_distributed = isdistributed(A)
b_distributed = isdistributed(Ad)

classAd = class(Ad)

%
% To create a regular array from a distributed array:
%

Ac = gather(Ad);
classAc = class(Ac)
whos

%

% to see all available operators and methods
methods('distributed')

%
%%

% create the distributed arrays A and B

tic
%compute the multiplication. In this case let's just use Matlab operators

toc

%% EXERCISE: In the code above
% Add code to create two random matrices (let's say dimension 10000*10-00)
% Add code to perform the Matrix multiplication
%
%%
%% How is the array actually distributed among the workers?
%
% Although conceptually regular arrays, distrubuted arrays are actually distributed on the workers. We can access them directly inside spmd.

spmd
    % to see how an array is distributed
    local = getLocalPart(Ad)
    
    % to see global indices along the dimensions
    global_rows = globalIndices(Ad,1)
    global_columns = globalIndices(Ad,2)
end

%
%
%% GPU computing
%
% GPU computing is getting more and more popular and utilizing GPUs can speedup certain codes dramatically.
% What exactly is a GPU
% Accelerator card
% Thousends of simple computing cores
% Dedicated high speed memory
% How to check if a gpu is available:

% check how many GPU's are available on computer/node
number_gpus = gpuDeviceCount()
% list all the gpus
if (number_gpus > 0)
    gpuDevice()
end

%
%% Matlab provides

methods('gpuArray')

%
%% How to use the GPU for computation.
%
% Matlab will automatically perform a computations on the the GPU is the opexecute a function on the GPU if the operands involved in the computation are GPU variables (i.e. reside on the GPU)
%
% Since a GPU has its own memnory and Matlab only utilizes the GPU when the operands are located on the GPU, we need a way to transfer variables to the GPU. Using Matlab, you can do the following
% copy variable to the GPU
% create variable directly on GPU
% copy variable back to client

% copy variable from client workspace to the GPU
N=100;
A=rand(N);
Ag=gpuArray(A);

% create variable directly on the GPU
Bg = gpuArray.ones(N);

% copying variable back to client workspace
Ac=gather(Ag);

whos

% This is similar to how Matlab deals with distributed arrays
%
%%

% create the gpu arrays A and B

tic
%compute the multiplication. In this case let's just use Matlab operators

toc

%% EXERCISE: In the code above
% Add code to create two random matrices (let's say dimension 10000*10000)
% Add code to perform the Matrix multiplication
%
%
%% Back to the Slides
%
%%
%% CASE STUDY: COMPUTE PI USING MONTE CARLO METHOD
%
%% Assume a circle and an enclosing square box . If you pick N random points, approximately N*\pi/4 of these points will be inside the circle. Assuming the number of random points is large enough, the following formala can be used to approximate \pi
%% \pi =4*\frac{\sum_1^N \left({x\left(i\right)}^{2\;} +{y\left(i\right)}^2 \le 1\right)}{N}
%
%%
%% MONTE CARLO PI: PARFOR
%
%% A straight forward approach to code this in Matlab is:
% generate N randon cooordinates x  (all between 0\le x\le 1), do the same for y
% for every random point (x,y) check if it's  inside the circle using the formula \sqrt{x^2 +y^2 \;}\le 1
% \pi can be approximated using the formula \pi =4*\frac{\textrm{#points}\;\textrm{inside}}{N}
%

N=1000000000
% step 1: generate large number of random points with X and Y coordinate
xcoord=rand(N,1);
ycoord=rand(N,1);

%delete(gcp)
%lc = parcluster('local');
%lc.NumThreads=1;
%lc.NumWorkers=8
%pool=parpool(lc)
% step 2: iterate over all the random points and check if it's inside the loop
tic;
total_inside=0;
parfor count = 1 : N
    if ((xcoord(count)^2 + ycoord(count)^2 <= 1))
       total_inside = total_inside + 1;%(xcoord(count)^2 + ycoord(count)^2 <= 1);
    end
end
% step 3: compute pi
pi = 4*total_inside/N
mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);

%
%% MONTE CARLO PI: SPMD
%
% Every worker will create a subset of X and Y coordinates
% Ever worker will compute number of points inside
% add the partial sums ( either client aggreates the composite or workers send their results to each other, i.e. list communication)
% Client computes pi

N=1000000000;
tic;
clear NS xcoord ycoord tc;
spmd
    % worker creates suset of X and Y coordinates
    NS = N/spmdSize;
    xcoord=rand(NS,1);
    ycoord=rand(NS,1);
    tc=0;
    % compute inside points
    for i=1:NS
        if (xcoord(i)^2 + ycoord(i)^2 <= 1)
            tc=tc+1;
        end
    end
end
% 3 client combines partial results
global_inside=0;
for i=1:length(tc)
    global_inside = global_inside+tc{i};
end
%client computes pi
piresult = 4*global_inside/N;
toc
fprintf('The computed value of pi is %8.7f.\n',piresult);

%
%
%%
%% MONTE CARLO PI: DISTRIBUTED ARRAYS
%
% in the previous section we computed pi, both using a for loop and using matlab operators. Now we will do it using distributed arrays. We will do it both using matlab oprators and directly one the workers

% directly
clearvars
N=8000000;
inside=0;
% step 1: generate large number of random points with X and Y coordinate

% step 2: Use the element wise operators to compute sum of points in circle
pi = 4*inside/N;

tic;

mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime)

%
%%
%
%% MONTE CARLO PI: GPU
%
% We will compute pi using the monte carlo method. This time we will compute it on the gpu and compare results with the cpu versions.


%PI
N=320000000;
% step 1: generate large number of random points with X and Y coordinate


% step 2: Use the element wise operators to compute
tic;


%mytime=toc;
%fprintf('Computation time is %8.7f seconds\n',mytime);

% WHAT HAPPENS TO TIME WHEN OUTPUT IS NOT SURPRESSED

% use gputimeit to get precise time measurement
f = @() pi_gpu(N);
gtime = gputimeit(f)

%
%
%
%%
%% Let's run the matix mult example from the multithreading section and compare run time

% the cpu version
% 
A=rand(16000);

for count = 1:3
    tic;
    C=A*A;
    mytime = toc;
    fprintf('Computation time is %8.7f seconds.\n',mytime);
end

%
%
%
%
%
%%
%



%
%% MONTE CARLO PI: MULTITHREADING
%
% Now let;sclook at multi threading. We saw before many matlab operators and build-in functions can run multi=threaded. So far we have used a for loop, which is not suitable for multi-threading. Let's redo the monte carlo method again, instead of using a for use only matlab operators.

X=rand(N,1);
Y=rand(N,1);
tic;

mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);

%
% Now let's run the code and vary the number of threads.
%
%%
%% MONTE CARLO PI: GPU
%
% We will compute pi using the monte carlo method. This time we will compute it on the gpu and compare results with the cpu versions.


%PI
N=320000000;
% step 1: generate large number of random points with X and Y coordinate


% step 2: Use the element wise operators to compute
tic;


%mytime=toc;
%fprintf('Computation time is %8.7f seconds\n',mytime);

% WHAT HAPPENS TO TIME WHEN OUTPUT IS NOT SURPRESSED

% use gputimeit to get precise time measurement
f = @() pi_gpu(N);
gtime = gputimeit(f)

%%
%
%%
%% THAT'S ALL FOLKS, BACK TO THE SLIDES
%
%
%
%
% THIS SECTION CONTAINS ALL THE FUNCTIONS CALLED IN MATLAB SNIPPETS ABOVE

function [pi_val ,mytime] = pi_fun(N , numthreads)
feature('NumThreads',numthreads);
X=rand(N,1);
Y=rand(N,1);
tic;
pi_val = 4 * sum(X.*X + Y.*Y <= 1)/N;
mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);
pause(20)
end


function res= pi_gpu(N)
% step 1: generate large number of random points with X and Y coordinate
X=gpuArray.rand(N,1);
Y=gpuArray.rand(N,1);

% step 2: Use the element wise operators to compute
res = 4 * sum(X.*X + Y.*Y <= 1)/N;
end

function mysave(output,count)
save(['output',num2str(count),'.mat'],'output');
end
