There is a much easier way to find the Gaussian curvature of ellipsoid without going through the tedious computation. I found this from Do Carmo's classic book of differential geometry (Ex.21, p.172-173).
First of all, we prove the following:
Let $M$ be a regular orientable surface with orientation $\mathbf{N}$ (i.e. a unit normal vector field). Let $f$ be a nowhere-zero smooth function on $M$. Let $p\in M$ and $\mathbf{v_1}$,$\mathbf{v_2}\in T_p(M)$ such that $\mathbf{v_1}$ and $\mathbf{v_2}$ are orthonormal, and that $\mathbf{v_1} \times\mathbf{v_2}=\mathbf{N}$. Then the Gaussian curvature of $M$ at $p$ is
\begin{align}
K=\frac{<d(f\mathbf{N})(\mathbf{v_1})\times d(f\mathbf{N})(\mathbf{v_2}),f\mathbf{N}>}{f^3}
\end{align}
where $d(f\mathbf{N})$ is the differential of $f\mathbf{N}$ at $p$.
To prove the above proposition, one has to start with the right-hand side of the equality and resort to the definition of differential, followed by the use of product rule of usual differentiation. After that, plug it into the expression and perform algebraic manipulation to finally cancel the function $f$. At last, consider $K$ as the determinant of the shape operator, and express every remaining vector in terms of the orthonormal basis $\mathbf{v_1}$,$\mathbf{v_2}$ in order to finish the proof.
Now, consider $M$ as the ellipsoid $\displaystyle g(x,y,z):=\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1$. Notice that the gradient of the function $h$ is normal to the tangent plane of $M$ at each point. So we have
\begin{align}
&\nabla{h}=2(\frac{x}{a^2},\frac{y}{b^2},\frac{z}{c^2}) &|\nabla{h}|=2\sqrt{\frac{x^2}{a^4}+\frac{y^2}{b^4}+\frac{z^2}{c^4}}=2f
\end{align}
where we let $f$ to be the restrction of the function
\begin{align}
\sqrt{\frac{x^2}{a^4}+\frac{y^2}{b^4}+\frac{z^2}{c^4}}
\end{align}
to $M$. Then we obtain
\begin{align}
N=\frac{\nabla{h}}{|\nabla{h}|}=\frac{1}{f} (\frac{x}{a^2},\frac{y}{b^2},\frac{z}{c^2})
\end{align}
It is easy to check that
\begin{align}
d(f\mathbf{N})=
\begin{pmatrix}
\frac{1}{a^2} & 0 & 0 \\
0 & \frac{1}{b^2} & 0 \\
0 & 0 & \frac{1}{c^2} \\
\end{pmatrix}
\end{align}
Let $\mathbf{u}$ and $\mathbf{v}$ be an orthonormal basis for $T_pM$ such that $\mathbf{u}\times\mathbf{v}=\mathbf{N}$. Then we obtain
\begin{align}
&d(f\mathbf{N})(\mathbf{u})=(\frac{u_1}{a^2},\frac{u_2}{b^2},\frac{u_3}{c^2}) &d(f\mathbf{N})(\mathbf{v})=(\frac{v_1}{a^2},\frac{v_2}{b^2},\frac{v_3}{c^2})
\end{align}
Using the formula,
\begin{align}
K&=\frac{1}{f^3a^2b^2c^2}\begin{vmatrix}
u_1 & u_2 & u_3 \\
v_1 & v_2 & v_3 \\
x & y & z \\
\end{vmatrix}\\
&=\frac{1}{f^3a^2b^2c^2}<\mathbf{N},\mathbf{p}> \\
&=\frac{1}{f^4a^2b^2c^2}(\frac{x^2}{a^2},\frac{y^2}{b^2},\frac{z^2}{c^2}) \\
&=\frac{1}{f^4a^2b^2c^2}
\end{align}
at point $\mathbf{p}=(x,y,z)\in{M}$
http://www.ann.jussieu.fr/~frey/papers/meshing/Goldman%20R.,%20Curvature%20formulas%20for%20implicit%20curves%20and%20surfaces.pdf
– sid Oct 26 '13 at 21:56