Consider some dataset

    N1 = 13;
    N2 = 130*10^3;
    N3 = 2*10^4;
    tabtest = 
      Join[Join[RandomReal[{0, 3}, {N1, 4}], Table[{1.}, N1], 2], 
       Join[RandomReal[{0, 3}, {N2, 4}], Table[{2.}, N2], 2], 
       Join[RandomReal[{0, 3}, {N3, 4}], Table[{7.}, N3], 2]];

I need to get a list of lists with the same last element. The obvious one is 

    SplitBy[tabtest,Last];

However, it is relatively slow, taking 0.15 s on my machine. I made another code:

    PhaseSpaceSplitter[phasespace_] := 
     Module[{grouped, lengths, positions, ranges, pso, pdgs},
      pdgs = phasespace[[All, -1]];
      If[Length[Union[pdgs]] != 1,
       (*Grouping consecutive identical elements*)
       grouped = Split[pdgs];
       (*Calculating the lengths of each group*)
       lengths = Length /@ grouped;
       (*Generating the start positions*)
       positions = Accumulate[Prepend[lengths, 0]];
       (*Forming ranges*)
       ranges = 
        Transpose[{Most[positions] + 1, Most[positions] + lengths}];
       pso = Take[phasespace, {#[[1]], #[[2]]}] & /@ ranges
       ,
       pso = {phasespace};
       ];
      pso
      ]

    splitted1 = SplitBy[tabtest, Last]; // AbsoluteTiming//First
    splitted2 = PhaseSpaceSplitter[tabtest]; // AbsoluteTiming//First
    splitted1 == splitted2

> 0.188216

> 0.0133451 

> True

I am okay with how fast the last approach works. However, the speed gain is lost once I apply `RandomSample` on tabtest:

    tabtest=tabtest//RandomSample;
    splitted1 = SplitBy[tabtest, Last]; // AbsoluteTiming//First
    splitted2 = PhaseSpaceSplitter[tabtest]; // AbsoluteTiming//First
    splitted1 == splitted2

> 0.152886
> 
> 0.243063
> 
> True

How to make an efficient code that makes the list of lists?

P.S. I cannot first make the backward sorting `SortBy[tabtest,#[[1]]&]` since `tabtest` is a set of columns of another table, and the ordering is important.