Parallel Computing Toolbox™ 3.3

Parallel Programming in MATLAB®

Distributed Arrays and Parallel MATLAB Functions

Distributed arrays are special arrays that store segments of data on participating labs (sets of MATLAB® workers that work together as peers in a parallel computation). You can perform operations such as transforms and decompositions directly on distributed arrays, just as you would with regular arrays.

Because distributed arrays store data on many labs, you can handle larger data sets than you can on a single MATLAB session. Support for more than 150 MATLAB functions enables you to quickly develop and prototype your parallel applications without having to manage low-level details of message passing. Some of the linear algebra functions are based on ScaLAPACK library. Refer to Using MATLAB Functions on Distributed Arrays for a full list of supported functions.

You can construct distributed arrays in several ways:

• Directly using constructor functions such as rand, ones, zeros, and true
• Concatenating variant arrays, which are arrays with the same name that contain different data on each lab
• Distributing replicated arrays, which are arrays with the same name and same data on all labs
• You can control how the data contained in these arrays is distributed over different labs.

% Create a distributed array with default distribution scheme. 
Number of workers is not assumed
dist = distributor();   
n = 1000;
on = ones(n-1, 1, dist);% A distributed vector of ones
 
% Construct a distributed tridiagonal matrix
DM = 2 * eye (n, dist); % EYE(..., dist) creates a distributed identity matrix 
DL = diag(on, -1  );    % A distributed matrix with lower diagonal set to ones
DU = diag(on,  1  );    % A distributed matrix with upper diagonal set to ones
 
A = DM + DL + DU;       % A is a distributed tridiagonal matrix
 
% Compute eigenvalue decomposition. 
[V, D] = eig(A);        % EIG works directly on distributed arrays
 
% Verify that the decomposition is valid. 
norm(A*V - V*D, 1);     % NORM, matrix multiply work directly on distributed 
		        % arrays
 
% Verify that the eigenvectors are orthonormal. 
norm(V'*V - eye(n, dist), 1);
 
% Find the maximum eigenvalue. 
maxD = max(diag(D));    % MAX and DIAG work directly on distributed arrays

Programming with distributed arrays. Distributed arrays and associated parallel algorithms let you create parallel applications with minimal changes to your serial code.

In general, very few changes are required to convert your serial program to a parallel program. Note that in the NAS Conjugate Gradient Benchmark example below, the only requirement for changing between serial and parallel code is the declaration of distributed arrays. Other operations work transparently on distributed arrays.

Refer to the Working with Distributed Arrays section in the Parallel Computing Toolbox documentation for more information.

function CG
%CG     NAS Conjugate Gradient benchmark
%   This function is similar to the NAS CG benchmark described in:
%   http://www.nas.nasa.gov/News/Techreports/1994/PDF/RNR-94-007.pdf
%   See code on page 19-20 for the pseudo code translated to m below.
 
n = 1400; nonzer = 7; lambda = 20; niter = 15;
nz = n * (nonzer + 1) * (nonzer + 1) + n * (nonzer + 2);
 
% Create A, a distributed sparse random matrix
% using default distribution scheme returned by distributor() 
A = sprand(n, n, 0.5 * nz/n^2, distributor());
A = 0.5 * (A + A');
 
% I is a distributed sparse identity matrix
I = speye(n, distributor()); 
A = A - lambda * I;             % Shift at lambda, an approx eigenvalue
 
% x is a vector of ones
x = ones(n, 1);  
%% Outer iteration: exactly same code will work for distributed and
%% non-distributed arrays. 
for iter = 1 : niter
   [z, rnorm] = cgit(A, x);     % Conjugate Gradient iteration. 
   zeta = lambda + 1/(x' * z);
   x = z / norm(z);
 
   if labindex == 1
       fprintf('%5d %13.4e %20.13f\n', iter, rnorm, zeta);
   end
end
 
%% ------------------------------------------------------------ %
function [z, rnorm] = cgit(A, x)
%CGIT  Conjugate Gradient iteration
%   [z, rnorm] = cgit(A, x) generates approximate solution to A*z = x.
z = zeros(size(x));
r = x;
rho = r' * r;
p = r;
 
for i = 1 : 15
   q = A * p;
   alpha = rho / (p' * q);
   z = z + alpha * p;
   rho0 = rho;
   r = r - alpha * q;
   rho = r' * r;
   beta = rho/rho0;
   p = r + beta * p;
end

rnorm = norm(x - A * z);

Working with distributed arrays: Implementation of the NAS Conjugate Gradient Benchmark. The declaration of distributed arrays, highlighted, is the only change required to convert the serial program to a parallel program.