Construct a matrix by accumulation of indexed values (C, C++, or Fortran)

INTRO: The function 'accumarray' in MATLAB is too slow for my codes (because they use large vectors and matrices), so I am seeking a FASTER version satisfying my needs. I am seeking the fastest possible implementation. GOAL: Write a C, C++, or Fortran code which performs part of the functionality of the MATLAB function "accumarray". Your code must then be made into a MEX function for calling from MATLAB. DETAILS: Your function must reproduce the following behavior (or similar): %% GIVEN EXAMPLE INPUTS (This is MATLAB code): numvalues = 10^6; % number of values to accumulate matSize = 288; % width of matrix vals = ones(numvalues,1); % values to be accumulated xn = uint32(randi(matSize,numvalues,1)); % x index for values yn = uint32(randi(matSize,numvalues,1)); % y index for values zn = uint32(randi(matSize,numvalues,1)); % z index for values pos3 = yn+(xn-1)*matSize+(zn-1)*matSize^2; % only needed for alternative code %% YOUR FUNCTION WILL GIVE THE SAME OUTPUT THAT MATLAB GIVES USING ACCUMARRAY AS FOLLOWS: newMat3 = accumarray([yn, xn, zn], vals, [matSize matSize matSize]); % a 3D evaluation newMat2 = accumarray([yn, xn], vals, [matSize matSize]); % a 2D evaluation newMat1 = accumarray(yn, vals, [matSize 1]); % a 1D evaluation %% The above form is preferred, but the following alternative is also acceptable if it improves speed in your function: newMat3 = accumarray(pos3,vals,[matSize^3 1]); If possible, please report your expected speedup times. I am also open to other ideas about how to speedup this function.