I have n-by-n correlation matrix $C_{ij}$, it's symmetric positive definite, and $C_{ii}=1$ and $C_{ij}\le1$ etc.

I want to generate a bunch of matrices $M_{ij}$ which are on the same distance from this initial matrix. The distance measure is Riemannian, like in this answer: $$d(C,M)=\sqrt{\sum_{k=1}^n \log^2 \lambda_k}$$ where $\lambda_k$ are eigenvalues of $CM^{-1}$

I'm looking for an algorithm that can find me the matrices $M$.