Skip to content

Commit 3b9fb87

Browse files
committed
Adding MATLAB post source code
1 parent 860e2a2 commit 3b9fb87

10 files changed

Lines changed: 678 additions & 0 deletions

MATLAB_arrayfun/ArrayfunArticle.m

Lines changed: 474 additions & 0 deletions
Large diffs are not rendered by default.

MATLAB_arrayfun/README.md

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
Companion Code for "High-Performance MATLAB with GPU Acceleration" by Joss Knight
2+
============================
3+
4+
This folder contains source code from the NVIDIA Parallel Forall Blog post on MATLAB by Joss Stone (The Mathworks) entitled [High-Performance MATLAB with GPU Acceleration](http://developer.nvidia.com/parallel-forall/high-performance-matlab-gpu-acceleration).
5+
6+
License
7+
-------
8+
9+
These examples are released under the BSD open source license. Refer to license.txt in this directory for full details.

MATLAB_arrayfun/drawAntenna.m

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
function drawAntenna(p, h, dir, b, nsections)
2+
% Draws a 3D representation of an antenna
3+
4+
% Draw box on the ground
5+
base = [p+[b;0;0], p+[0;b;0], p+[-b;0;0], p+[0;-b;0], p+[b;0;0]];
6+
drawline(base);
7+
8+
% Draw edges of pyramid
9+
top = p+[0;0;h];
10+
pyramid = [top, base(:,1), top, base(:,2), top, base(:,3), top, base(:,4)];
11+
drawline(pyramid);
12+
13+
% Draw cross on base
14+
crx = [p+[b;0;0], p+[-b;0;0], p+[0;b;0], p+[0;-b;0]];
15+
drawline(crx);
16+
17+
% Draw segments
18+
if nargin < 5
19+
nsections = floor(h/b);
20+
end
21+
bottom = base;
22+
sh = h/nsections;
23+
for s = 1:nsections-1
24+
up = p + s*[0;0;sh];
25+
r = ((h-(s*sh))/h) * b;
26+
top = [up+[r;0;0], up+[0;r;0], up+[-r;0;0], up+[0;-r;0], up+[r;0;0]];
27+
28+
% Draw top of segment
29+
drawline(top);
30+
31+
% Draw crosses on each face
32+
for f = 1:4
33+
crx = [bottom(:,1), top(:,2), top(:,1), bottom(:,2)];
34+
drawline(crx);
35+
top = circshift(top,1,2);
36+
bottom = circshift(bottom,1,2);
37+
end
38+
39+
bottom = top;
40+
end
41+
42+
% Draw antenna direction
43+
l = norm(dir);
44+
45+
% Direction perp to dir horizontally
46+
perph = cross(dir, [0;0;1])/2;
47+
% Direction perp to dir vertically
48+
perpv = cross(dir, perph)/l/2;
49+
50+
% Draw 'Image' plane
51+
top = p+[0;0;h];
52+
c = top + dir;
53+
im = [c+perph+perpv, c+perph-perpv, c-perph-perpv, c-perph+perpv];
54+
drawline(im);
55+
56+
% Draw view pyramid
57+
pyramid = [top, im(:,1), top, im(:,2), top, im(:,3), top, im(:,4)];
58+
drawline(pyramid);
59+
60+
end
61+
62+
function drawline(coords)
63+
line(coords(1,:), coords(2,:), coords(3,:));
64+
end

MATLAB_arrayfun/drawMap.m

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
function drawMap(map, masts, dirn, H)
2+
% Draws a 3D map of the topological data and antenna positions
3+
4+
close all
5+
N = size(map,1);
6+
N2 = sqrt(N);
7+
M = size(masts,1);
8+
mapX = reshape(map(:,1), [N2 N2]);
9+
mapY = reshape(map(:,2), [N2 N2]);
10+
mapZ = reshape(map(:,3), [N2 N2]);
11+
mesh(gather(mapX), gather(mapY), gather(mapZ));
12+
view(-41,41);
13+
hold on;
14+
for a = 1:M
15+
mastBase = masts(a,:)' - [0;0;H];
16+
scale = 20;
17+
baseline = 10;
18+
drawAntenna(mastBase, H*scale, dirn(:,a)*scale*baseline, scale*baseline, 5);
19+
end
20+
axis off
21+
set(gca, 'Position', [0 0 1 1], 'ActivePositionProperty', 'position');
22+
end
Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
function drawMultipleSignalMaps(map, masts, AntennaDirection, H, signalMap)
2+
% Draws signal maps stored pagewise along dimension 3 of the input
3+
% signalMap
4+
5+
close;
6+
numMaps = size(signalMap,2);
7+
for m = 1:numMaps
8+
subplot(1, numMaps, m);
9+
drawSignalMap(map, masts, AntennaDirection, H, signalMap(:,m));
10+
opos = get(gca, 'OuterPosition');
11+
opos(1) = (m-1)/numMaps; % Left
12+
opos(3) = 1/numMaps; % Width
13+
set(gca, 'OuterPosition', opos);
14+
end
15+
set(gcf, 'MenuBar', 'none');
16+
oldpos = get(gcf, 'OuterPosition');
17+
newpos = [oldpos(1:2) [oldpos(3)*1.5 oldpos(4)*0.65]*0.8];
18+
set(gcf, 'OuterPosition', newpos);
19+
set(gcf, 'Position', newpos);
20+
end

MATLAB_arrayfun/drawSignalMap.m

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
function drawSignalMap(map, masts, dirn, H, signalMap)
2+
% Draws a 2D map of received power at each gridpoint as well as a
3+
% representation of the antennae.
4+
5+
hold off
6+
N = size(map,1);
7+
N2 = sqrt(N);
8+
M = size(masts,1);
9+
mapX = reshape(map(:,1), [N2 N2]);
10+
mapY = reshape(map(:,2), [N2 N2]);
11+
mapZ = reshape(map(:,3), [N2 N2]);
12+
contour(mapX, mapY, mapZ); hold on;
13+
for a = 1:M
14+
mastBase = masts(a,:)' - [0;0;H];
15+
scale = 20;
16+
baseline = 10;
17+
drawAntenna(mastBase, H*scale, dirn(:,a)*scale*baseline, scale*baseline, 5);
18+
end
19+
colormap hot;
20+
caxis(gather([min(signalMap) max(signalMap)]));
21+
surf(mapX, mapY, reshape(signalMap, size(mapX)));
22+
shading flat;
23+
grid off;
24+
axis off;
25+
end
Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
% Given the setup at the start of the article, this script divides the
2+
% antennae between 3 networks and computes a signal map for each one
3+
% independently. While this could be done using a loop, we prefer to use
4+
% vectorized operations to process all the data together.
5+
6+
% Assign network 1-3 (at random) to each antenna at every point
7+
Network = gpuArray.randi(3, 1, M);
8+
NetworkReplicated = repmat(Network, [N 1]);
9+
% Sort power data and networks by power
10+
[signalPowerSorted, I] = sort(signalPowerDecibels, 2, 'descend');
11+
NetworkSorted = NetworkReplicated(sub2ind([N M], RowIndex(:), I(:)));
12+
% Sort again to group by network
13+
[~, J] = sort(reshape(NetworkSorted, [N M]), 2);
14+
signalPowerGrouped = signalPowerSorted(J);
15+
% Sorting original list and diff finds group boundaries as an index array
16+
NetworkGrouped = sort(Network);
17+
groupStartColumn = find(diff([0, NetworkGrouped]));
18+
% Accumulated sum will allow us to find means without loops
19+
signalPowerCum = cumsum([zeros(N,1), signalPowerGrouped], 2);
20+
% Pick out the sum of the three most powerful antennae to get the mean
21+
signalMap = (signalPowerCum(:,groupStartColumn+3) - ...
22+
signalPowerCum(:,groupStartColumn)) / 3;
23+
drawMultipleSignalMaps(map, masts, AntennaDirection, H, signalMap)
Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
function signalPower = powerCalculationWithGpuArray(X, mastPos, Power, Frequency, dirn)
2+
3+
pathToAntenna = bsxfun(@minus, mastPos, X);
4+
distanceSquared = sum(pathToAntenna.^2, 3);
5+
distance = sqrt(distanceSquared);
6+
pathLoss = (4 .* pi .* distance * Frequency ./ 3e8).^2;
7+
signalPowerWatts = bsxfun(@rdivide, Power, pathLoss);
8+
signalPower = 30 + 10 * log10(signalPowerWatts);
9+
directionDotProduct = sum(bsxfun(@times, pathToAntenna, dirn), 3);
10+
isFacing = directionDotProduct >= 0;
11+
signalPower(~isFacing) = -inf;
12+
13+
end

MATLAB_arrayfun/powerKernel.m

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
function signalPower = powerKernel(X, Y, Z, mastX, mastY, mastZ, Power, Frequency, dirnX, dirnY, dirnZ)
2+
% Arrayfun kernel for computing received power at [X;Y;Z] from antenna at
3+
% [mastX; mastY; mastZ]
4+
5+
pathX = mastX - X; pathY = mastY - Y; pathZ = mastZ - Z;
6+
distanceSquared = pathX.*pathX + pathY.*pathY + pathZ.*pathZ;
7+
distance = sqrt(distanceSquared);
8+
pathLoss = (4 .* pi .* distance .* Frequency ./ 3e8).^2;
9+
signalPowerWatts = Power ./ pathLoss;
10+
signalPower = 30 + 10 * log10(signalPowerWatts);
11+
directionDotProduct = pathX*dirnX + pathY*dirnY + pathZ*dirnZ;
12+
if directionDotProduct < 0
13+
signalPower = -inf;
14+
end
15+
16+
end

MATLAB_arrayfun/reportTime.m

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
function tnow = reportTime(t, messageString, scale, tbefore)
2+
3+
if nargin < 3
4+
scale = 1;
5+
end
6+
wait(gpuDevice);
7+
tnow = t * scale;
8+
disp([messageString ' = ' num2str(tnow*1e3) 'ms']);
9+
if nargin > 3
10+
disp(['Speedup of ' num2str(tbefore/tnow) 'x']);
11+
end
12+
end

0 commit comments

Comments
 (0)