I have $n$ arbitrary $p\times 1$ vectors $x_i$, and $p\times k$ matrices $A_i$, and $n$ $p \times p$ positive semidefinite matrices $S_i$, where some (often most) of the $S_i$'s are same (for example only two different $S$ matrices, one positive definite which applies to $i=1,\ldots n-1$ and semidefinite $S$ for $i=n$). I want to make following linear transformations for all $x_i$ and $A_i$:
$x_i^* = L_i^{-1}x_i$ and $A_i^* = L_i^{-1}A_i$, where
$L_i$ is a lower triangular matrix so that $L_i L'_i=S_i$ (or even better $L_i D_i L'_i=S_i$ where $D_i$ is diagonal matrix.)
$n$ can be range from few to thousands or even more, $p$ is usually less than 10 and $k$ is usually less than 100. $S$ can contain rows and columns with zeros, and it is originally formed as $BB'$ where $B$ is a lower tridiagonal matrix with nonnegative diagonal elements.These $B$'s are not available though.
I'm wondering what would be the optimal way of accomplishing this in terms of speed and accuracy, with more weight on accuracy?
Currently I'm using my self written LDL decomposition function for different $S$, invert $L$ and compute the transformations, as I have beem dealing with cases with small number of different $S$ and large $n$. If I have understood correctly, it would be wiser to just solve the linear equations without explicitly inverting the matrix in terms of accuracy, but for speed it seems to be better option?
I'm using Fortran, and I would rather use some LAPACK function for the factorization instead of my own LDL decomposition which I cannot guarantee to be stable etc., but I'm not sure which would be the good way.
edit: This may be better suited for stackoverflow, so this can be removed by whoever has the power to do so.
