Logarithm of singular matrices

I wanted to calculate following quantity:

$X=Tr[\rho_1 \log[\rho_2]]$,

as in the relative entropy. Here, $\rho_1, \rho_2$ are positive semidefinite matrices with non-orthogonal support (so that the thing does not diverge - terms like $0 \log[0]$ should be taken to zero, which is standard assumption) and $Tr$ is trace.

There is a similar SE question, the given function

MatrixLogSafe[x_] := MatrixFunction[Piecewise[{{Log[#1], #1 > 0}}] &, x] 

which should deal with the matrix logarithm, it behaves, however, strange.

For e.x. I assume that $\rho_1 =\rho_2 =\{\{0.33,0,0\},\{0,0,0\},\{0,0,0.66\}\}$. The quantity $X$ should then be

$X=0.33 \log[0.33] + 0.66 \log[0.66]= -0.640099$.

However, using MatriLogSafe in the definition gives different output:

In[402]:= Tr[{{0.33, 0, 0}, {0, 0, 0}, {0, 0, 0.66}}.MatrixLogSafe[{{0.33, 0, 0}, {0, 0, 0}, {0, 0, 0.66}}]] Out[402]= -0.731717 

The problem is, that MatrixLogSafe sometimes "switch the eigenvectors",

In[403]:= MatrixLogSafe[{{0.33, 0, 0}, {0, 0, 0}, {0, 0, 0.66}}] Out[403]= {{0., 0., 0.}, {0., -0.415515, 0.}, {0., 0., -1.10866}} 

(so $\log[0.33]= -1.10866$ and $\log[0.66]=-0.415515$, but the output should be { { -1.10866, 0, 0.}, {0., 0., 0.}, {0., 0.,-0.415515}}).

(Somehow I think the problem is that the I use numerical values, but I want that the function works for both numerical and "exact" (?) numbers)

How one can fix it?

---

I have consider the answer given by Carl Woll, however still something is not working. In particular, consider two matrices, $\rho_1$:

{{1/4, 1/4 E^(-((I \[Pi])/10)), 1/4 E^(-((I \[Pi])/10)), 1/ 4}, {1/4 E^((I \[Pi])/10), 1/4, 1/4, 1/4 E^((I \[Pi])/10)}, {1/4 E^((I \[Pi])/10), 1/4, 1/4, 1/4 E^((I \[Pi])/10)}, {1/4, 1/4 E^(-((I \[Pi])/10)), 1/4 E^(-((I \[Pi])/10)), 1/4}} 

and $\rho_2$:

{{1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E]), 1/4 E^(-(1/2) - (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10)), 1/4 E^(-(1/2) - (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10)), 1/(4 Sqrt[ E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10))^2}, {1/4 E^(-(1/2) + (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10)), 1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E]), 1/(4 Sqrt[ E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10))^2, 1/4 E^(-(1/2) + (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10))}, {1/4 E^(-(1/2) + (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10)), 1/(4 Sqrt[ E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10))^2, 1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E]), 1/4 E^(-(1/2) + (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10))}, {1/( 4 Sqrt[E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10))^2, 1/4 E^(-(1/2) - (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10)), 1/4 E^(-(1/2) - (I \[Pi])/10) + 1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) + 1/4 E^((I \[Pi])/10)), 1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E])}} 

which are maybe ugly, but they both are Hermitian and positive-semi-definite.

I get however complex $X$ ($X=densityTrace[\rho_1,\rho_2]=-0.0019613 + 0.393667 I$), which cannot be!

However, if I calculate firstly the numerical matrices I get real output $densityTrace[N[\rho_1],N[\rho_2]]=-0.0432473$ (?!)

I have slightly changed the definition of densityTrace

densityTrace[a_, b_] := Module[{λ, S, d}, {λ, S} = Eigensystem[b]; S = Transpose[S]; d = Diagonal[Inverse[S]. a. S]; Total @ MapThread[If[Chop[#1]==0,0,Chop[#1] Log[Chop[#2]]]&, {d, λ}] 

]

(adding Chop), to get rid of some very small imaginary "waste", is it a good idea?