(* 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 ] ]