code: @Override
      public void solve(DMatrixRMaj B, DMatrixRMaj X) {
          if( B.numRows != numRows) {
              throw new IllegalArgumentException("Unexpected matrix size");
          }
          X.reshape(numCols,B.numCols);
  
          if( B != X ) {
              SpecializedOps_DDRM.copyChangeRow(pivot,B,X);
          } else {
              throw new IllegalArgumentException("Current doesn't support using the same matrix instance");
          }
  
          // Copy right hand side with pivoting
          int nx = B.numCols;
          double[] dataX = X.data;
  
          // Solve L*Y = B(piv,:)
          for (int k = 0; k < numCols; k++) {
              for (int i = k+1; i < numCols; i++) {
                  for (int j = 0; j < nx; j++) {
                      dataX[i*nx+j] -= dataX[k*nx+j]*dataLU[i* numCols +k];
                  }
              }
          }
          // Solve U*X = Y;
          for (int k = numCols -1; k >= 0; k--) {
              for (int j = 0; j < nx; j++) {
                  dataX[k*nx+j] /= dataLU[k* numCols +k];
              }
              for (int i = 0; i < k; i++) {
                  for (int j = 0; j < nx; j++) {
                      dataX[i*nx+j] -= dataX[k*nx+j]*dataLU[i* numCols +k];
                  }
              }
          }
      }
docs: An other implementation of solve() that processes the matrices in a different order.
  It seems to have the same runtime performance as {@link #solve} and is more complicated.
  It is being kept around to avoid future replication of work.
  
  @param B A matrix that is n by m.  Not modified.
  @param X An n by m matrix where the solution is writen to.  Modified.