docs: function luDecomposition (A, inPlace) {
    const n = A.dimension;
    const AA = inPlace ? A : A.clone();
    const a = AA.elements; // indexed via a_ij = a[i + j * n]
    const P = [];
  
    const rowScalers = [];
    for (let i = 0; i < n; ++i) {
      P.push(i);
      const col = findLargestInRow(A, i);
      const scaler = a[col * n + i]; // row scaling
      if (scaler === 0) {
        throw new Error('luDecomposition(): matrix is singular and cannot be LU factorized.');
      }
      rowScalers[i] = scaler; // don't actually want to scale the matrix or else (PA != LU).
    }
    // iterate over columns to reduce the matrix
    // implicitly reduce the matrix
    for (let j = 0; j < n; ++j) {
      for (let i = 0; i < j; ++i) {
        // compute upper tri, computes values across the row
        let sum = 0;
        for (let k = 0; k < i; ++k) {
          sum += a[i + k * n] * a[k + j * n]; // avoid big + small roundoffs
        }
        a[i + j * n] = a[i + j * n] - sum;
      }
      let pivotLrgElem = 0;
      let pivotIndex = j;
      for (let i = j; i < n; ++i) {
        // compute lower tri and diagonal, computes values down the column
        let sum = 0;
        for (let k = 0; k < j; ++k) {
          sum += a[i + k * n] * a[k + j * n]; // avoid big + small roundoffs
        }
        a[i + j * n] = a[i + j * n] - sum;
  
        // find the pivot element
        const pivotTest = rowScalers[i] * Math.abs(a[i + j * n]);
        if (pivotTest > pivotLrgElem) {
          pivotLrgElem = pivotTest;
          pivotIndex = i;
        }
      }
      swapRowsInPlace(A, j, pivotIndex);
      swapValuesInArray(P, j, pivotIndex);
      swapValuesInArray(rowScalers, j, pivotIndex);
      if (j < n - 1) {
        const pivotScale = a[j + j * n];
        for (let i = j + 1; i < n; ++i) {
          a[i + j * n] /= pivotScale;
        }
      }
    }
    return {
      P,
      A
    };
  }