@Rouben Rostamian  Many thanks for the quick answer. 

I need further help to understand how the package works.  In the example below (pseudocode) I did all the calculations by going back or forward to avoid the time dependence and I wonder whether the Physics package will do the same elegantly. 

f:=[f1(x(t),y(t),f2(x(t),y(t)];
sys:={diff(x,t)=f[1],diff(y,t)=f[2]);
newsys:=subs({x=phi[1],y=phi[2]},sys);  # A transformation of Coordinates - Can we use the Physics package here?
# And another transformation of coordinates
z__1:=x;
z__2   # Lf(z__1) - Lie derivative of z__1 in f
z__3  # Lf(z__2) 
.....

Many thanks

@Carl Love and @mmcdara Many thanks.   

@Carl Love I remember reading in one of your answers only one of the time outputs of Usage matters when using Grid or Threads but I don't remember which one.  Could you remind me?

I did not have time to run a memory demanding example but a small (not that small) example showed me that Carl's solution is blazing fast.   

@Carl Love Since your solution is very concise, could you be so kind to open it up to show the details (curry, splitranks,..)?  

@Carl Love Many thanks.   Before reading your reply I tried exactly what you suggested (first option) and it worked.   

@tomleslie Many thanks.   Grid:-Set(UtilsIORelation) did the trick for removing the error message.  Now I have to check how to deal with the other issues.  

@Carl Love Many thanks.  I have just sent you an email with a link to the file with only two variables (parms10 and parms11).  

The Grid version of my procedure is running so far without any problem.   I hope Grid does not choke (as it does in other calculations for no reason).  

@Carl Love Many thanks.  Yes, "kernel panic" meant "maple: fatal error, lost connection to kernel". 

By changing the line, I got a different message - "Execution stopped: Stack limit reached."  - The limit was changed but the same msg came out.  kernelopts adjustments made no difference.  

@Carl Love Many thanks.  I am using Maple 2020.2 on both machines:  the server and the MacBook pro.  I sent the examples from the MacBook pro because I was running a much larger example on the server using your code without any modification (my attempt to modify your not-threaded code is not worth spending time). The outcome was "kernel panic". 

parse/sprintf were used to select the examples without having to modify the worksheet (I use batch files to run mpl files). As you know I had to save the intermediate results in files due to the size of the computations.  

As I said the use of the threaded version ended up with a kernel panic on the server.   

ProbSize: [[1122436, 1639162], [11224360, 17995320], [29, 29]]
 

@Carl Love Many thanks.  Here are the results for a less computationally demanding case.    SubsetPairsSparse2 is a procedure with the last line changed to [Threads:-Seq](`intersect`(seq(T[x], x= a)),a in A).

Questions:

1) Why does Maple return numcpus when kernelopts is called for the first time?

2) Why does Maple return 2 instead of 4 in some of the subsequent calls of kernelopts?  

@Carl Love The result of issuing the command shown above and running SubsetPairsSparse is

memory used=16.37GiB, alloc change=0 bytes, cpu time=5.53h, real time=25.90m, gc time=5.48h
 

There is certainly an improvement over the last calculation but did not reach 17 minutes when using Grid (as I said Grid is not a good option because I don't know when it is going to hang for no reason and memory demand).  

@vv I am not sure what you meant by all computations. No, it is no a CPU test, although if I use Grid:-Map, it shows that there is a problem when Grid is used for large computations (reported in another post to MaplePrimes).   Yes, I could use C(or Fortran, sorry but I am old school) but I am trying to have everything under one platform, that is, Maple.  I used to have parts in Mathematica, Matlab, R, and Fortran.   

@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.