QR constructor

Implementation

QR(Array2d A) {
  // Initialize.
  _QR = A.copy();
  _m = A.row;
  _n = A.column;
  _Rdiag = Array.fixed(_n);

  // Main loop.
  for (var k = 0; k < _n; k++) {
    // Compute 2-norm of k-th column without under/overflow.
    var nrm = 0.0;
    for (var i = k; i < _m; i++) {
      nrm = hypotenuse(nrm, _QR[i][k]);
    }

    // A matrix wider than tall has no diagonal element in this column, so
    // R has a zero there and the matrix is not full rank.
    if (k >= _m) {
      _Rdiag[k] = 0.0;
      continue;
    }

    // As in LAPACK, no reflection is applied when there is nothing below
    // the diagonal: R keeps the original sign of the diagonal element.
    var belowDiagonalIsZero = true;
    for (var i = k + 1; i < _m; i++) {
      if (_QR[i][k] != 0.0) {
        belowDiagonalIsZero = false;
        break;
      }
    }
    if (belowDiagonalIsZero) {
      _Rdiag[k] = _QR[k][k];
      _QR[k][k] = 0.0;
      continue;
    }

    if (nrm != 0.0) {
      // Form k-th Householder vector. The sign of nrm follows the diagonal
      // element to avoid cancellation in _QR[k][k] += 1.0 below.
      if (_QR[k][k] < 0) {
        nrm = -nrm;
      }
      for (var i = k; i < _m; i++) {
        _QR[i][k] /= nrm;
      }
      _QR[k][k] += 1.0;

      // Apply transformation to remaining columns.
      for (var j = k + 1; j < _n; j++) {
        var s = 0.0;
        for (var i = k; i < _m; i++) {
          s += _QR[i][k] * _QR[i][j];
        }
        s = -s / _QR[k][k];
        for (var i = k; i < _m; i++) {
          _QR[i][j] += s * _QR[i][k];
        }
      }
    }
    _Rdiag[k] = -nrm;
  }
}