INTRODUCTION TO MATLAB PARALLEL COMPUTING TOOLBOX Kadin Tseng Boston University Scientific Computing and Visualization MATLAB Parallel Computing Toolbox What is the PCT ? The Parallel Computing Toolbox is a MATLAB tool box. This tool box provides parallel utility functions to enable users to run MATLAB operations or procedures in parallel to speed up processing time. 2 MATLAB Parallel Computing Toolbox 3 Where To Run The PCT ? Run on a desktop or laptop MATLAB must be installed on local machine Starting with R2011b, up to 12 processors can be used;

up to 8 processors for older versions Must have multi-core to gain speedup. The thin client you are using has a dual-core processor Requires BU userid (to access MATLAB, PCT licenses) Run on a Katana Cluster node (as a multi-cored desktop) Requires SCV userid Run on multiple Katana nodes ( for up to 32 processors) Requires SCV userid MATLAB Parallel Computing Toolbox Types of Parallel Jobs ? There are two types of parallel applications. Distributed Jobs task parallel Parallel Jobs data parallel 4 MATLAB Parallel Computing Toolbox 5 Distributed Jobs This type of parallel processing is classified as: Multiple tasks running independently on multiple workers with

no information passed among them. On Katana, a distributed job is a series of single-processor batch jobs. This is also known as task-parallel, or embarrassingly parallel , jobs. Examples of distributed jobs: Monte Carlo simulations, image processing Parallel utility function: dfeval MATLAB Parallel Computing Toolbox 6 Parallel Jobs A parallel job is: A single task running concurrently on multiple workers that may communicate with each other. On Katana, this results in one batch job with multiple processors. This is also known as a data-parallel job. Examples of a parallel job include many linear algebra applications : matrix multiply; linear algebraic system of equations solvers; Eigen solvers. Some may run efficiently in parallel and others may not. It depends on the underlining algorithms and operations. This also include jobs that mix serial and parallel processing. Parallel utility functions: spmd, drange, parfor, . . . MATLAB Parallel Computing Toolbox

7 How much work to parallelize my code ? Amount of effort depends on code and parallel paradigm used. Many MATLAB functions you are familiar with are overloaded to handle parallel operations based on the variables data type. MATLAB Parallel Computing Toolbox 8 Distributed Jobs dfeval Example: Run a dfeval job interactively Computes 1x4, 3x2 random arrays on 2 cores, locally or on Katana application mrows or SGE columns Parallel config. file >> y = dfeval(@rand, {1 3}, {4 2}, 'Configuration',local) run on 2 cores Submitting task 1 Job output will be written to: /usr1/scv/kadin/Job6_Task1.out QSUB output: Your job 211908 ("Job6.1") has been submitted Submitting task 2

Job output will be written to: /usr1/scv/kadin/Job6_Task2.out QSUB output: Your job 211909 ("Job6.2") has been submitted y = [1x4 double] [3x2 double] Job ran in batch or background, output returns to client workspace. MATLAB Parallel Computing Toolbox 9 Distributed Jobs SCV script For task-parallel applications on the Katana Cluster, we strongly recommend the use of an SCV script instead of dfeval. This script does not use the PCT. For details, see http://www.bu.edu/tech/research/computation/linux-cluster/katana-clu ster/runningjobs/multiple_matlab_tasks/ If your application fits the description of a distributed job, you dont need to know any more beyond this point . . . MATLAB Parallel Computing Toolbox 10 How Do I Run Parallel Jobs ? Two ways to run parallel jobs: pmode or matlabpool.

Use these procedure turns on (or off) parallelism and allocate (or deallocate) resources. pmode is a special mode of application; useful for learning the PCT and parallel program prototyping interactively. matlabpool is the general mode of application; it can be used for interactive and batch processing. MATLAB Parallel Computing Toolbox 11 pmode >> pmode start local % assuming 4 workers are available Above is a MATLAB window on a Katana node. A separate Parallel Command Window is spawned. (Similar for Windows) In Room B27, each thin client has 2 cores (workers). MATLAB Parallel Computing Toolbox 12 pmode

>> pmode start local % use 4 workers in SPMD (parallel) mode PCT terminologies: worker = processor labindex: processor number numlabs: Number of processors Worker number Command entered Enter command here Any command issued at the P>> prompt is executed on all workers. Enter labindex to query for workers ID. Use if conditional with labindex to issue instructions to specific workers, like this: if labindex==1, numlabs, end pmode Replicate array P>> A = magic(3); % A is replicated on every worker

Variant array P>> A = magic(3) + labindex 1; LAB 1 | |8 |3 |4 1 5 9 LAB 2 | 6|9 7|4 2|5 2 6 10 Private array P>> if labindex==2, LAB 1

| | | |9 |undefined|4 | |5 Spring 2012 | 7|10 8| 5 3| 6 3 7 11 LAB 4 | 8|11 9| 6 4| 7 4 8

12 9 10 5 A = magic(3) + labindex 1; end LAB 2 2 6 10 LAB 3 % labindex=1,2,3,4 LAB 3 LAB 4 | | 7| | 8|undefined|undefined

3| | 13 Switch from pmode to matlabpool If you are running MATLAB pmode, exit it. P>> exit MATLAB allows only one parallel environment at a time. pmode may be started with the keyword open or start. matlabpool can only be started with the keyword open. You can also close pmode from the MATLAB window: >> pmode close Now, open matlabpool >> matlabpool open local Spring 2012 % rely on default worker size 14 matlabpool matlabpool is the essential tool for requesting processors for MATLAB parallel computing. Usage: >> matlabpool open Config Nworkers % italic optional item >> % Use default value if optional item is omitted

>> % Config is local or SGE on Katana >> % Nworkers is no. of cores; if omitted, default to 2 in Rm B27 >> % >> % perform parallel tasks inside matlabpool region >> % . . . >> matlabpool close % ends matlabpool Parallel method choices within matlabpool region: parfor parallel for-loop; cant mix with spmd spmd single program multiple data parallel region Arrays are distributed among the cores with codistributed and distributed commands drange parallel for-loop within spmd region; parfor-like Spring 2012 15 Parallel Methods parfor It is a for-loop that executes on multiple workers in parallel parfor starts and ends on Client Number of workers is determined by matlabpool Without matlabpool activation, parfor reverts to for Computation for all loop indices must be independent Results must be independent of loop index execution order

parfor is responsible for work load distribution via loop indexing More rules governing parfor operations in the PCT doc. Spring 2012 16 Parallel Methods parfor (2) matlabpool open % by default, say, 2 workers s = 0; parfor i=1:10 x(i) = sin(2*pi*i/10); % i=1:5 by 1st worker; i=6:10 by 2nd s = s + i; % not all reduction operations allowed end matlabpool close parfor can do reductions provided that the operator satisfies associative rule (e.g., +, *). Example. x + (y + z) = (x + y) + z Subtract () and divide (/) operators fail rule -- indeterministic s=1000; parfor i=1:500, s=s*i, end, s % do a few times s=1000; parfor i=1:500, s=s/i, end, s % do a few times % Example of loop dependency

% matlab editor will tell you this parfor usage is wrong a = zeros(1,100); parfor i = 2:100 Springa(i) 2012 = myfct(a(i-1)); end 17 Integration Example An integration of the cosine function between 0 and / 2 Integral sum of areas of rectangles (height x width) Several parallel methods will be demonstrated. p b n aij h cos( x)dx a

i 1 j 1 Worker 1 aij p cos( x )dx cos(x) h Worker 3 Worker 4 Spring 2012 i 1 mid-point of increment Worker 2 x=a

n h cos( aij 2 )h j 1 a = 0; b = pi/2; % range m = 8; % # of increments h = (b-a)/m; % increment p = numlabs; n = m/p; % inc. / worker ai = a + (i-1)*n*h; aij = ai + (j-1)*h; x=b 18 Integration Example Serial Integration % serial integration (with for-loop) tic m = 10000; a = 0; % lower limit of integration

b = pi/2; % upper limit of integration dx = (b a)/m; % increment length intSerial = 0; % initialize intSerial for i=1:m x = a+(i-0.5)*dx; % mid-point of increment i intSerial = intSerial + cos(x)*dx; end toc dx X(1) = a + dx/2 Spring 2012 X(m) = b - dx/2 19 Integration Example parfor This example performs parallel integration with parfor. matlabpool open 4 tic

m = 10000; a = 0; b = pi/2; dx = (b a)/m; % increment length intParfor = 0; parfor i=1:m intParfor = intParfor + cos(a+(i-0.5)*dx)*dx; end toc matlabpool close Spring 2012 20 Parallel Methods spmd In spmd parallel environment: A MATLAB client process manages multiple workers Each worker, or lab, has its own independent memory Each worker runs the same program (Single Program) but on its own data (hence Multiple Data)

Can perform point-to-point communications between workers (labsend and labreceive) Collective communications among workers (gather, gop) Workers can access Client data by simply referencing it Client can create/modify/access data on workers Two very essential functions within spmd region: numlabs returns number of labs available for job labindex returns current lab ID; this is essential to control what to do on specific labs >> spmd, labindex, end % prints local worker ID (1:numlabs) spmd region is enclosed between spmd & end Spring 21 following 2012 Arrays created in one region are preserved for Integration Example spmd n = 500; % number of integral increments per lab a = 0; b = pi/2; % lower and upper limits of integration spmd deltax = (b - a)/numlabs; % integral range length per lab ai = a + (labindex - 1)*deltax;

% local integral limits, [ai, bi] bi = a + labindex*deltax; dx = deltax/n; % increment length for each rectangle x = ai+dx/2:dx:bi-dx/2; % mid-points of n increments per worker local_int = sum(cos(x)*dx); % integral sum for current worker intSPMD = gplus(local_int,1); % total sum saved in lab 1 only end % spmd local_int{2} % prints local_int of lab 2 (from Client) L = [local_int{:}]; % saves local_int from labs to L on Client; needs [ ] Variables defined in spmd are class composite - check with whos Composite array intSPMD is NULL on labs 2:end; set by gplus n=2; spmd, a=rand(n), end, A=[a{:}]; % local a is 2D; A on Client is 2x8 B = a{1}; for j=2:4, B(:,:,j) = a{j}; end % B is 2x2x4 3D double array

Spring 2012 22 Integration Example drange Similar to parfor but used within spmd n = 500; % number of integral increments per lab a = 0; b = pi/2; % lower and upper limits of integration spmd deltax = (b - a)/numlabs; % increment length per worker for i=drange(1:numlabs) ai = a + (i - 1)*deltax; % local integral limits, [ai, bi] bi = a + i*deltax; dx = deltax/n; % increment length per ractangle x = ai+dx/2:dx:bi-dx/2; % mid-points of n increments local_int = sum(cos(x)*dx); end intDrange = gplus(local_int, 1); % send total sum to lab 1 end % spmd myint = intDrange{1}; Spring 2012

% collect integral from lab 1 23 Integration Example drange Similar to parfor but used within spmd m = 10000; % number of integral increments for all labs a = 0; b = pi/2; % lower and upper limits of integration spmd dx = (b - a)/m; % increment length per rectangle local_int = 0; for i=drange(1:m) local_int = local_int + cos(a+(i-0.5)*dx)*dx; end intDrange = gplus(local_int, 1); % send total sum to lab 1 end % spmd myint = intDrange{1}; % send integral to client Advantage: let drange do work load distribution Disadvantage: code not vectorized

Spring 2012 24 Array Distributions Objective is to distribute data, hence workload, to workers for shorter walkclock time. Data distribution utilities: distributed Data distributed and accessible readily from Client. Data distributed along the last dimension codistributed Run within spmd to distribute data on labs along arbitrary dimension composite - Non-distributed variables and arrays created within spmd belongs to this class Spring 2012 25 How To Distribute Arrays Ways to distribute data: Partitioning a larger array

Often, this larger array already exist on Client. Distributing it to labs will take up further memory as well as run time. Building from smaller arrays Array sections are constructed directly on individual labs resulting in smaller memory footprint. Created with MATLAB constructor function (rand, zeros, . . .) Supported by distributed and codistributed Example. >> d = distributed.rand(n) % distributes nxn rand Spring 2012 26 matlabpool open 'local' 4 Array Distribution A = magic(8); % A is on Client; Usage accessibleExamples from labs ad = distributed(A); % on Client, divide A to 8x2s along last dim N = size(A,2)/matlabpool(size); % columns distributed per lab spmd

n = size(A,2)/numlabs; % N is double; n is composite al = getLocalPart(ad); % al is local to the current lab al(:,1) = 0; % modify locally on lab; use local index ad(:,1+(labindex-1)*N) = 0; % not local; must use global index aa = gcat(al); % concatenate variant array al to replicated aa ab = codistributed(aa); % co-distribute along last dim (default) br = codistributed(A,codistributor1d(1)); % distribute by rows % partition from larger array b = codistributed.rand(8); % built from smaller arrays end Spring 2012 B = [al{:}]; % combine labs 8x2 arrays to 8x827array on Client Array Distribution Usage Examples (2) >> whos Name A B N aa ab

ad al b br n Spring 2012 Size Bytes 8x8 8x8 1x1 1x8 8x8 8x8 1x8 8x8 8x8 1x8 512 512 8

2041 2077 2077 2041 2077 2077 2041 Class Attributes double double double Composite distributed distributed Composite distributed distributed Composite 28 Data Parallel Example Matrix Multiply

Computations of distributed arrays need not be in an spmd region. The operations, however, are governed by spmd paradigm. >> >> >> >> >> >> >> >> >> matlabpool open 4 n = 3000; A = rand(n); C = A * B; % maxNumCompThreads(1);% C1 = A * B; % a = distributed(A); % b = distributed(B); % c = a * b; % matlabpool close B = rand(n);

run with 4 threads set threads to 1 run on single thread distributes A, B along columns on client a, b on workers; accessible from client run on workers; c is distributed Wall clock time, in seconds, for the above operations C1 = A * B C=A*B a = distribute(A) c=a*b (1 thread) (4 threads) b = distribute(B) (4 workers) 12.06 3.25 2.15 3.89 The cost for distributing the matrices is recorded separately as

this is incurred only once over the life of the distributed matrices. Time required to distribute matrices is not significantly affected by matrix size. Spring 2012 29 Additional Ways to Distribute Matrices More efficient to distribute arrays within spmd: All workers participate in distribution (co-distribute) No need to decompose array on client, then copy to individual workers as with distributed. >> matlabpool open 4 >> A = rand(3000); B = rand(3000); >> spmd p = rand(n, codistributor1d(1)); % 2 ways to directly create q = codistributed.rand(n); % distributes rand array by col. s = p * q; % run on workers; s is codistributed % distribute matrix after it is created u = codistributed(A, codistributor1d(1)); % by row v = codistributed(B, codistributor1d(2)); % by column w = u * v; % run on workers; w is distributed end >> matlabpool close

Spring 2012 30 A Case Of parfor vs spmd % First, consider the serial code N = 40; B = rand(N); C = rand(N); p = ones(N,1)*0.5; v = p; a = rand(N,1); for i = 1:1000 % iterative loop for j=1:N vn(j)=foo(v,p(j),a,B(j,:),C(j,:)); % compute new v end if max(abs(vn-v))< 1e-7 % check for convergence fprintf('****** Converged at iteration %5d\n',i); break; end v = vn; % update v end Spring 2012

31 The parfor version matlabpool open 4 N = 40; B = rand(N); C = rand(N); p = ones(N,1)*0.5; v = p; a = rand(N,1); for i = 1:1000 % iterative loop parfor j=1:N vn(j)=foo(v,p(j),a,B(j,:),C(j,:)); % compute new v end if max(abs(vn-v))< 1e-7 % check for convergence fprintf('****** Converged at iteration %5d\n',i); break; end v = vn; % update v end matlabpool close Parfor and other parallel paradigms have high overhead cost to start. Avoid using parfor within an iterative loop. Spring 2012

32 spmd version matlabpool open 4 N = 40; spmd % Single Program Multiple Data parallel paradigm B=codistributed.rand(N,codistributor('1d',1)); C=codistributed.rand(N,codistributor('1d',1)); p=codistributed.ones(N,1,codistributor())*0.5; v = p; % initial guess; v is distributed by virtue of p for i = 1:1000 % iterative loop vn=foo(v,p,a,B,C); % compute new v if max(abs(vn-v))< 1e-7 % check for convergence printf('****** Converged at iteration %5d\n',i); break; end v = vn; % update v end end % spmd matlabpool close

SPMD need be started once. Distributed arrays remain in the worker space until spmd region ends. Spring 2012 33 Preferred Way to Distribute Matrices ? For matrix-matrix multiply, there are 4 combinations on how to distribute the 2 matrices (by row or column). While all 4 ways lead to the correct solution, some perform better than others. n = 3000; A = rand(n); B = rand(n); spmd ar = codistributed(A, codistributor1d(1)) % distributed by row ac = codistributed(A, codistributor1d(2)) % distributed by column br = codistributed(B, codistributor1d(1)) % distributed by row bc = codistributed(B, codistributor1d(2)) % distributed by column crr = ar * br; crc = ar * bc; ccr = ac * br; ccc = ac * bc; end Wall clock times of the four ways to distribute A and B C (row x row) C (row x col) C (col x row) C (col x col) 2.44

2.22 3.95 3.67 The above is true for on-node or across-nodes runs. Across-nodes array distribution and runs are more likely to be slower than on-node ones due to slower communications. Specifically for MATLAB, Katana uses 100-Mbits Ethernet for across-nodes communications instead of Gigabit InfiniBand. Spring 2012 34 Linear algebraic system Example: Ax = b matlabpool open 4 % serial operations n = 3000; M = rand(n); x = ones(n,1); [A, b] = linearSystem(M, x); u = A\b; % solves Au = b; u should equal x clear A b % parallel operations in spmd spmd m = codistributed(M, codistributor('1d',2)); % by column y = codistributed(x, codistributor(1d,1)); % by row

[A, b] = linearSystem(m, y); v = A\b; end clear A b m y % parallel operations from client m = distributed(M); y = distributed(x); [A, b] = linearSystem(m, y); W = A\b; matlabpool close function [A, b] = linearSystem(M, x) % Returns A and b of linear system Ax = b A = M + M'; % A is real and symmetric b = A * x; % b is the RHS of linear system Spring 2012 35 Task Parallel vs. Data Parallel matlabpool open 4 n = 3000; M = rand(n); x = ones(n,1); % Solves 4 cases of Ax=b sequentially, each with 4 workers for i=1:4 m = distributed(M); y = distributed(x*i); [A, b] = linearSystem(m, y); % computes with 4 workers u = A\b;

% solves each case with 4 workers end clear A b m y % solves 4 cases of Ax=b concurrently (with parfor) parfor i=1:4 [A, b] = linearSystem(M, x*i); % computes on 1 worker v = A\b; % solves with 1 worker end % solves 4 cases of Ax=b concurrently (with drange) spmd for i=drange(1:4) [A, b] = linearSystem(M, x*i); % computes on 1 worker w = A\b; % 1 worker function [A, b] = linearSystem(M, x) end % Returns A and b of linear system Ax = b end A = M + M'; % A is real and symmetric matlabpool close b = A * x; % b is the RHS of linear system Spring 2012 36

7. How Do I Parallelize My Code ? 1. Profile serial code with profile. 2. Identify section of code or function within code using the most CPU time. 3. Look for ways to improve code section performance. 4. See if section is parallelizable and worth parallelizing. 5. If warranted, research and choose a suitable parallel algorithm and parallel paradigm for code section. 6. Parallelize code section with chosen parallel algorithm. 7. To improve performance further, work on the next most CPU-time-intensive section. Repeats steps 2 7. 8. Analyze the performance efficiency to know what the sweet-spot is, i.e., given the implementation and platforms on which code is intended, what is the minimum number of workers for speediest turn-around (see the Amdahls Law page). Spring 2012 37 8. How Well Does PCT Scales ? Task parallel applications generally scale linearly. Data parallel applications parallel efficiency depend on individual code and algorithm used. My personal experience is that the runtime of a

reasonably well tuned MATLAB program running on single or multiprocessor is at least an order of magnitude slower than an equivalent C/FORTRAN code. Additionally, the PCTs communication is based on the Ethernet. MPI on Katana uses Infiniband and is faster than Ethernet. This further disadvantaged PCT if your code is communication-bound. Spring 2012 38 Speedup Ratio and Parallel Efficiency S is ratio of T1 over TN , elapsed times of 1 and N workers. f is fraction of T1 due to code sections not parallelizable. S T1 TN T1 1 as N

1 f f (f )T1 N Amdahls Law above states that a code with its parallelizable component comprising 90% of total computation time can at best achieve a 10X speedup with lots of workers. A code that is 50% parallelizable speeds up two-fold with lots of workers. The parallel efficiency is E = S / N Program that scales linearly (S = N) has parallel efficiency 1. A task-parallel program is usually more efficient than a dataparallel program. Parallel codes can sometimes achieve super-linear behavior due to efficient cache usage per worker. Spring 2012 39 Example of Speedup Ratio & Parallel Efficiency Spring 2012 40

Batch Processing On Katana MATLAB provides various ways to submit and run jobs in the background or in batch. SCV recommends a simple and effective method for all PCT batch processing on Katana. Cut-and-paste the below into a file, say, mbatch #!/bin/csh -f # MATLAB script for running serial or parallel background jobs # For MATLAB applications that contain MATLAB Parallel Computing Toolbox # parallel operations request (matlabpool or dfeval), a batch job # will be queued and run in batch (using the SGE configuration) # Script name: mbatch (you can change it) # Usage: katana% mbatch # : name of m-file to be executed, DONOT include .m ($1) # : output file name; may include path ($2) nohup matlab nodisplay -r $1 exit >! $2 & katana% chmod +x mbatch give mbatch execute attribute Katana% mbatch my_mfile myOutput Do not start MATLAB with -nojvm. PCT requires Java. Spring 2012

41 Use local Config. on Katana You can use the Katana Cluster as if it is local with 4 workers for interactive usages: 1. Request 4 processors on the same node (You need x-win32 on your client computer) Katana% qsh pe omp 4 2. In the new X-window, run matlab Katana% matlab Request a node with 4 processors. 3. Request workers in MATLAB window >> matlabpool open local 4 or >> pmode start local 4 . Spring 2012 Only need a PCT

license; no worker licenses needed. 42 Communications Among Workers, Client Communications between workers and MATLAB client pmode: lab2client, client2lab matlabpool: A = a{1} % from lab 1 to client Collective communications among workers gather, gop, gplus, gcat Spring 2012 43 MPI Point-to-Point Communications MPI point-to-point communication among workers labSend and labReceive % Example: each lab sends its lab # to lab 1 and sum data on lab 1 matlabpool open 4 spmd % MPI requires spmd a = labindex; % define a on workers If labindex == 1 % on worker 1 . . . % lab 1 is designated to accumulate sum on a

s = a; % sum s starts with a on worker 1 for k = 2:numlabs % loop over remaining workers s = s + labReceive(k); % receive a from worker k, then add to s end else % for all other workers . . . labSend(a,1); % send a on workers 2 to 4 to worker 1 end end % spmd indexSum = s{1}; % copy s on lab 1 to client matlabpool close It is illegal to communicate with itself. Spring 2012 44 MATLAB Parallel Computing Toolbox 45 Useful SCV Info SCV home page (www.bu.edu/tech/about/research/scv)

Resource Applications www.bu.edu/tech/accounts/special/research/ Help System [email protected], bu.service-now.com Web-based tutorials ( www.bu.edu/tech/about/research/training/) (MPI, OpenMP, MATLAB, IDL, Graphics tools) HPC consultations by appointment Kadin Tseng ([email protected]) Yann Tambouret ([email protected])