How to do Lowdin symmetric orthonormalisation?
You obtain the diagonalization of $S$ as defined
$$
S = U \cdot s  \cdot U^{\dagger}
$$
(which holds for the case of symmetric $S$, which it is) and apply the power needed
$$
S^{-\frac{1}{2}} = U \cdot s^{-\frac{1}{2}}  \cdot U^{\dagger}
$$

For more background, try the application section of Wikipedia's article on diagonalizable matrices.
As far as understanding it goes, I'm afraid not much beyond doing some model calculations on 2D or 3D vectors can be done.

I do not have a strong background in parallel computing, but all operations involved are well-known and your linear algebra framework (BLAS, anyone?) will document how to run any of this in parallel and what storage may be needed for it. Note that, typically, diagonalization is not the most expensive step of a self-consistent field calculation (Hartree-Fock or DFT) anyway, so parallelization may not be critical for this step.