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");
  
          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.data[i*2]   = B.data[indexB];
                  a.data[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^T)*b = b - u*(gamma*U^T*b)
              for( int n = 0; n < numCols; n++ ) {
                  double []u = QR[n];
  
                  double realVV = u[n*2];
                  double imagVV = u[n*2+1];
  
                  u[n*2]   = 1;
                  u[n*2+1] = 0;
  
                  QrHelperFunctions_ZDRM.rank1UpdateMultR(a, u,0, gammas[n], 0, n, numRows, temp.data);
  
                  u[n*2]   = realVV;
                  u[n*2+1] = imagVV;
              }
  
              // solve for Rx = b using the standard upper triangular solver
              TriangularSolver_ZDRM.solveU(R.data, a.data, numCols);
  
              // save the results
              for( int i = 0; i < numCols; i++ ) {
                  int indexB = (i*BnumCols + colB)*2;
                  X.data[indexB]   = a.data[i*2];
                  X.data[indexB+1] = a.data[i*2+1];
              }
          }
      }