in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/SquareMatrixSupport.java [235:274]
private static int multiply(double[] a, int sa, double[] b, int sb, double[] col, double[] out) {
final int m = col.length;
// Rows are contiguous; Columns are non-contiguous
int k;
for (int c = 0; c < m; c++) {
// Extract column from b to contiguous memory
k = c;
for (int i = 0; i < m; i++, k += m) {
col[i] = b[k];
}
// row * col
k = 0;
for (int r = 0; r < m; r++) {
double sum = 0;
for (int i = 0; i < m; i++, k++) {
sum += a[k] * col[i];
}
out[r * m + c] = sum;
}
}
int s = sa + sb;
// Overflow protection. Ideally we would check all elements but for speed
// we check the central one only.
k = m >> 1;
if (out[k * m + k] > SCALE_THRESHOLD) {
// Downscale
// We could downscale by the inverse of SCALE_THRESHOLD.
// However this does not account for how far above the threshold
// the central element is. Here we downscale so the central element
// is roughly 1 allowing other elements to be larger and still protected
// from overflow.
final int exp = Math.getExponent(out[k * m + k]);
final double downScale = Math.scalb(1.0, -exp);
s += exp;
for (int i = 0; i < out.length; i++) {
out[i] *= downScale;
}
}
return s;
}