(* Mathematica package *)

(* aux *)
Clear[xName, pName]
pName[i_, j_] :=
    ToExpression[StringJoin["p", ToString[i], ToString[j]]];
xName[i_] :=
    ToExpression[StringJoin["x", ToString[i]]];
Clear[pObjective, pConstraints]
pObjective[mx_, xs_] :=
    Range[Length[xs]].mx.xs;
pConstraints[mx_] :=
    Module[ {rowC, colC, n},
        n = Length[mx];
        rowC = Flatten /@ (mx . Table[1, {n}]);
        colC = Flatten /@ (Transpose[mx].Table[1, {n}]);
        Join[rowC, colC,-rowC,-colC]
    ];
Clear[progM, progB, progC]
progM[pCon_, mx_] :=
    Outer[Coefficient, pCon, Flatten[mx]];
progB[pCon_, mx_] :=
    Join[Table[1, {Length[pCon]/2}],Table[-1, {Length[pCon]/2}]];
progC[pObj_, mx_] :=
    First[progM[{pObj}, mx]];

(* main *)
lpSort[xIn_,givePerm_:False] :=
    Module[ {n,mx,xs,pObj,pCon,m,b,c,lpP,sorted},
        n = Length[xIn];
        mx  = Table[pName[i, j], {i, n}, {j, n}];
        xs = Table[xName[i], {i, n}];
        pObj = pObjective[mx, xs];
        pCon = pConstraints[mx];
        m = progM[pCon, mx];
        b = progB[pCon, mx];
        c = progC[pObj, mx];
        c = (c /. Thread[xs -> xIn]);
        lpP = Partition[LinearProgramming[-c, -m, -b], Length[xs]];
        sorted = lpP . xIn;
        If[ givePerm,
            {sorted,lpP},
            sorted
        ]
    ]