Do commuting matrices with real and complex eigenvectors still share eigenvectors? Main Question
Say that $K$ and $S$ are two commuting matrices with a full set of real eigenvectors, where the dimension of the eigenspaces of the matrices may be greater than one.
$$KS=SK$$
It is a common theorem in linear algebra that $S$ and $K$ will share eigenvectors. The top answer in this question explains it beautifully:
Matrices commute if and only if they share a common basis of eigenvectors?
However, say that the eigenvectors of $K$ are real, while the eigenvectors of $S$ are complex.
In this case, will the matrices still share eigenvectors?
Or more specifically, will we allways be able to linearly combine the complex eigenvectors of $S$ to get the real eigenvectors of $K$?

Reason for Asking:
(Ignore the physics part if this doesn't ring a bell).
In physics textbooks deriving the normal-modes of coupled oscillators, its common practice to introduce a "symmetry" matrix transformation. Instead of finding the eigenvectors of the matrix describing the coupling of the system $K$, the eigenvectors of the symmetry matrix $S$ are found, and those correspond to the normal-modes of the actual system.
The theorem that's used to justify the fact that $K$ and $S$ share eigenvectors is that they always commute.
However, often the eigenvectors of $S$ are complex, while the eigenvectors of $K$ are real. In this case, the complex eigenvectors of $S$ are linearly combined to give real vectors, and those turn out to be eigenvectors of $K$.
Can someone explain/prove why this is? Or point me to some useful theorems about commuting matrices when the eigenvectors of one are complex, while the eigenvectors of the other are real?
Thank you!

Example:
Consider the two matrices:
$$K=\left(\begin{array}{ccc}
-2 & 1 & 1 \\
1 & -2 & 1 \\
1 & 1 & -2
\end{array}\right) S=\left(\begin{array}{lll}
0 & 1 & 0 \\
0 & 0 & 1 \\
1 & 0 & 0
\end{array}\right)$$
$$$$
These two matrices commute, to give:
$$SK=KS=\left(\begin{array}{lll}
1 & -2 & 1 \\
1 & 1 & -2 \\
-2 & 1 & 1
\end{array}\right)$$
The eigenvectors of $S$ would be complex. We can solve for them by first solving for the eigenvalues, $\lambda_k$, of which there should be three.
$$\begin{array}{c}
\operatorname{det}\left|\mathrm{S}-\mathrm{\Lambda}_{\mathrm{k}}\right|=0 \\
\operatorname{det}\left[\begin{array}{ccc}
-\lambda_{k} & 1 & 0 \\
0 & -\lambda_{k} & 1 \\
1 & 0 & -\lambda_{k}
\end{array}\right]=0 \\
-\lambda_{k}^{3}=-1 \\
\lambda_{k}^{3}=1
\end{array}$$
And we need to solve for the cubic roots of $1 .$ There are three of those:
$$
\lambda_{0}=1, \lambda_{1}=\frac{2 \pi}{3}, \lambda_{2}=\frac{4 \pi}{3}
$$
If we assume that the first component of each of the eigenvectors is $1$, that gives us three distinct eigenvectors:
$$\lambda_{0}=1 \Rightarrow \overrightarrow{A^{(0)}}=\left[\begin{array}{l}1 \\ 1 \\ 1\end{array}\right], \lambda_{1}=e^{i \frac{2 \pi}{3}} \Rightarrow \overrightarrow{A^{(1)}}=\left[\begin{array}{c}1 \\ e^{i \frac{2 \pi}{3}} \\ e^{i \frac{4 \pi}{3}}\end{array}\right], \lambda_{2}=e^{i \frac{4 \pi}{3}} \Rightarrow \overrightarrow{A^{(2)}}=\left[\begin{array}{c}1 \\ e^{i \frac{4 \pi}{3}} \\ e^{i \frac{8 \pi}{3}}\end{array}\right]$$
If we write out the complex and real parts of $S$, we get:
$$\overrightarrow{A^{(0)}}=\left[\begin{array}{l}
1 \\
1 \\
1
\end{array}\right]$$
$$\overrightarrow{A^{(1)}}=\left[\begin{array}{c}
1 \\
-0.5 \\
-0.5
\end{array}\right] + \left[\begin{array}{c}
0 \\
.866i\\
-.866i
\end{array}\right]$$
$$\overrightarrow{A^{(2)}}=\left[\begin{array}{c}
1 \\
-0.5 \\
-0.5
\end{array}\right] + \left[\begin{array}{c}
0 \\
-.866i\\
+.866i
\end{array}\right]$$
And, I'll leave it up to you to verify that the complex and real parts of each of those eigenvectors are eigenvectors of $K$ as well.
 A: Your problem is virtually trivial when your write it properly. S is the celebrated cyclic permutation matrix,
$$S=\left(\begin{array}{lll}
0 & 1 & 0 \\
0 & 0 & 1 \\
1 & 0 & 0
\end{array}\right)~~\leadsto S^T=S^{-1},$$
with the famed eigenvalues and eigenvectors
$$
1, ~~\left[\begin{array}{l}
1 \\
1 \\
1
\end{array}\right]; ~~~~~ \omega, ~~\left[\begin{array}{l}
1 \\
\omega \\
\omega^2
\end{array}\right]; ~~~~~ \omega^2 , ~~\left[\begin{array}{l}
1 \\
\omega^2 \\
\omega
\end{array}\right],
$$
where $\omega=\exp (i2\pi/3)$, so $\omega^3=1$, and $1+\omega+ \omega^2=0$.

*

*The last two are of course complex conjugate to each other, both eigenvalues and eigenvectors!

Now,
$$
K= -2 {\mathbb I} +S + S^2,
$$
so of course it manifestly commutes with S, and shares its 3 eigenvectors, now with eigenvalues  0, -3, -3, respectively.
The first is real, and you may take real and imaginary parts by considering the sum and difference of the latter two, and chuck the i of the imaginary piece. For this to work you need degeneracy, as you have here.
A: 2 different answers depending on what OP actually wants to ask here.

will the matrices still share eigenvectors?
Or more specifically, will we always be able to linearly combine the
complex eigenvectors of  to get the real eigenvectors of ?

First Answer
There are some significant ambiguities in the original question.  OP's question "will we always be able to linearly combine the complex eigenvectors of  to get the real eigenvectors of ?"
is answered yes but in a somewhat trivial way if S is diagonalizable -- by definition $S$'s eigenvectors form a basis for the n dim vector space so any eigenvector of $K$ can be written as a linear combination of the eigenvectors of $S$.  This has nothing to do with commutativity.

I suspect OP wants a 'sharper' answer that leans on commutativity which is what I give below.  In short, subject to a couple of mild assumptions, any eigenvector of $K$ may be written as a linear combination of eigenvectors of $S$ that have $real(\lambda) = a$ for some $a\in \mathbb R$.
Second Answer
I assume that $S$ is normal; then the answer is yes.
For simplicity assume that $K$ is diagonalizable.  (Diagonalizability of $K$ doesn't seem crucial though introducing defective matrices leads to clumsier results.  On the other hand normality of a $S$ however plays a key role in the below.)
1.) $SK=KS\longrightarrow S^*K=KS^* $
where $^*$ denotes conjugate transpose.
Proof:  Fuglede's theorem in finite-dimensional vector space
2.) consider the symmetrization map
$S\mapsto \frac{1}{2}\Big(S+S^*\Big):= H$
$H$ is real symmetric and has eigenvectors that may be selected to be real orthonrmal.
Note: if $\lambda = a+bi$ is an eigenvalue of $S$ then $\lambda\mapsto a$.
3.) Let $H$ have distinct eigenvalues $\{\lambda_1,\lambda_2,...,\lambda_r\}$, and specialize to some $\lambda_i$ in that set.
Since $S$ is normal, it is unitarilly diagonalizable.
$Q^*SQ = D$ which tells us
$m=\dim\ker\Big(H-\lambda_iI\Big)= \dim\ker\Big(Q^*\big(H-\lambda_iI\big)Q\Big)= \dim\ker\Big(\frac{1}{2}\big(D+D^*\big)-\lambda_iI\Big)$
The sum of diagonal matrices is diagonal, so this tells us that that we have
$\Big(\frac{1}{2}\big(D+D^*\big)-\lambda_iI\Big)$
with exactly $m$ zeros on the diagonal. WLOG (i.e. up to re-labelling) assume the first $m$ entries are zero.  Thus the first m standard basis vectors span the above kernel, and after 'undoing' the unitary similarity transform, we see that $\{\mathbf q_1, \mathbf q_2,...,\mathbf q_m\}$ spans $\ker\Big(H-\lambda_iI\Big)$. The $\mathbf q_k$ are orthonormal eigenvectors of $S$ that are not entirely real, yet $H$ being real symmetric has eigenvectors that may be chosen to be entirely real hence there must be some linear combination of those $\mathbf q_k$ that results in $m$ real linearly indendent vectors that spans $\ker\Big(H-\lambda_iI\Big)$.
4.) By (1) we know $HK=KH$ and since both matrices are diagonalizable they are simultaneously diagonalizable.  This means we may select $n$ linearly independent eigenvector $\{\mathbf z_1,\mathbf z_2,...,\mathbf z_n\}$ where $\mathbf z_k$ is an eigenvector of $K$ and $\mathbf z_k \in \ker\big(H-\lambda_iI\big)$ for some $i$.  By (3) each $\mathbf z_k$ may be written as a linear combination of the complex eigenvectors of $S$ which completes the proof.
