function hpccLinpack( n )
%HPCCLINPACK An implementation of the HPCC Global HPL benchmark
%
%  hpccLinpack(m) creates a random codistributed real matrix A of size 
%    n-by-n and a real random codistributed vector B of length n. It then
%    measures the time to perform the matrix division of A into B (X = A\B,
%    which is the solution to the equation A*X = B) in a parallel way using
%    the currently available resources (MATLAB pool). This time indicates
%    the performance metric. Finally the function computes the scaled
%    residuals to ensure that the error on the computation is within
%    acceptable bounds.
%
%    If you do not specify m, the default value is that returned from
%    hpccGetProblemSize('hpl'), which assumes that each process in the pool
%    has 256 MB of memory available. This is expected to be smaller than
%    the actual memory available. 
%
%    Details of the HPC Challenge benchmarks can be found at
%    www.hpcchallenge.org and the specific Class 2 specs are linked off
%    that page. (At the time of writing, the specs are linked at
%    www.hpcchallenge.org/class2specs.pdf.)
%
%    Examples:
%
%      % Without a parallel pool open
%      tic; hpccLinpack; toc
%      Data size: 0.108665 GB
%      Performance: 16.351622 GFlops
%      Elapsed time is 2.791896 seconds.
%
%      % With a local parallel pool of size 4
%      tic; hpccLinpack; toc
%      Data size: 0.434774 GB
%      Performance: 18.650758 GFlops
%      Elapsed time is 21.647003 seconds.
%
%    See also: hpccGetProblemSize, parpool

%   Copyright 2008-2018 The MathWorks, Inc.

% If no size provided then get a default size
if nargin < 1
    n = hpccGetProblemSize( 'hpl' );
end

spmd
    % Create a distributed matrix in the 2d block cyclic distribution and a
    % distributed column vector in 1d
    A = codistributed.randn(n, n, codistributor2dbc);
    b = codistributed.rand(n, 1);
    
    % Time the solution of the linear system   
    tic
    x = A\b;
    t = toc;    

    % Need to convert to a 1d distribution for the checking code below
    A = redistribute(A, codistributor1d);
    % Compute scaled residuals
    r = norm(A*x-b,inf)/(n*eps*(norm(A,inf)*norm(x,inf)+norm(b,inf)));
    % This test is specified in the benchmark definition
    if r > 16
        error('Failed the HPC HPL Benchmark');
    end
end
% Performance in gigaflops
perf = (2/3*n^3 + 3/2*n^2)/max([t{:}])/1.e9;

fprintf('Data size: %f GB\nPerformance: %f GFlops\n', 8*n^2/(1024^3), perf);

