The other solution, involving the replacement x2->x1 at the end, should work but it's rather cumbersome for the applications I have in mind, where the number and entity of repeated terms varies.

That's easy to handle

In[12]:= Clear[deNicolaPairs]
deNicolaPairs[l_List?VectorQ] := Block[{l0 = Array[f, Length[l]]},
   Select[Subsets[Subsets[l0, {2}], {Length[l]/2}], Sort[Flatten[#]] == l0 &] /. (Rule @@@ Transpose[{l0, l}])
   ] /; EvenQ[Length[l]] && DisjointQ[l, Array[f, Length[l]]]

In[14]:= deNicolaPairs[{x1, x2, x3, x4}]
Out[14]= {{{x1, x2}, {x3, x4}}, {{x1, x3}, {x2, x4}}, {{x1, x4}, {x2, x3}}}

In[15]:= deNicolaPairs[{x1, x1, x1, x1}]
Out[15]= {{{x1, x1}, {x1, x1}}, {{x1, x1}, {x1, x1}}, {{x1, x1}, {x1, x1}}}

but the performance problem is fatal: deNicolaPairs crashes with 8 GB RAM on a list of length 12.

The proposal in the mentioned forum using ArrayReshape fixed that by generating first an index (it depends on the length of the input only) and then distributing the input accordingly.

This doesn't quite do what I mean: for example, I changed it for a case with 6 variables (one is repeated, which is the case in my actual problem): {x1, x1, x3, x4, x5, x6}

Select[Subsets[Subsets[{x1, x1, x3, x4, x5, x6}, {2}], {3}], 
 Sort[Flatten[#]] == {x1, x1, x3, x4, x5, x6} &]

The output has 27 elements, while it should have 15 (5x3).

Perhaps I haven't formulated the problem clearly... What I'm looking for is all the possible partitions of the initial list in pairs of elements, i.e. a list of such partitions, as given in my example, rather than just the subsets. e.g. the partitions of {x1,x2,x3,x4} are {{x1,x2},{x3,x4}}, {{x1,x3},{x2,x4}} and {{x1,x4},{x2,x3}} (the order does not matter), so the output I want is: {{{x1,x2},{x3,x4}}, {{x1,x3},{x2,x4}}, {{x1,x4},{x2,x3}}}.

@Tim Mayes, thank you but what I'm looking for is not just finding all possible pairs but find all possible partitions of the list into pairs, as in the example I gave