I am doing some simulations that requires me to take partial trace over a three qubit density matrix $\rho_{ABC}$. I find the mixed state density matrix of one qubit by tracing out the other two,
$\rho_A = Tr_{BC}(\rho_{ABC})$.
How can I check if my computations are right ? In other words how can I check to see if $\rho_A , \rho_B, \rho_C$ are really the mixed states coming from $\rho_{ABC}$ ?
Your problem is in some sense similar to the reverse problem of the partial trace: Given a selection of density matrices, does there exist a state (and which is it) with these marginals? This problem is known as the quantum marginal problem - a very hard problem that is (for example if the marginals are overlapping) not completely solved.
It is easy to see that the quantum marginal problem with non-overlapping marginals (for example, when you are given the one-body reduced densities) only depends on the eigenvalues of the reduced density matrices since you can locally diagonalize any state. So what you really search are eigenvalue inequalities. One such inequality is clear: The trace of the reduced density matrices and therefore the sum of all eigenvalues must be equal to one for normalized states.
In principle, you can give a number of eigenvalue equations for any dimension, but finding them is extremely hard. See Klyatchko's solution for a proof and an appendix listing all inequalities for certain low-dimensional systems (e.g. three qubits, where the complete solution consits of 40 independent inequalities). Also see Walter's thesis for a more modern approach and discussion of related problems. You will immediately see that it is insane to use the solution to the quantum marginal problem as a sanity check, not only because calculating the complete spectra will be much more prone to errors than just taking partial traces, but also because it doesn't actually tell you whether your specific reduced states are the reduced density matrices of your specific state.
So using the quantum marginal problem is out of the question.
Your problem is much more specific than that. However, except for obvious sanity checks (trace, positivity), the only thing I see that you could possibly do is check the definition: Take local observables $X$ and calculate e.g. $\operatorname{tr}(\rho_{ABC}X\otimes \operatorname{id} \otimes \operatorname{id})$ and $\operatorname{tr}(\rho_AX)$. If you take a basis of $X$, this gives you a second way to calculate the reduced density matrices, however, it is much more clumsy then the usual way. This is, because the partial trace is among the simplest calculations immaginable (from a computational viewpoint).
To put it differently, how do you check that you calculated the trace correctly? You could of course calculate the characteristic polynomial and check the second coefficient, but wouldn't it be far more likely you made a mistake calculating this?
Therefore, I doubt that there is any good way to check your partial trace calculations...
