luDecomposition method

@override List<RealMatrix> luDecomposition()

Factors the matrix as the product of a lower triangular matrix L and an upper triangular matrix U. The matrix must be square.

The returned list contains L at index 0 and U at index 1.

A MatrixException object is thrown if the matrix isn't square.

Implementation

@override
List<RealMatrix> luDecomposition() {
  // Making sure that the matrix is squared
  if (!isSquareMatrix) {
    throw const MatrixException(
      'LU decomposition only works with square matrices!',
    );
  }

  // Creating L and U matrices
  final L = List<double>.generate(
    rowCount * columnCount,
    (_) => 0.0,
    growable: false,
  );
  final U = List<double>.generate(
    rowCount * columnCount,
    (_) => 0.0,
    growable: false,
  );

  // Computing L and U
  for (var i = 0; i < rowCount; ++i) {
    for (var k = i; k < rowCount; k++) {
      // Summation of L(i, j) * U(j, k)
      var sum = 0.0;
      for (var j = 0; j < i; j++) {
        sum += _getDataAt(L, i, j) * _getDataAt(U, j, k);
      }

      // Evaluating U(i, k)
      _setDataAt(U, i, k, this(i, k) - sum);
    }

    // Lower Triangular
    for (var k = i; k < rowCount; k++) {
      if (i == k) {
        _setDataAt(L, i, i, 1);
      } else {
        // Summation of L(k, j) * U(j, i)
        var sum = 0.0;
        for (var j = 0; j < i; j++) {
          sum += _getDataAt(L, k, j) * _getDataAt(U, j, i);
        }

        // Evaluating L(k, i)
        _setDataAt(L, k, i, (this(k, i) - sum) / _getDataAt(U, i, i));
      }
    }
  }

  return [
    RealMatrix.fromFlattenedData(
      rows: rowCount,
      columns: rowCount,
      data: L,
    ),
    RealMatrix.fromFlattenedData(
      rows: rowCount,
      columns: rowCount,
      data: U,
    ),
  ];
}