Good question: I think the tutorial/finiteT/ancilla.cc code you referenced would be a good starting point for doing what you are asking, but for other site types.
That code has a portion like this:
auto psi = MPS(sites);
for(int n = 1; n <= 2*N; n += 2)
auto s1 = sites(n);
auto s2 = sites(n+1);
auto wf = ITensor(s1,s2);
psi.ref(n) = ITensor(s1);
psi.ref(n+1) = ITensor(s2);
psi.ref(n) *= D;
The main parts that would need to be generalized are:
instead of defining a singlet state |up,dn> - |dn,up> as in the above code, you'll probably want to define some other kind of maximally entangled state. As you probably know, the ancilla/purification method doesn't depend on the precise initial state used, as long as it's maximally entangled (unitarily equivalent to a Bell-pair state, meaning up to a unitary acting on the ancilla sites). A more general state that's good to define is one where the quantum numbers of the ancilla sites are different from the ones on the physical site in such a way that the overall entangled state has well-defined total quantum numbers. And you can just make all the coefficients the same sign. For the spin case above, the resulting state would be the Sz=0 triplet state rather than the singlet state.
you'll need to add two additional loops inside the code, where instead of writing each wf.set line by hand, you write something like:
for(int n : range1(d))
wf /= norm(wf); //normalize to 1.0
where above the d+1-n argument to s2 is because of the idea I mentioned in step #1 above about pairing physical-ancilla states together to make the total quantum number always the same (assuming here that your sites have the property that state n and state d+1-n together have a fixed number of particles).
Also note above d is the site dimension.
Even if d is very large, as you said is the case for you, the above code should run rather quickly. If it doesn't then you probably have a bug or we'll need to have a closer look.