docs: @Override
      public void solve(DMatrixRMaj B, DMatrixRMaj X) {
          if( B.numRows != numRows )
              throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numRows);
          X.reshape(numCols,B.numCols);
  
          int BnumCols = B.numCols;
          
          // solve each column one by one
          for( int colB = 0; colB < BnumCols; colB++ ) {
  
              // make a copy of this column in the vector
              for( int i = 0; i < numRows; i++ ) {
                  a.data[i] = B.data[i*BnumCols + colB];
              }
  
              // Solve Qa=b
              // a = Q'b
              // a = Q_{n-1}...Q_2*Q_1*b
              //
              // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
              for( int n = 0; n < numCols; n++ ) {
                  double []u = QR[n];
  
                  double vv = u[n];
                  u[n] = 1;
                  QrHelperFunctions_DDRM.rank1UpdateMultR(a, u, gammas[n], 0, n, numRows, temp.data);
                  u[n] = vv;
              }
  
              // solve for Rx = b using the standard upper triangular solver
              TriangularSolver_DDRM.solveU(R.data,a.data,numCols);
  
              // save the results
              for( int i = 0; i < numCols; i++ ) {
                  X.data[i*X.numCols+colB] = a.data[i];
              }
          }
      }