invert function Null safety
- Matrix m
Invert m in place.
Implementation
void invert(final Matrix m) {
if (!m.isLU) {
decomposeLU(m);
}
// Clone the matrix to a native double matrix, for better internal precision.
// This seems to give results closer to the real 15C than the version that
// did the internal math using Value's precision, in a quick test. I suspect
// that the 15C may be using a more clever algorithm, but brute force works,
// too!
final dm = List<List<double>>.generate(m.rows,
(row) => List<double>.generate(m.columns, (col) => m.getF(row, col)));
/// Now use A^-1 = U^-1 * l^-1 * P, as per HP 15C Advanced Functions p. 83
// Calculate U^-1. Adapted from dtri2.f in LAPACK from www.netlib.org.
for (int j = 0; j < m.rows; j++) {
final ajj = -1 / dm[j][j];
dm[j][j] = -ajj;
// Compute elements 0..j-1 of the jth column
// DTRMV call:
for (int jj = 0; jj < j; jj++) {
final temp = dm[jj][j];
for (int i = 0; i < jj; i++) {
dm[i][j] = dm[i][j] + temp * dm[i][jj];
}
dm[jj][j] = dm[jj][j] * dm[jj][jj];
}
// DSCAL call:
for (int i = 0; i < j; i++) {
dm[i][j] = dm[i][j] * ajj;
}
}
// Calculate L^-1, adapted from dtri2.f.
for (int j = m.rows - 2; j >= 0; j--) {
const ajj = -1;
// DTRMV call:
for (int jj = m.rows - 2 - j; jj >= 0; jj--) {
final temp = dm[j + jj + 1][j];
for (int i = m.rows - 2 - j; i > jj; i--) {
dm[j + 1 + i][j] = dm[j + 1 + i][j] + temp * dm[j + i + 1][j + jj + 1];
}
}
// DSCAL call:
for (int i = j + 1; i < m.rows; i++) {
dm[i][j] = dm[i][j] * ajj;
}
}
// Calculate m = U^-1 dot L^-1 in-place:
for (int r = 0; r < m.rows; r++) {
for (int c = 0; c < m.columns; c++) {
double v = 0;
for (int k = max(r, c); k < m.columns; k++) {
assert(r <= k); // Otherwise U is zero
assert(c <= k); // Otherwise L is zero;
final uv = dm[r][k];
final lv = (k == c) ? 1 : dm[k][c];
v += uv * lv;
}
dm[r][c] = v;
}
}
// Now copy back into m...
m.visit((r, c) => m.setF(r, c, dm[r][c]));
m.dotByP();
m.isLU = false;
}