If matrix $A$ comes in Gram-factored form $A=X^T X$, then pivoted Cholesky can be done through QR like in this answer. QRDecomposition has a Pivoted->True option.
(* given X produces R,p such that p\[Transpose].X\[Transpose].X.p==R\
\[Transpose].R *)
PivotedCholeskyDecompositionOfXTX[X_] := Module[{q, R, p},
{q, R, p} = QRDecomposition[X, Pivoting -> True];
{Map[If[Negative[#], -1, 1] &, Diagonal[R]]*R, p}
];
(* from https://mathematica.stackexchange.com/a/271404/217 *)
CholeskyDecompositionOfXTX[X_] :=
With[{R = Last[QRDecomposition[X]]},
Map[If[Negative[#], -1, 1] &, Diagonal[R]]*R];
d = 5;
numSamples = 5;
evals = Table[1/i, {i, 1, d}];
sigma = DiagonalMatrix[evals];
dist = MultinormalDistribution[sigma];
X = RandomVariate[dist, numSamples];
R = CholeskyDecompositionOfXTX[X];
Print["Unpivoted diagonal is decreasing: ",
Diagonal[R] == ReverseSort@Diagonal@R] (* False *)
{R, p} = PivotedCholeskyDecompositionOfXTX[X];
Print["Pivoted reconstruction works: ",
Norm[p\[Transpose] . X\[Transpose] . X . p - R\[Transpose] . R] <
10^-10]; (* True *)
Print["Pivoted diagonal is decreasing: ",
Diagonal[R] == ReverseSort@Diagonal@R] (* True *)
If it's not factored, then one approach is this (needs Mathematica implementation, one issue is that QRDecomposition gives reduced decomposition, but that approach needs the full version)