AM-GM-HM Triplets
One way to look at this problem is to look for the extremes of the geometric mean given fixed arithmetic and harmonic means.
The following inequalities can be verified by taking $\log$s and using Jensen: $$ \frac{1}{\int\frac{1}{x}\;\mathrm{d}\mu}\le\exp\left(\int\log(x)\;\mathrm{d}\mu\right)\le\int x\;\mathrm{d}\mu\tag{1} $$ That is, $HM\le GM\le AM$.
Scaling the data scales all three means by the same factor, so let us assume that $$ \int x\;\mathrm{d}\mu=1\quad\text{(arithmetic mean }=1\text{)}\tag{2} $$ and $$ \int\frac{\mathrm{d}\mu}{x}=\frac{1}{h}\quad\text{(harmonic mean }=h\text{)}\tag{3} $$ and attempt to find the critical values for the $\log$ of the geometric mean: $$ \int\log(x)\;\mathrm{d}\mu\tag{4} $$ Taking variances of $(2)$, $(3)$, and $(4)$, we get that for all variations so that the $AM$ is stationary $$ \int\delta x\;\mathrm{d}\mu=0\tag{5} $$ and the reciprocal of the $HM$ is stationary, $$ \int\frac{\delta x}{x^2}\mathrm{d}\mu=0\tag{6} $$ we should also have that the $\log$ of the geometric mean is stationary: $$ \int\frac{\delta x}{x}\mathrm{d}\mu=0\tag{7} $$ Linearity considerations say that $(5)$, $(6)$, and $(7)$ imply that there are constants $a$ and $b$ so that $$ \frac{1}{x}=a+\frac{b}{x^2}\tag{8} $$ Multiplying $(8)$ by x and combining with $(2)$ and $(3)$ yields $$ 1=a+\frac{b}{h}\tag{9} $$ Solving $(8)$ for $x$ and $(9)$ for $b$, we get $$ x=\frac{1\pm\sqrt{1-4ha(1-a)}}{2a}\tag{10} $$ To maintain $(2)$, we need a combination of the two values of $x$: $$ \lambda\frac{1+\sqrt{1-4ha(1-a)}}{2a}+(1-\lambda)\frac{1-\sqrt{1-4ha(1-a)}}{2a}=1\tag{11} $$ Solving $(11)$ for $\lambda$, we get $$ \lambda=\frac{1}{2}\left(1+\frac{2a-1}{\sqrt{1-4ha(1-a)}}\right)\quad\text{and}\quad a=\frac{1}{2}\left(1+\frac{2\lambda-1}{\sqrt{1+4\frac{h}{1-h}\lambda(1-\lambda)}}\right)\tag{12} $$ Note that $(12)$ shows that $\lim\limits_{a\to0}\frac{\lambda}{a}=\lim\limits_{a\to1}\frac{1-\lambda}{1-a}=1$.
For a given harmonic mean, $h$, and split between values of $x$, $\lambda$, $(10)$ and $(12)$ yield that the stationary geometric mean is $$ \begin{align} g(h,\lambda) &=\left(\frac{1+\sqrt{1-4ha(1-a)}}{2a}\right)^\lambda\left(\frac{1-\sqrt{1-4ha(1-a)}}{2a}\right)^{1-\lambda}\\ &=\frac{\left(\sqrt{1+4\frac{h}{1-h}\lambda(1-\lambda)}+1\right)^\lambda\left(\sqrt{1+4\frac{h}{1-h}\lambda(1-\lambda)}-1\right)^{1-\lambda}}{\sqrt{1+4\frac{h}{1-h}\lambda(1-\lambda)}+2\lambda-1}\tag{13} \end{align} $$ For each $h\in[0,1]$, $g(h,\lambda)$ is monotonic in $\lambda$ (verified below) and $$ \lim_{\lambda\to0}g(h,\lambda)=h\qquad\text{and}\qquad\lim_{\lambda\to1}g(h,\lambda)=1\tag{14} $$ With a countinuous selection of $\lambda$, we can attain any value of $g(h,\lambda)$ between $h$ and $1$. However, When dealing with $n$ numbers, $\lambda$ can only take the values $\frac{1}{n}\dots \frac{n-1}{n}$. Thus, with a harmonic mean of $h$, the geometric mean can only vary between $g(h,\frac{1}{n})$ and $g(h,\frac{n-1}{n})$.
$\mathbf{n}$ objects with arithmetic mean $\mathbf{=1}$:
Proof that for any $\mathbf{h}$, $\mathbf{g(h,\lambda)}$ is monotonic in $\mathbf{\lambda}$:
Fix $h$. Define $$ \Delta=\sqrt{1+4\frac{h}{1-h}\lambda(1-\lambda)}\tag{15} $$ Using $(13)$ and $(15)$, we get that $$ g(h,\lambda)=\frac{(\Delta+1)^\lambda(\Delta-1)^{1-\lambda}}{\Delta+2\lambda-1}\tag{16} $$ Equation $(15)$ also implies that $$ (1-h)\Delta^2+h(2\lambda-1)^2=1\tag{17} $$ In consideration of $(17)$, define $\eta$ and $\theta$ by $$ \begin{align} \sin(\eta)&=\sqrt{h}\\ \cos(\eta)&=\sqrt{1-h} \end{align}\tag{18} $$ and $$ \begin{align} \sin(\theta)&=\sin(\eta)\;(2\lambda-1)\\ \cos(\theta)&=\cos(\eta)\;\Delta \end{align}\tag{19} $$ where $0\le\lambda\le1$ and thus $-\eta\le\theta\le\eta$. Then using the identities $$ \begin{array}{ccccc} \cos(\eta)(\Delta+1)&=&\cos(\theta)+\cos(\eta)&=&2\cos\left(\frac{\eta+\theta}{2}\right)\cos\left(\frac{\eta-\theta}{2}\right)\\ \cos(\eta)(\Delta-1)&=&\cos(\theta)-\cos(\eta)&=&2\sin\left(\frac{\eta+\theta}{2}\right)\sin\left(\frac{\eta-\theta}{2}\right)\\ \sin(\eta)2\lambda&=&\sin(\eta)+\sin(\theta)&=&2\sin\left(\frac{\eta+\theta}{2}\right)\cos\left(\frac{\eta-\theta}{2}\right)\\ \sin(\eta)2(1-\lambda)&=&\sin(\eta)-\sin(\theta)&=&2\cos\left(\frac{\eta+\theta}{2}\right)\sin\left(\frac{\eta-\theta}{2}\right) \end{array}\tag{20} $$ we get that $$ \begin{align} &\left(\frac{1}{\sin(\eta)}\frac{(\Delta+1)^\lambda(\Delta-1)^{1-\lambda}}{\Delta+2\lambda-1}\right)^{2\sin(\eta)}\\ &\vphantom{\Huge{\dfrac{A}{A}}}=\frac{(\cos(\theta)+\cos(\eta))^{\sin(\eta)+\sin(\theta)}(\cos(\theta)-\cos(\eta))^{\sin(\eta)-\sin(\theta)}}{\left(\sin(\eta)\cos(\theta)+\cos(\eta)\sin(\theta)\right)^{2\sin(\eta)}}\\ &\vphantom{\Huge{\dfrac{A}{A}}}=\frac{(2\cos(\frac{\eta+\theta}{2})\cos(\frac{\eta-\theta}{2}))^{\sin(\eta)+\sin(\theta)}(2\sin(\frac{\eta+\theta}{2})\sin(\frac{\eta-\theta}{2}))^{\sin(\eta)-\sin(\theta)}}{\sin(\eta+\theta)^{2\sin(\eta)}}\\ &\vphantom{\Huge{\dfrac{A}{A}}}=\frac{(\sin(\eta+\theta)\sin(\eta-\theta))^{\sin(\eta)}(\cot(\frac{\eta+\theta}{2})\cot(\frac{\eta-\theta}{2}))^{\sin(\theta)}}{\sin(\eta+\theta)^{2\sin(\eta)}}\\ &\vphantom{\Huge{\dfrac{A}{A}}}=\left(\frac{\sin(\eta-\theta)}{\sin(\eta+\theta)}\right)^{\sin(\eta)}\left(\frac{\cos(\theta)+\cos(\eta)}{\cos(\theta)-\cos(\eta)}\right)^{\sin(\theta)}\tag{21} \end{align} $$ The logarithmic derivative of $(21)$ is $$ \begin{align} &\frac{\mathrm{d}}{\mathrm{d}\theta}\left(\sin(\eta)\log\left(\frac{\sin(\eta-\theta)}{\sin(\eta+\theta)}\right)+\sin(\theta)\log\left(\frac{\cos(\theta)+\cos(\eta)}{\cos(\theta)-\cos(\eta)}\right)\right)\\ &\vphantom{\Huge{A}}=\sin(\eta)(-\cot(\eta-\theta)-\cot(\eta+\theta))\\ &+\cos(\theta)\log\left(\frac{\cos(\theta)+\cos(\eta)}{\cos(\theta)-\cos(\eta)}\right)\\ &+\sin(\theta)\left(\frac{-\sin(\theta)}{\cos(\theta)+\cos(\eta)}+\frac{\sin(\theta)}{\cos(\theta)-\cos(\eta)}\right)\\ &\vphantom{\Huge{\dfrac{A}{A}}}=-2\cos(\eta)+\cos(\theta)\log\left(\frac{\cos(\theta)+\cos(\eta)}{\cos(\theta)-\cos(\eta)}\right)\tag{22} \end{align} $$ Dividing $(22)$ by $\cos(\theta)$ and setting $t=\dfrac{\cos(\eta)}{\cos(\theta)}$, yields that $(22)$ vanishes precisely when $$ -2t+\log\left(\frac{1+t}{1-t}\right)=0\tag{23} $$ However, $(23)$ vanishes only at $t=0$, but because $-\eta\le\theta\le\eta$, we have $t\ge\cos(\eta)=\sqrt{1-h}$. Therefore, $(22)$ doesn't vanish, and thus, $(21)$ is monotonic.
I feel like I'm bringing rather heavy machinery in this problem, but here goes.
Rahul Narain's guess is correct and follows from the Equal Variable Theorem as proved here
The complete theorem is much more general, but for your question Corollary 1.9, Case 1 suffices:
Let $a_1, a_2, \ldots, a_n$ $(n \geq 3)$ be fixed non-negative numbers, let $0 \leq x_1 \leq x_2 \leq \cdots \leq x_n$ such that $$ x_1+x_2+\cdots+x_n = a_1+a_2+\cdots+a_n,$$ $$x_1^p + x_2^p +\cdots+x_n^p = a_1^p + a_2^p +\cdots+ a_n^p,$$ and let $E = x_1^q + x_2^q +\cdots+x_n^q$.
Case 1. $p \leq 0$ ($p=0$ yields $x_1x_2\cdots x_n = a_1a_2\cdots a_n > 0$).
(a) For $q \in (p,0) \cup (1, \infty)$, $E$ is maximal when $0 \lt x_1 = x_2 = \cdots = x_{n-1} \leq x_n$, and is minimal when $0 < x_1 \leq x_2 = x_3 = \cdots = x_n$.
(b) For $q \in (-\infty,p) \cup (0,1)$, $E$ is minimal when $0 \lt x_1 = x_2 = \cdots = x_{n-1} \leq x_n$, and is maximal when $0 < x_1 \leq x_2 = x_3 = \cdots = x_n$.
Now take $p=0$ for the geometric mean, $q = -1$ for the harmonic mean and use part b) from above to arrive at the claim.
OK, I have the right proof now. Once I finish this one, I'll delete the others. I am proving Rahul's intuition, that the extreme value of the product occurs when all the $x_i$ but one are equal. More specifically, I will prove that the product is maximized when the multiple value is larger than the single value.
Lemma: Let $a$, $b$ and $c$ be positive reals; define $n = a+b+c$; let $P$ and $Q$ be positive reals with $PQ > n^2$. Then, on the curve $$\begin{matrix} P &=& ax+by+cz \\ Q &=& ax^{-1} + b y^{-1} + c z^{-1} \\ x,y,z &>& 0 , \end{matrix}$$ the maximal value of $x^a y^b z^c$ occurs at a point where two of $(x,y,z)$ are equal and the third is smaller.
Remark: If $PQ < n^2$, then the above equations have no real solutions; if $PQ=n^2$ then the only real solution is $x=y=z=P/n$. So the lemma discusses the interesting case.
Proof: First, notice that the equations above have all of their solutions within the box $[a/Q, P/a] \times [b/Q, P/b] \times [c/Q, P/q]$, and form a closed subset of that box. So the solution space is compact (closed subset of a bounded set), and the maximum must exist.
Let $(x,y,z)$ be a point where the two equations are satisfied. If we had $x=y=z$ then $PQ=n^2$, so it is not possible that all of $x$, $y$ and $z$ are equal.
The gradient of the function $f:=ax+by+cz$ is $(a,b,c)$. The gradient of $g:=a x^{-1} + b y^{-1} + c z^{-1}$ is $(-a x^{-2}, -b y^{-2}, -c z^{-2})$. Since $x$, $y$ and $z$ are not all equal, these vectors are not parallel. Thus, the two surfaces $f=P$ and $g=Q$ meet transversely, in a smooth curve. A tangent vector to that curve is $$(a,b,c) \times (-a x^{-2}, -b y^{-2}, -c z^{-2}) = {\LARGE (} bc (y^{-2} - z^{-2}), ac (z^{-2} - x^{-2}), ab(x^{-2} - y^{-2}) {\LARGE )} \quad (\ast)$$
Rather than optimizing $x^a y^b z^c$, we will optimize $h:= \log(x^a y^b z^c) = a \log x + b \log y + c \log z$. The derivative of $h$ along the vector $(\ast)$ is $$abc (x^{-1} y^{-2} - x^{-1} z^{-2} + y^{-1} z^{-2} - y^{-1} x^{-2} + z^{-1} x^{-2} - z^{-1} y^{-2} ) = \frac{abc}{x^2 y^2 z^2} (x-y)(x-z)(y-z). \quad (\ast\ast).$$
If $x$, $y$ and $z$ are all distinct, then $(\ast \ast)$ is nonzero, so moving along the curve $f=P$, $g=Q$ in one direction or the other will cause $h$ to increase. A more detailed analysis (omitted) shows that we have a local maximum when the largest number of $(x,y,z)$ occurs twice. $\square$
We now prove the result. Let $(x_1, \ldots, x_n)$ be any $n$-tuple of positive reals. We will show that there are $y$ and $z$ with $y+(n-1) z = \sum x_i$, $y^{-1} + (n-1) z^{-1} = \sum x_i^{-1}$, $y<z$ and $y z^{n-1} \geq \prod x_i$. In other words, $(y,z,z,z,\ldots, z)$ maximizes $\prod x_i$ subject to fixed values of $\sum x_i$ and $\sum x_i^{-1}$.
The proof proceeds in $n-2$ steps. At each step, we increase the number of times the maximal value occurs among $(x_1, x_2, \ldots, x_n)$ and increase $\prod x_i$, all while holding $\sum x_i$ and $\sum x_i^{-1}$ constant. So, suppose that $x_1=x_2=\cdots = x_k > x_{k+1}, \ldots, x_n$ for some $k$. If $k=n-1$ then we are done. If not, hold $x_{k+3}$ through $x_n$ fixed, while replacing $(x_1, \ldots, x_{k+2})$ by the values which maximize $\prod_{i=1}^{k+2} x_i$ while holding $\sum_{i=1}^{k+2} x_i$ and $\sum_{i=1}^{k+2} x_i^{-1}$ fixed. By the Lemma, the new maximal value will now occur $k+1$ times. After finitely many steps of this sort, the maximum value occurs $n-1$ times and we have reached the maximum.
If you want an explicit formula for the maximal value of the geometric mean, then you need to actually solve the equations $$\begin{matrix} y+(n-1) z &=& P \\ y^{-1} + (n-1) ^{-1} &=& Q \end{matrix}$$ This isn't too bad: Clearing denominators in the second equation and substituting in from the first gives $z+(n-1) (P - (n-1) z) = Q z(P-(n-1) z)$. This is a quadratic in $z$, so it is explicitly solvable, but I must say I don't find the solution illuminating. Anyway, take that solution and compute $y^{1/n} z^{(n-1)/n}$, and you will have computed the optimal value of the geometric mean when the arithmetic mean is $P/n$ and the harmonic is $n/Q$.