function out = hpccDGEMM(m)
%hpccDGEMM HPCC DGEMM Benchmark computation
%   hpccDGEMM(m) computes C = beta * C + alpha * A * B where
%   A, B, and C are m-by-m random real matrices, and alpha and
%   beta are random real scalars.

% Copyright 2018 The MathWorks, Inc.

if nargin < 1
    sz = hpccDataSizes();
    m = sz.DGEMM;
end

dataSizeBytes = 8*m^2;
hpccBenchmarkStarting("DGEMM", dataSizeBytes);

spmd
    % Create distributed matrices in the default 1d distribution scheme
    A   = codistributed.rand(m, m);
    B   = codistributed.rand(m, m);
    Cin = codistributed.rand(m, m);

    alpha = A(ceil(m / 2));
    beta  = B(ceil(m / 2));
    
    % Time the solution of the matrix product
    tic
    Cout = beta * Cin + alpha * A * B;
    t = toc;

    % Verify the solution. Use a sequence of matrix-vector products to compute an effective scaled
    % residual.
    x = codistributed.rand(m, 1);
    
    cnrm = norm(Cout, Inf);
    xnrm = norm(x, Inf);
    
    y = Cout * x;
    z = B * x;
    y = alpha * A * z - y;
    y = beta * Cin * x + y;

    sres = norm(y, Inf) / cnrm / xnrm / m / eps;
    if sres > 16
        error('Failed the HPC DGEMM Benchmark');
    end
end

% Performance in gigaflops
tmax = max([t{:}]);
perf = 2e-9 * m^3 / tmax;
out  = hpccBenchmarkEnd("DGEMM", dataSizeBytes, tmax, perf, "GFlops");
end
