NumCosmoMathMatrixchol_chi2_cols

Computes the squared Mahalanobis distance of each column of cm from theta under the covariance $C = U^\intercal U$, $$\chi^2_p = (x_p - \theta)^\intercal C^{-1} (x_p - \theta),$$ and stores it in the first $n_p$ entries of chi2. The triangular solve $y_p = (x_p - \theta) U^{-1}$ and the norm $|y_p|^2$ are done in one pass, on blocks of $n_b$ columns, with work as the only scratch.

Only the upper triangle of U is read, its diagonal taken as stored. Since work belongs to the caller, several threads can share cm and U, each with its own work. The four arguments must be distinct objects that do not overlap.