So I have a rather involved question here (and some ideas of my own about how to do this):
What is an efficient method to find order 16 subgroups of a given finite group?
Obviously there are many ways to do this, including brute force computation, but I am mainly interested in efficient methods, that exploit easily computed information. Note that I already have plenty of information on groups of order 2, 4, 8, and 16 to draw from and help me in determining the isomorphism type of the subgroups in question. So I am basically looking for a systematic procedure that could be used to find and identify the order 16 subgroups of a given group, exploiting knowledge of the order 16 groups and easily computed information from the given finite group.
As an example of what I mean, here is a pretty effective method for order $8$ subgroups:
-Any order 8 elements will generate a cyclic subgroup(s) of order 8, each one will have four order 8 elements (which can only be in that one cyclic order 8 subgroup).
-To find copies of $\mathbb{Z}_2\times\mathbb{Z}_4$, look for pairs of order 4 cyclic subgroups with common order 2 element, take one order 4 element from each and compute the product $a*b=c$ and $b*a$, if $b*a=c$ and if that element is order 2, compute $a^2*c$ and $c*a^2$ (note that $a^2=b^2$), if this shows that $a^2$ and $c$ commute, then the element $a^2*c$ is also order 2 and forms the third element of a subgroup isomorphic to $\mathbb{Z}_2\times\mathbb{Z}_2$, and you have found a copy of $\mathbb{Z}_2\times\mathbb{Z}_4$ (made up of the order 4 cycles and the klein four subgroup you found). In doing the above, if any of the conditions were violated, you can end the process immediately knowing that the given order 4 cycles are not part of a subgroup isomorphic to $\mathbb{Z}_2\times\mathbb{Z}_4$.
-To find copies of the dihedral group $D_8$, look for order 4 cycles whose order 2 element is contained in at least 2 separate subgroups isomorphic to $\mathbb{Z}_2\times\mathbb{Z}_2$ (these can be found by multiplying pairs of order 2 elements, in both orders, to see if they generate the same order 2 element). Test the unshared order 2 elements of these subgroups to see if they fail to commute (and generate the either order 4 element from the given order 2 cycle). If so, you have found a copy of $D_8$, if anything fails along the way, you know you haven't, (note that if there are enough klein four groups with the that share that same order 2 element, then you may have quite a few pairs to check before concluding that the given order 4 cycle is not part of any copy of $D_8$).
-To find copies of the quaternion group, $Q_8$, look for trios of 4 cycles that share a single order 2 element, choose a generator for each cycle and check that it does not commute with the other generators and that the products are always one of the order 4 elements from the trio. If so, you have found a copy of $Q_8$, otherwise, check a different trio (if more than 3 share a common order 2 element, then you can rule out any pair of order 4 gens that commute immediately from being in the same $Q_8$ copy).
-And finally for copies of the elementary abelian group, $\mathbb{Z}_2\times\mathbb{Z}_2\times\mathbb{Z}_2$, look for at least 7 order 2 elements, if there are, then try to find 2 that commute, they will generate a 3rd and form a klein 4 subgroup, next check order 2 elements no in the klein 4 subgroup again 2 of the 3 that are. If it commutes with both, you have found a copy of $\mathbb{Z}_2\times\mathbb{Z}_2\times\mathbb{Z}_2$ (if there are lots of order 2 elements, be sure to compute which 7 are in this subgroup).
I know something about the way Magma computes subgroups of finite groups. I expect similar considerations apply to GAP, but I know less about that. In Magma, you can ask for subgroups of a group $G$ of a specific size. For example:
This actually returns representatives of the conjugacy classes of subgroups of order $16$ in the group $G$ rather than all such subgroups.
The algorithm used does not know anything about the structure of the individual groups of order $16$, so it might not be what you are looking for. But it is not clear from your post whether your interest is just theoretical or whether you really want to exceute this algorithm, and I would be surprised if using information about the groups of order $16$ would make much difference to its performance in practice.
A quick summary of the algorithm is that if first finds a series
$$1 = N_0 < N_1 < \cdots N_r = L \le G$$
of normal subgroups of $G$, where $L$ is the largest normal solvable subgroup of $G$ and each $L_i/L_{i-1}$ is elementary abelian. It uses stored information about the subgroups of the simple groups to find the subgroups of $G/L$ (so it is not a general purpose algorithm - if the group has a nonabelian composition factor that is too large then it will give up, but this is not really a serious drawback because any algorithm of this kind will fail in practice if the group is too big).
Then if lifts the subgroups of $G/L$ through the elementary abelian layers $L_i/L_{i-1}$ in turn eventaul finding the subgroups of $G$. If you ask for subgroups of a specific order, then it will only lift subgroups of $G/L$ with order dividing the specified order.