docs: public double[] solveInplace(double[] b) {
      if(b.length != m) {
        throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
      }
      if(!this.isFullRank()) {
        throw new ArithmeticException(ERR_MATRIX_RANK_DEFICIENT);
      }
      // Compute Y = transpose(Q)*B
      for(int k = 0; k < n; k++) {
        double s = 0.0;
        for(int i = k; i < m; i++) {
          s += QR[i][k] * b[i];
        }
        s /= QR[k][k];
        for(int i = k; i < m; i++) {
          b[i] -= s * QR[i][k];
        }
      }
      // Solve R*X = Y;
      for(int k = n - 1; k >= 0; k--) {
        b[k] /= Rdiag[k];
        for(int i = 0; i < k; i++) {
          final double QRik = QR[i][k];
          b[i] -= b[k] * QRik;
        }
      }
      return b;
    }