private static int multiply()

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;
        }