LU constructor

LU(Array2d A)

LU Decomposition Structure to access L, U and piv.

• A Rectangular matrix

Implementation

LU(Array2d A) {
  // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
  _LUMatrix = A.copy();
  _m = A.row;
  _n = A.column;
  _piv = Array.fixed(_m);
  for (var i = 0; i < _m; i++) {
    _piv[i] = i.toDouble();
  }
  _pivsign = 1;
  Array LUrowi;
  var LUcolj = Array.fixed(_m);

  // Outer loop.

  for (var j = 0; j < _n; j++) {
    // Make a copy of the j-th column to localize references.

    for (var i = 0; i < _m; i++) {
      LUcolj[i] = _LUMatrix[i][j];
    }

    // Apply previous transformations.

    for (var i = 0; i < _m; i++) {
      LUrowi = _LUMatrix[i];

      // Most of the time is spent in the following dot product.

      var kmax = math.min(i, j);
      var s = 0.0;
      for (var k = 0; k < kmax; k++) {
        s += LUrowi[k] * LUcolj[k];
      }

      LUrowi[j] = LUcolj[i] -= s;
    }

    // Find pivot and exchange if necessary.

    var p = j;
    for (var i = j + 1; i < _m; i++) {
      if (LUcolj[i].abs() > LUcolj[p].abs()) {
        p = i;
      }
    }
    if (p != j) {
      for (var k = 0; k < _n; k++) {
        var t = _LUMatrix[p][k];
        _LUMatrix[p][k] = _LUMatrix[j][k];
        _LUMatrix[j][k] = t;
      }
      var k = _piv[p].toInt();
      _piv[p] = _piv[j];
      _piv[j] = k.toDouble();
      _pivsign = -_pivsign;
    }

    // Compute multipliers.
    if (j < _m && _LUMatrix[j][j] != 0.0) {
      for (var i = j + 1; i < _m; i++) {
        _LUMatrix[i][j] /= _LUMatrix[j][j];
      }
    }
  }
}