for (int k = 0; k < n; k++) {
for (int j = 0; j < nx; j++) {
BigDecimal s = BigDecimal.ZERO;
for (int i = k; i < m; i++) {
int scale = MatrixBigDecimal.GLOBAL_SCALE;
s = s.add(QR[i][k].multiply(X[i][j]), new MathContext(scale));
}
int scale = Math.max(s.scale(), QR[k][k].scale());
s = s.negate().divide(QR[k][k], scale, roundingMode);
for (int i = k; i < m; i++) {
int scale2 = MatrixBigDecimal.GLOBAL_SCALE;
X[i][j] = X[i][j].add(s.multiply(QR[i][k], new MathContext(scale2)));
}
}
}
// Solve R*X = Y;
for (int k = n-1; k >= 0; k--) {
for (int j = 0; j < nx; j++) {
int scale = MatrixBigDecimal.GLOBAL_SCALE;
X[k][j] = X[k][j].divide(Rdiag[k], scale, roundingMode);
}
for (int i = 0; i < k; i++) {
for (int j = 0; j < nx; j++) {
int scale = MatrixBigDecimal.GLOBAL_SCALE;
X[i][j] = X[i][j].subtract(X[k][j].multiply(QR[i][k], new MathContext(scale)));
}
}
}
return (new JamaMatrixBigDecimal(X,n,nx).getMatrix(0,n-1,0,nx-1));
}