@Carl Love Many thanks.  Sorry, I meant nops (old matlab habit). 

For the largest example (not discussed here), the output of Probsize is: [[1905271, 1891223], [34294878, 32150791], [30, 30]]

For the example shown in Parallelization: [[52190, 121802], [365330, 974416], [29, 29]]

I run SubsetPairsSparse with the same example and obtained the following metrics:

memory used=16.37GiB, alloc change=66.22MiB, cpu time=7.42h, real time=33.65m, gc time=7.35

approximately 7 times faster than SubsetPairs1a. If you can, please send the threaded version of it.

Since I am not sure if the procedure will work for larger problems (three trillion - a bunch of yikes), I was thinking that perhaps the following problem is feasible: given lists of sets A and B return a list L of sets in B that does not contain any set of A. 

@Carl Love @vv  @tomleslie Thank you for the procedures.  I am testing them with L1 (length=52190) and L2(lenth=121802).   Some of the procedures are taking more than 3 hours to finish.   I have modified my procedure to use Grid:-Map (Unfortunately, I don't know how to use infolevel[Threads:-Map] to check the progress, and therefore I left the use of Threads for later), and by doing that, the time is down to 17 minutes.   I guess all the procedures will need some tuning to use parallelization tools.   I can't imagine how long it will take to run the procedure when the lengths are over one million.   

@Carl Love Many, many thanks for the implementation of the bipartite graph.  I would appreciate if you could show me how to use embellishment options.  

Updates:

a) @tomleslie 's second routine

memory used=170.20GiB, alloc change=4.80MiB, cpu time=77.26h, real time=3.04h, gc time=76.90h
 

b) @Carl Love 's procedure 

memory used=0.65TiB, alloc change=0 bytes, cpu time=20.05h, real time=3.26h, gc time=17.38h
 

c) @vv 's procedure

memory used=0.74TiB, alloc change=23.93GiB, cpu time=57.76h, real time=5.09h, gc time=54.58h
 

Bipartite option returns "Error, (in GraphTheory:-DrawGraph) input graph is not bipartite" using the example shown in the post above.  

@Carl Love what I meant was to check if the calculations could somehow demand less time and memory when the replacement is used.  That is definitely the case. Many many thanks. 

@Carl Love Many thanks.  I am doing some tests right now with the replacement for abs. I will let you know how it goes. 

If @ is what you suspect it is, I will have to learn to use it efficiently.  

I agree with you.  As far as I could see, the solutions come out with the same general form.   

The initial equations  (the ones that will go through SolveTools:-PolynomialSystem to get the solutions you see in the examples in this post) are polynomial functions with rational coefficients. I hope this is what you are referring to.  

@vv What shown here is just part of the problem.  There is an expression in which y and z are replaced by the solutions you see in the posts.  I have noticed that, if there is no difference in the expression when I replaced y by -y and z by -z,   it makes no difference to use either a root or its symmetric counterpart.  The number of such systems is very small, but unfortunately, I came across some of them. 

@Carl Love You are right.  abs seems to consume all the available memory in some cases, bringing the server almost to a halt. Yes, the idea is to put e and -e in the same equivalence class (the symmetric roots included).   Unfortunately, there are some radicals involved in the calculations as well as polynomial and rational functions.  I thought I would not need to be worried about the solutions with radicals but some simple models proved me wrong.   Below is an example:

ans := {{y = sqrt(alpha[1, 7]*(X2 - alpha[1, 0]))/alpha[1, 7], z = -(2*X1*X2*alpha[2, 5] - 2*X1*alpha[1, 0]*alpha[2, 5] - X3)/(2*alpha[2, 3]*sqrt(alpha[1, 7]*(X2 - alpha[1, 0])))}, {y = -sqrt(alpha[1, 7]*(X2 - alpha[1, 0]))/alpha[1, 7], z = (2*X1*X2*alpha[2, 5] - 2*X1*alpha[1, 0]*alpha[2, 5] - X3)/(2*alpha[2, 3]*sqrt(alpha[1, 7]*(X2 - alpha[1, 0])))}}

@Carl Love Thanks.  Last question - I modified  RemoveSymmetrySols to accommodate another condition, that is, class[abs~(rhs~(e))][e] and it works as intended.  To compare it to the previous solution I need to add something to LIstTools:-Classify to do the same. Although I had it modified, I confess I am not sure if my solution is efficient as far as speed is concerned.   Can you help me again?   Many thanks.

op~(1, {entries}(ListTools:-Classify(w->abs~(rhs~(w)), Sol), 'nolist'));

@Carl Love Many thanks.  What I did not get was the use and the reason for ':-nolist' (especially the apostrophe and :~).

The use of table was very clever. I need to incorporate that into my procedures.  

@Carl Love Many thanks.  As soon as I get the results of a job that is running right now and I will have something to compare both procedures.  
 

Would you be so kind to explain what the line with op does? Specially the reason of the use of '. 
 

@Carl Love Sorry for using here as a reply to three earlier posts.  1)  Thanks for pointing out the use [],   2) I will definitely take into account what you said about the help pages and 3)  many thanks.   

Is ListTools:-Classify optimized for speed?   You know that I will have to use it for a large number of system of equations

@Carl Love Many thanks.  I am confused with what the help page says and that is the reason why I posted the question.  What does "partial solutions" mean?   If I add a notzero as {x,y} all solutions with a zero (x and y) are eliminated.  Does it mean that {x=0.y=0} is considered as partial solution?  

As for your solution I had to remove the last [] before it works.  

I know if I use set as you did if there are identical solutions (not in the case shown above) only one is kept. Can I do something similar with symmetric solutions?  That is, only one is kept.  

@vv Thank you for the answer. I am not sure if the answer will help test millions of equations but it surely helps me to understand some of the possibilities (without using assuming).  

@Carl Love Many thanks.  I will give a try on simplify(..., symbolic) and also check assuming real. 

Acer's answer helps me with another problem - part b) of my original post.  

@acer Thanks.  The methods are not top-secret at all: 1st Method - solve the general equation with unknown symbolic coefficients for the monomials (beta, for instance).  Then, by using the relation between beta coefficients and alpha coefficients (such as the ones you see in my first post), find the solution.  2nd Method -solve the system with original alpha coefficients (they have an intricate relationship that depends on other calculations).  
My expectation was to reduce the number of systems to be solved using the general equation.  Of course, I have no choice but to use Method 2 in some cases (again due to the relationships in the alpha coefficients). 

I hope this sheds some light on the problem.  

@acer Many thanks.   I will avoid evalb in future calculations.  Since I am away from my notebook, I wonder whether assuming is the key to all alternatives.  If so, that will make the comparison almost impossible since I have millions of such a system to be tested and I can't see how I will write a specific assuming for each one of them.