docs: private double[][] solveInplace(double[][] B) {
      int rows = B.length, nx = B[0].length;
      if(rows != 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++) {
        for(int j = 0; j < nx; j++) {
          double s = 0.0;
          for(int i = k; i < m; i++) {
            s += QR[i][k] * B[i][j];
          }
          s /= QR[k][k];
          for(int i = k; i < m; i++) {
            B[i][j] -= s * QR[i][k];
          }
        }
      }
      // Solve R*X = Y;
      for(int k = n - 1; k >= 0; k--) {
        final double[] Xk = B[k];
        for(int j = 0; j < nx; j++) {
          Xk[j] /= Rdiag[k];
        }
        for(int i = 0; i < k; i++) {
          final double[] Xi = B[i];
          final double QRik = QR[i][k];
          for(int j = 0; j < nx; j++) {
            Xi[j] -= Xk[j] * QRik;
          }
        }
      }
      return B;
    }