eigh delivers the eigenvectors as the columns of the returned array. Essentially, the current code results in a wrong eigendecomposition of the form E^T Diag(evals) E instead of E Diag(evals) E^T, with E the matrix containing the evecs in its columns. This completely jumbles up the samples.

In order to correct this, we just need to transpose the _factor array. I'll put up a PR momentarily.

Reproducing code example:

import numpy as np
from numpy.random import default_rng
rg = default_rng()

mu = np.array([1., 1., 1.2])
Sigma = np.array([[1., 1., 0.8], [1., 1., 0.8], [0.8, 0.8, 1.]])

rg.multivariate_normal(mu, Sigma)
array([1.36821232, 1.36821231, 1.6910961 ])  # SVD works fine

rg.multivariate_normal(mu, Sigma, method="eigh")
array([1.00000002, 0.72847392, 1.02539934])  # this is obviously wrong.

UPDATE:
Looks like samples from the 'cholesky' method were also incorrect, see #15872 (comment)

Numpy/Python version information:

1.18.1 3.7.7 (default, Mar 26 2020, 10:32:53)
[Clang 4.0.1 (tags/RELEASE_401/final)]