Let $PJP^{-1}$ denote the Jordan decomposition of $M$. The matrix $J$ is a direct sum of Jordan blocks; it is unique up to permutation of the Jordan blocks. The matrix $P$ is not unique.
There are two questions here, which are similar enough that they deserve to be asked together:
- Can the Jordan decomposition be computed in a backwards stable way?
- The Jordan decomposition is often avoided in numerical linear algebra, even though it is a very convenient decomposition for developing a lot of theory. My question is: Why?
One unconvincing answer to 2 is that it's discontinuous: The matrices $P$ and $J$ vary discontinuously as $M$ varies. But this is a bad explanation because it also applies to the SVD (note that the factors $U$ and $V$ vary discontinuously in the SVD), but the SVD still gets used in numerical linear algebra.
[EDIT: Please do not tell me that in the SVD, the map $M \mapsto \Sigma$ is Lipschitz continuous (and not just $\frac1 2$-Hölder as one answer claims) because note that people in numerics often want to compute the whole SVD, including the matrices $U$ and $V$. And that's discontinuous.]
I somehow suspect that the problem might lie with $P$ being non-unitary.