all representations by ternary quadratic form with nonnegative coefficients

Hello,

Is there an efficient way to find all solutions in nonnegative integers
x,y,z to
xy + xz + yz = n
for a (large) given integer n?
Under efficiency I understand anything that is noticeably better than just
fixing the value of one (smallest) of the variables and solving the
resulting bivariate quadratic equation by factorization. This approach
requires O(sqrt(n)) iterations.

Would it be a binary form, there is a handy qfbsolve() whose flags enable
finding all solutions, but for a more generic qfsolve() supports finding
just a single solution. For the given example, finding one solution is
trivial - it's (0, 1, n).
There is a way to parametrize all the solutions via qfparam(), but I don't
see how to get just nonnegative integer ones out of it.
Any hints would be appreciated.

Regards,
Max