People have mentioned so far how dimensional analysis is fundamental in many inequalities in analysis, especially estimates (and even formulas) which come from partial differential equations. I want to elaborate on this point. For example, in ${\mathbb R}^{n + 1}$ ($n$ spatial dimensions), look at an equation like the wave equation
$(-\partial_t^2 + \Delta) u = f $ --
for which we have the formal units $U/T^2 \sim U/X^2 \sim f$ where $T$ is the unit for the time variable and $X$ is the unit for the space variable. or the Schrodinger equation
$(i \partial_t + \Delta)u = f $
(for which $U/T \sim U / X^2 \sim f$). The rigorous meaning of these units for Schrodinger is that for $\lambda \neq 0$, we can rescale a given solution to construct another solution $u(t/\lambda^2, x/\lambda)$ whose forcing term is $\lambda^{-2} f(t/\lambda^2, x/\lambda)$ (here we are imagining $\lambda$ has the units of $X$). (Similarly, scaling in the $x$ variable alone, we can change the equation to remove or reinsert the physical constants which usually appear in these equations.) Consider a large forcing term, e.g. $f$ belongs some $L^p$ space or mixed space-time $L^p$ space; you use dimensional analysis to figure out to what space we could possibly guarantee the solution $u$ lives in -- the goal here is to make rigorous the idea that $u$ depends continuously on $f$, and thanks to some abstract facts of functional analysis, proving this continuous dependence for a linear equation can essentially be achieved by no other means than by proving a (Strichartz) estimate such as
$$ \left(\int |u|^r \, dt \, dx\right)^{1/r} \leq C \left(\int |f|^p \,dt \, dx\right)^{1/p}$$
for some constant $C$ independent of $f$. Together with linearity, this estimate when applied to the difference of nearby $f$'s, says, if true, that as $f$ varies in $L^p$, $u$ varies continuously in $L^r$.
The function space norms themselves ``have units'': e.g. $\|u\|_{L^2} = (\int |u|^2 \,dt \,dx)^{1/2}$ has units of $u$ times $T^{1/2}X^{n/2}$ where $T$ are the time units and $X$ are the spatial units ($T \sim X$ for the wave equation, but $T \sim X^2$ for Schrodinger). One can use these "units" to figure out what kind of estimates might possibly be possible. E.g. for Schrodinger,
$$ \left(\int |u|^r \,dt \,dx\right)^{1/r} \leq C \left(\int |f|^p \,dt \,dx\right)^{1/p}$$
would be impossible if the units were compatible, the left hand side is like
$$U (TX^n)^{1/r} \sim UX^{(n+2)/r}$$
whereas the term with $f \sim (U/X^2)$ has dimensions $(U/X^2)\cdot(TX^n)^{1/p} = U(X^{-2+(n+2)/p})$.
So dimensional analysis tells us that we have no chance of proving this kind of estimate unless $(n+2)/p - 2 = (n+2)/r$; otherwise one could rescale any solution to produce a counterexample. Similar considerations apply to the problem of determining in what sense $u$ can depend continuously on its Cauchy data.
The dimensional analysis is very far from a proof of the estimate, which is not always true and requires a real understanding of, say, the dispersion properties of the specific equation at hand (you could have been looking at a different equation with the same scaling properties but which is qualitatively very different). Denis Serre mentioned how these dispersion properties can be read off from the curvatures of the parabola (for Schrodinger) and cone (for wave) in frequency space, although the proofs also bring in complex interpolation theory to also take into account the energy, and to get the most sharp results you need an interpolation scheme which treats individual frequency scales differently. You also need to look at norms which are different in the space and time variables, because the equation also distinguishes these variables. (The space of functions whose "energy" or "mass" is bounded in time gives a particularly natural spacetime norm.) Suffice to say, dimensional analysis is far from enough to justify a bound.
For nonlinear problems with some kind of scaling symmetry, the units of the solution and the time or space dimensions can be tied together. E.g. consider the incompressible Euler equation without force for an unknown velocity field $\partial_t u + (u \cdot \nabla) u = - \nabla p$, $\nabla \cdot u = 0$. Then formally $U/T \sim U^2/X \sim P/X$ or $U \sim X/T$ (which is good for a ``velocity'' field), which means that given a solution $u(t,x)$ you can rescale to obtain a new solution $\lambda^{b-a} u(t/\lambda^a, x/\lambda^b)$ with pressure $\lambda^{2(b - a)} p(t/\lambda,x/\lambda)$. After all, if $\lambda^b \sim X$ and $\lambda^a \sim T$, then this rescaled $u$ looks like has units of $U$ and $P$ appears to have units of energy $U^2$. Any estimate for the Euler equations has to be consistent with this scaling, so it has consequences for studying the equation (of course, not so many a priori estimates exist, which is another problem...). On the other hand, like Deane Yang said, it's very good for dummy-checks.
I guess all I've said is that the equations of physics often have a scaling symmetry, and due to this symmetry dimensional analysis is of use for rigorous mathematical treatment of the equation. (But that's no surprise since everyone's talking about the same equation.) But in mathematics the idea that inequalities must be consistent with a scaling symmetry goes a long way. For example, consider the Hölder inequality
$$\int f(x) g(x) \, d\mu(x) \leq \left( \int |f(x)|^p \, d\mu(x)\right)^{1/p} \left(\int |g(x)|^q \, d\mu(x)\right)^{1/q}$$
This inequality is valid for any measure $\mu$; so you could be integrating over a surface or just taking a discrete weighted sum or something and it's still true. In particular, it is true for $\mu$ and true for the renormalized measure $\mu/\lambda$. If you compare "units", the left hand side is like $f g \mu$ and the right hand side is like $f g \mu^{1/p + 1/q}$. Thus, for Hölder's inequality to be true for arbitrary measure spaces, we see we need $1/p + 1/q = 1$. Here I am talking about abstract measure spaces -- something which has nothing to do with a physical problem or Euclidean spaces. Terry Tao has harmonic analysis notes where he uses this symmetry to renormalize $\mu$ and reduce the proof of Hölder's inequality to the case where $g$ is not even present, which allows one to view Hölder's inequality as an interpolation statement. The point of view I've expressed in great part derives from and is elaborated in those notes, which can be found here.
But from a broader mathematical viewpoint, dimensional analysis is probably only one example of paying attention to a group of symmetries (not just scaling symmetries). But at the moment, I cannot think of a particularly good example to illustrate this point. Probably one can also be found in the linked notes.
