function dec = addMonomialMatrixDec(dec, nMax)

%% generate monomial matrices
mat = cell(nMax, 1);
% mat{1} = [];%zeros(1, dec.nC);
nMonomial = 0;

for iExpo = 1:nMax
    % obtain indices for exponent
    mat{iExpo} = obtainOrderedSetMatrix(dec.nC, iExpo, 1);
    nRow = size(mat{iExpo}, 1);
    nMonomial = nMonomial + nRow;
    
%     % exponent index matrix
%     mat{iExpo+1} = zeros(nRow, dec.nC);
%     for iRow = 1:nRow
%         mat{iExpo+1}(iRow
%     
%     mat = [mat; 
end

nMonomial = nMonomial + 1;
dec.monomialMatrix = zeros(nMonomial, nMax);
iRow = 0;
for iExpo = nMax:-1:1
    dec.monomialMatrix(iRow + (1:size(mat{iExpo}, 1)), 1:iExpo) = mat{iExpo};
    iRow = iRow + size(mat{iExpo}, 1);
end


%% construct derivative matrix

% will be: monomialDerivativeMatrix
matIndex = zeros(nMonomial, nMonomial);
matFactor = zeros(nMonomial, nMonomial);

for iExpo = 1:nMax
    % we get the partial derivative monomial matrix
    derivativeMonomialMatrix = dec.monomialMatrix(:, ...
        [(1:(iExpo-1)) ((iExpo+1):end)]);
    % we add a zero column to obtain the same configuration as in
    % monomialMatrix
    derivativeMonomialMatrix = [derivativeMonomialMatrix zeros(nMonomial, 1)];
    
    % now we need to find the index to monomialMatrix that matches each
    % line in derivativeMonomialMatrix
    colVector = zeros(nMonomial, 1);
    for iMono = 1:nMonomial
        filterMatrix = bsxfun(@eq, dec.monomialMatrix, ...
            derivativeMonomialMatrix(iMono, :));
        filterVector = sum(filterMatrix, 2) >= size(dec.monomialMatrix, 2);
        if sum(filterVector) > 1 || ~any(filterVector)
            error('elk:decMonomial:internalError', ['The construction of the ' ...
                'monomial derivatives failed. The partial derivative should ' ...
                'be unique in the monomialMatrix but is not. This should not ' ...
                'happen.']);
        end
        colVector(iMono) = find(filterVector);
    end
    % row vector is obvious
    rowVector = (1:nMonomial)';
    
    % but a partial derivative where monomialMatrix has a zero element
    %  should not be considered
    thisFilter = dec.monomialMatrix(:, iExpo) == 0;
    colVector(thisFilter) = [];
    rowVector(thisFilter) = [];
    
    % convert to indices into derivative matrix data
    matIndexVector = sub2ind(size(matIndex), rowVector, colVector);
    
    % add required indices - note that indices might be written multiple
    %   times but only identical indices can be written for each row
    matIndex(matIndexVector) = dec.monomialMatrix(~thisFilter, iExpo);
    
    % accumulate factor
    matFactor(matIndexVector) = matFactor(matIndexVector) + 1;
end

dec.monomialDerivativeIndexMatrix = matIndex;
dec.monomialDerivativeFactorMatrix = matFactor;
