I'm calculating Eigenvalues for Dense Matrices using the .eigenvalues() function of MatrixXd. However running the program multiple times I start getting different results for matrices of size 1024x1024 and larger. Is there a way to increase the precision of the eigenvalues function?

If you get different results, this means the input matrix is different. There is no randomness in the computation of the eigenvalues. Please, be more specific about your issue.

Ok that is interesting I'll have to double check my code. It doesn't appear to contain randomness on my end either that is why it is strange that apparently the eigenvalues I get are different. Is the ordering of the eigenvalues random?

Eigenvalues are not sorted, but the same matrix must output the eigenvalues in same order.

BTW, if your matrix is symmetric you should use SelfAdjointEigenSolver<MatrixXd>(A).eigenvalues() which is faster and produces sorted eigenvalues because in this case we know they are real.