decomposeLU function Null safety

void decomposeLU(
  1. Matrix m
)

Do an LU decomposition with row permutations, and with perturbations, if needed, to avoid a singular matrix. This is hoped to be compatible with the 15C's LU decomposition using the Doolittle method, as mentioned in the HP 15C Advanced Functions book, page 83. It's a port of la4j's RawLUCompositor.decompose() (in Java), cross-checked against Thomas Manz's ::matrix::dgetrf (TCL), which appears to trace back to LAPACK (Fortran).

Implementation

void decomposeLU(Matrix m) {
  assert(!m.isLU);
  m.isLU = true;
  for (int j = 0; j < m.columns; j++) {
    for (int i = 0; i < m.rows; i++) {
      int kMax = min(i, j);
      double s = 0;
      for (int k = 0; k < kMax; k++) {
        s += m.getF(i, k) * m.getF(k, j);
      }
      m.setF(i, j, m.getF(i, j) - s);
    }

    int pivot = j;

    for (int i = j + 1; i < m.rows; i++) {
      if (m.getF(i, j).abs() > m.getF(pivot, j).abs()) {
        pivot = i;
      }
    }
    if (pivot != j) {
      m.swapRowsLU(pivot, j);
    }
    if (j < m.rows) {
      final double vj = m.getF(j, j);
      if (vj.abs() > 0) {
        for (int i = j + 1; i < m.rows; i++) {
          m.setF(i, j, m.getF(i, j) / vj);
        }
      }
    }
  }

  // Avoid a singular matrix by perturbing the pivots, if needed, so they fall
  // within the 15C's precision.  See Advanced Functions, 98-99.
  double maxPivot = 0;
  for (int i = 0; i < m.rows; i++) {
    maxPivot = max(maxPivot, m.getF(i, i).abs());
  }
  final int minExp = max(-99, Value.fromDouble(maxPivot).exponent - 10);
  for (int i = 0; i < m.rows; i++) {
    final v = m.get(i, i);
    if (v.exponent < minExp) {
      if (v.asDouble < 0) {
        m.setF(i, i, -pow(10.0, minExp).toDouble());
      } else {
        m.setF(i, i, pow(10.0, minExp).toDouble());
      }
    }
  }
}