invert function Null safety

void invert(
  1. 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;
}