Efficiently sampling a uniform random matrix

I need to sample a huge random matrix, whose size is 499 x 15500, i.e. 7734500 elements. For this reason, I would like the sampling procedure to be as efficient as possible. At the moment, in a cpp file I am doing this:

std::random_device rd; // seed generator

std::mt19937_64 generator{rd()}; // generator initialized with seed from rd


std::uniform_int_distribution<> initialize(unsigned long long int modulus)
{
    std::uniform_int_distribution<> unifDist{0, (int)(modulus-1)};
    return unifDist;
}


Matrix<unsigned int> uniformRandomMatrix
    (unsigned int rows, unsigned int columns, unsigned long long int modulus)
{
    std::uniform_int_distribution<> dist = initialize(modulus);

    // Declare and allocate the matrix
    Matrix<unsigned int> matrix(rows, columns);
    // this constructor just does a resize on a std::vector

    // Fill the matrix with random elements
    for(unsigned int i = 0; i < rows; ++i)
    {
        for(unsigned int j = 0; j < columns; ++j)
        {
            matrix.setElementAt(i, j, dist(generator));
            // setElementAt just does matrix[somePosition] = newElement
        }
    }

    return matrix;
}

Note that a Matrix is implemented as a 1D std::vector for efficiency.

Can I do better than this? Right now, sampling this huge matrix takes approximately 160000 microseconds.