The residual MK*V-MM*V*w2 is expected to be numerically close to zero, however, you have to take into account the scaling of the whole system. Here norm(MK) is 4.2241e+20 and the residual will also be scaled with this factor, which gets us to a residual of around 1e7. Here's what I get for a scaled residual:
>> [V, w2]=eig(MK,MM);
>> max(vecnorm(MK*V - MM*V*w2)) / norm(MK)
About matching eigenvectors between MATLAB's and Abaqus's result - you definitely have to expect different scaling of each individual eigenvector (you could try applying V ./ vecnorm(V) to both matrices, which will given each column a 2-norm of 1). At that point, the only effect of scaling should be the sign of each column.
But in addition to this, here's a plot of the eigenvalues w2:
As you can see, the eigenvalues are (A) not sorted, so the order of the eigenvalues returned by Abaqus might be different, and (B) decrease rapidly, so the eigenvectors related to the relatively very small eigenvalues (< 1e7 let's say) probably don't have matching eigenvectors, because at that point they are mostly tracking numerical noise.