docs: @Override
      public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
          if( X.numRows != numCols )
              throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
          else if( B.numRows != numRows || B.numCols != X.numCols )
              throw new IllegalArgumentException("Unexpected dimensions for B");
  
          U = decomposer.getR(U,true);
          final double gammas[] = decomposer.getGammas();
          final double dataQR[] = QR.data;
  
          final 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++ ) {
                  int indexB = (i*BnumCols + colB)*2;
                  a[i*2]   = B.data[indexB];
                  a[i*2+1] = B.data[indexB+1];
              }
  
              // Solve Qa=b
              // a = Q'b
              // a = Q_{n-1}...Q_2*Q_1*b
              //
              // Q_n*b = (I-gamma*u*u^H)*b = b - u*(gamma*U^H*b)
              for( int n = 0; n < numCols; n++ ) {
                  int indexU = (n*numRows + n + 1)*2;
  
                  double realUb = a[n*2];
                  double imagUb = a[n*2+1];
  
                  // U^H*b
                  for( int i = n+1; i < numRows; i++ ) {
                      double realU = dataQR[indexU++];
                      double imagU = -dataQR[indexU++];
  
                      double realB = a[i*2];
                      double imagB = a[i*2+1];
  
                      realUb += realU*realB - imagU*imagB;
                      imagUb += realU*imagB + imagU*realB;
                  }
  
                  // gamma*U^T*b
                  realUb *= gammas[n];
                  imagUb *= gammas[n];
  
                  a[n*2]   -= realUb;
                  a[n*2+1] -= imagUb;
  
                  indexU = (n*numRows + n + 1)*2;
                  for( int i = n+1; i < numRows; i++) {
                      double realU = dataQR[indexU++];
                      double imagU = dataQR[indexU++];
  
                      a[i*2]   -= realU*realUb - imagU*imagUb;
                      a[i*2+1] -= realU*imagUb + imagU*realUb;
                  }
              }
  
              // solve for Rx = b using the standard upper triangular solver
              TriangularSolver_ZDRM.solveU(U.data, a, numCols);
  
              // save the results
  
              for( int i = 0; i < numCols; i++ ) {
                  int indexX = (i*X.numCols+colB)*2;
  
                  X.data[indexX]   = a[i*2];
                  X.data[indexX+1] = a[i*2+1];
              }
          }
      }