What is the usefulness of the Wigner-Eckart theorem?
I will not get into theoretical details -- Luboš ad Marek did that better than I'm able to.
Let me give an example instead: suppose that we need to calculate this integral:
$\int d\Omega (Y_{3m_1})^*Y_{2m_2}Y_{1m_3}$
Here $Y_{lm}$ -- are spherical harmonics and we integrate over the sphere $d\Omega=\sin\theta d\theta d\phi$.
This kind of integrals appear over and over in, say, spectroscopy problems. Let us calculate it for $m_1=m_2=m_3=0$:
$\int d\Omega (Y_{30})^*Y_{20}Y_{10} = \frac{\sqrt{105}}{32\sqrt{\pi^3}}\int d\Omega \cos\theta\,(1-3\cos^2\theta)(3\cos\theta-5\cos^3\theta)=$
$ = \frac{\sqrt{105}}{32\sqrt{\pi^3}}\cdot 2\pi \int d\theta\,\left(3\cos^2\theta\sin\theta-14\cos^4\theta\sin\theta+15\cos^6\theta\sin\theta\right)=\frac{3}{2}\sqrt{\frac{3}{35\pi}}$
Hard work, huh? The problem is that we usually need to evaluate this for all values of $m_i$. That is 7*5*3 = 105 integrals. So instead of doing all of them we got to exploit their symmetry. And that's exactly where the Wigner-Eckart theorem is useful:
$\int d\Omega (Y_{3m_1})^*Y_{2m_2}Y_{1m_3} = \langle l=3,m_1| Y_{2m_2} | l=1,m_3\rangle = C_{m_1m_2m_3}^{3\,2\,1}(3||Y_2||1)$
$C_{m_1m_2m_3}^{j_1j_2j_3}$ -- are the Clebsch-Gordan coefficients
$(3||Y_2||1)$ -- is the reduced matrix element which we can derive from our expression for $m_1=m_2=m_3=0$:
$\frac{3}{2}\sqrt{\frac{3}{35\pi}} = C_{0\,0\,0}^{3\,2\,1}(3||Y_2||1)\quad \Rightarrow \quad (3||Y_2||1)=\frac{1}{2}\sqrt{\frac{3}{\pi}}$
So the final answer for our integral is:
$\int d\Omega(Y_{3m_1})^*Y_{2m_2}Y_{1m_3}=\sqrt{\frac{3}{4\pi}}C_{m_1m_2m_3}^{3\,2\,1}$
It is reduced to calculation of the Clebsch-Gordan coefficient and there are a lot of, tables, programs, reduction and summation formulae to work with them.
First, W-E theorem is just a simple (just bear with me, I know that the theorem can appear formidable if not explained properly) statement about the decomposition of the tensor product of representations into its irreducible components.
Suppose we have a group $G$ and a tensor operator $A_r$ transforming under a(n irreducible) representation $\Gamma^1$ of it ($r$ counting the number of components of the operator, e.g. 3 for the usual angular momentum in 3 dimensions). Also suppose that you have additional (irreducible) representation $\Gamma^2$ with basis $\left\{ \left| \psi_n \right>\right\}$. It is easy to check that then vectors $\Psi_{rn} \equiv A_r \left|\psi_n\right>$ transform as a tensor product $\Gamma^1 \otimes \Gamma^2$.
In the following we use that the group $G$ is compact in order to be able to decompose its representations into irreducible components. This can be weakened but in general representations of non-compact groups don't have to be reducible to their irreducible components; there's pretty subtle mathematics behind representation theory of non-compact groups; even seemingly simple ones as $SL(2, \mathbb{R})$.
So, we'll decompose the representation as $$\Gamma^1 \otimes \Gamma^2 = \bigoplus_{\alpha} \Gamma^{\alpha}$$ where $\alpha$ runs over irreducible representations of $G$ (possibly repeated). This amounts to finding a more suitable basis for vectors $\Psi_{rn}$. We will write $\Phi_{\alpha m}$ for that basis (with $m$ indexing the components of representation $\Gamma^{\alpha}$). Then we can write $$\Psi_{rn} = \sum_{\alpha, m} U^{\alpha m}_{rn} \Phi_{\alpha m}$$.
Now, back to the problem: we are interested in computing some element such as $\left<\omega_k \right| A_r \left | \psi_n\right>$ with $\omega_k$ transforming in some $\Gamma^3$ representation. Thanks to Schur orthogonality relations and the fact that we decomposed $A_r \left | \psi_n\right>$ into its irreducible components we can see that for this element to be non-zero, there has to $\Gamma^3$ representation in the $\Gamma^{\alpha}$ decomposition. That already saves us a lot of computation. If there is no $\Gamma^3$ present there then we're finished, all those elements will be zero. And if it is present, we only have to carry out calculations corresponding to just the $\Gamma^3$ part of the decomposition (which will usually be a small part of the total).
Okay, I think the above might have been a little confusing, so let's try an example. The canonical one is with $SO(3)$ of course. Suppose we want to compute something like $\left<j m\right| \mathbf X \left|j' m'\right>$ with $\mathbf X$ being the position operator (which transforms under $SO(3)$ vector irrep). So we are interesting of doing a tensor product ${\mathbf 3} \otimes (\mathbf{2j +1})$ which can be shown to be equal to $(\mathbf{2j-1}) \oplus (\mathbf {2j+1}) \oplus (\mathbf {2j+3})$ (supposing $j$ is high enough for simplicity). So $j'$ has to be equal to one $j-1, j, j+1$ and moreover $m = m_X + m'$ (if $\mathbf X$ is written in a diagonal basis of the vectorial representations one can decompose it into operators having eigenvalues of m_X = $-1, 0, 1$) will need to hold (this is again a consequence of orthogonality relations).
Anyway, the most important point (at least for me) is about the decomposition of the tensor product into irreducible components. Supposing you can carry out this decomposition (which can be pretty hard at times), the calculation of matrix elements will simplify greatly and you can immediately decide stuff such as whether some molecule with a given symmetry will radiate in IR (which amounts to computing matrix elements of a dipole operator between eigenstates of that molecule). I am sure you can imagine lots of other similar stuff you may compute yourself. Applications of this theorem are pretty much limitless.
The Wigner-Eckart Theorem
http://en.wikipedia.org/wiki/Wigner_Eckart_theorem
is a formula that tells us about all "simple constraints" that group theory - the mathematical incarnation of the wisdom about symmetries, especially in the $SO(3)\approx SU(2)$ case (rotations in a three-dimensional space) - implies about matrix elements of tensor operators - those that transform as some representation of the same symmetry.
The dependence on the indices labeling basis vectors of the representations - $m_1, m_2, m_3$ in the $SU(2)$ case - is totally determined. Instead of thinking that some matrix elements depend on many variables, physicists may realize that the symmetry guarantees that the matrix elements only depend on a few labels labeling the whole "multiplets" of the states and operators rather than on all the labels identifying the individual components. It's always critically important to know how much freedom or how much uncertainty there is about some observables - for example, experimenters don't want to repeat their experiment $(2J+1)^3$ times without a good reason - and we would get a totally wrong idea without this theorem.
To see that the theorem is used all the time, check e.g. these 5700 papers
http://scholar.google.com/scholar?q=wigner-eckart+theorem&hl=en&lr=&btnG=Search
many of which are highly cited ones. The topics of the papers include optics, nanotubes, X-rays, spectroscopy, condensed matter physics, mathematical physics involving integrable systems, quantum chemistry, nuclear physics, and virtually all other branches of physics that depend on quantum mechanics. In many other cases, the theorem is being used without mentioning its name, and its generalizations are being used all the time in advanced theoretical physics, e.g. in the contexts with groups that are much more complicated than $SO(3)\approx SU(2)$.
It's interesting to note that in some sense, quantum mechanics allows the symmetry to impose as many constraints as classical physics. In classical physics, we could start with an initial state labeled by numbers $I_i$, apply some operations depending on parameters $O_j$, and we would obtain a final state described by parameters $F_k$. Classical physics would tell us Yes/No - whether we can get from $I_i$ via $O_j$ to $F_k$: that's the counterpart of the quantum probability amplitude.
The rotational $SO(3)$ symmetry - which has 3 parameters (3 independent rotations - or the latitude; longitude of the axis, and the angle) would only tell us that if we rotate all objects $I_i,O_j,F_k$ by the same rotation, we obtain a valid proposition again (Yes goes to Yes, No goes to No). So the dependence on 3 parameters - corresponding to the 3 rotations - is eliminated. In quantum mechanics, we also eliminate the dependence on 3 parameters - in this case $m_1,m_2,m_3$, the projections $j_z$ for the two state vectors and for the operator sandwiched in between them. In some proper counting, this is true for any $d$-dimensional group of symmetries, I think.