src/math-helpers.js
function luSolve (A, P, b, X) {
  // Since PA = LU, then L(U x) = Pb
  const a = A.elements;
  const n = A.dimension;
  // L * y  = P * b, solve for y.
  // Implicit 1's on the diagonal.
  for (let i = 0; i < n; ++i) {
    let sum = 0;
    for (let j = 0; j < i; ++j) {
      sum += a[i + j * n] * X.getComponent(j);
    }
    const xi = b.getComponent(P[i]) - sum;
    X.setComponent(i, xi);
  }

  // U * x = y  ==> i = n-1 -> 0, j = 0
  for (let i = n - 1; i >= 0; --i) {
    let sum = 0;
    for (let j = i + 1; j < n; ++j) {
      sum += a[i + j * n] * X.getComponent(j);
    }
    const scale = a[i + i * n];
    // if (scale === 0) {
    //   consoleWarning(`luSolve(): x[${i}] is free.`);
    // }
    X.setComponent(i, (X.getComponent(i) - sum) / scale);
  }
  return X;
}
function luSolve (A, P, b, X) {
  // Since PA = LU, then L(U x) = Pb
  const a = A.elements;
  const n = A.dimension;
  // L * y  = P * b, solve for y.
  // Implicit 1's on the diagonal.
  for (let i = 0; i < n; ++i) {
    let sum = 0;
    for (let j = 0; j < i; ++j) {
      sum += a[i + j * n] * X.getComponent(j);
    }
    const xi = b.getComponent(P[i]) - sum;
    X.setComponent(i, xi);
  }

  // U * x = y  ==> i = n-1 -> 0, j = 0
  for (let i = n - 1; i >= 0; --i) {
    let sum = 0;
    for (let j = i + 1; j < n; ++j) {
      sum += a[i + j * n] * X.getComponent(j);
    }
    const scale = a[i + i * n];
    // if (scale === 0) {
    //   consoleWarning(`luSolve(): x[${i}] is free.`);
    // }
    X.setComponent(i, (X.getComponent(i) - sum) / scale);
  }
  return X;
}
function luSolve ( A , P , b , X ) { // Since PA = LU, then L(U x) = Pb const a = A . elements ; const n = A . dimension ; // L * y  = P * b, solve for y. // Implicit 1's on the diagonal. for ( let i = 0 ; i < n ; ++ i ) { let sum = 0 ; for ( let j = 0 ; j < i ; ++ j ) { sum += a [ i + j * n ] * X . getComponent ( j ) ; } const xi = b . getComponent ( P [ i ] ) - sum ; X . setComponent ( i , xi ) ; } // U * x = y  ==> i = n-1 -> 0, j = 0 for ( let i = n - 1 ; i >= 0 ; -- i ) { let sum = 0 ; for ( let j = i + 1 ; j < n ; ++ j ) { sum += a [ i + j * n ] * X . getComponent ( j ) ; } const scale = a [ i + i * n ] ; // if (scale === 0) { //   consoleWarning(`luSolve(): x[${i}] is free.`); // } X . setComponent ( i , ( X . getComponent ( i ) - sum ) / scale ) ; } return X ; }
Solves the Ax=b problem with an LU decomposition.
@memberof Helpers
@param {Matrix2|Matrix3|Matrix4} A The LU decomposed matrix.
@param {number[]} P The permutation array.
@param {Vector2|Vector3|Vector4} b The right side vector in the Ax=b problem.
@param {Vector2|Vector3|Vector4} X The solution to the problem.
Solves the Ax = b problem with an LU decomposition .