When should matrices have units of measurement? As a mathematician I think of matrices as $\mathbb{F}^{m\times n}$, where $\mathbb{F}$ is a field and usually $\mathbb{F} = \mathbb{R}$ or $\mathbb{F} = \mathbb{C}$. Units are not necessary.
However, some engineers have told me that they prefer matrices to have units such as metres, kg, pounds, dollars, etc. Assigning a unit of measurement to each entry to me seems restrictive (for instance if working with dollars then is $A^2$ allowed?).
Here are a few things that I would like to understand more deeply:


*

*Are there examples where it is more appropriate to work with matrices that have units?

*If units can only restrict the algebra, why should one assign units at all?

*Is there anything exciting here, or is it just engineers trying to put physical interpretations on to matrix algebra?   
Also, see: 
https://stackoverflow.com/questions/11975658/how-do-units-flow-through-matrix-operations
 A: The point of having units in computations is to restrict what one can do, such that does doesn't end up inadvertently writing down a computation whose result depends on one's choice of units in an unexpected way.
In mathematics this is not done a lot, but mathematics can certainly do it if we want to. What you do is to work with matrices (etc) over the field of rational functions of $n$ variables, and assign one of the variables to stand for each of your fundamental units. That is, we're working in the field of fractions of the polynomial ring $\mathbb R[\rm m, s, kg, A, \ldots]$.
The fraction field contains all of the (unitless) real numbers, as well as all of the unitful measuremens such as $2\,{\rm s}$ or $4000\,{\rm m}$ or $35\,\rm{\frac m{s^2}}$. It also contains a lot of nonsensical elements such as $\frac{42\,\rm{kg}}{3\,{\rm m}+8\,{\rm s}}$, but that doesn't really bother us, because we know that the field-of-fractions construction still keeps the entire system consistent.
We can then define that the computation we're speaking of is well-formed if we can prove that the output always can be written as a real vector (or matrix or whatever) times a scalar such as $1\,{\rm \frac ms}$, which defines the dimension we expect (on physical grounds) that the answer must have.
Once this proof (known to physicists and engineers as dimensional analysis) has been carried out, we can then use an evaluation homomorphism to map all of the unit-variables to $1$, so they disappear from the formulas and the actual calculation can be done on ordinary real numbers.
Whether this is exciting or not is subjective.
A: Changing the units for a vector of measurements, e.g.
$$
x=\pmatrix{1\,{\rm ft}\\1\,{\rm lb}\\1\,{\rm hr}}
\implies x'=\pmatrix{0.3048\,{\rm m}\\453.59\,{\rm g}\\3600.0\,{\rm s}}
$$
is accomplished by multiplying the vector
by a diagonal matrix of scaling factors, i.e.
$${F={\rm Diag}\big(0.3048,\,453.59,\,3600\big),\qquad x' = Fx}$$
Suppose that you have a linear model $\,y=Ax\,$ and wish to convert the units of both the input and output vector
$$\eqalign{
y' &= Ey,\qquad x' = Fx,\qquad A' = EAF^{-1} \\
}$$
So that the model equation in the new system of units becomes
$$\eqalign{
y' &= A'x' \\
}$$
If $A^{-1}$ exists, then the answer will be the same no matter which system is solved.
It is however, extremely common when fitting such a model for the system to be over-determined (rectangular) and consequently, it is solved in a least-squares sense. In such situations, you will obtain very different answers depending on which units are being used.
That's because numerically $(3600\,{\rm sec})^2$ generates a much larger residual than $(1\,{\rm hr})^2$ despite the fact that they are identical when units are taken into account.
To address such situations, Jeffrey Uhlmann has proposed a unit-consistent inverse. This is a generalized inverse which satisfies the first two Penrose conditions
$$\eqalign{
AGA = A,\qquad GAG=G
}$$
as well as ensuring that $G$ is invariant with respect to
arbitrary scalings of the rows/columns of $A$.
A: The book by George W. Hart mentioned in the comments is actually a quite good reference on the topic.
Furthermore, I also gave a talk on how a C++ library can be designed that is able to solve the topic for the general case of non-uniform physical units in vectors and matrices:
https://www.youtube.com/watch?v=J6H9CwzynoQ
A: Think about $Ax=b$ (not solving it, just the equation itself). If $x_j$ and $b_i$ have units, then $a_{ij} x_j$ has the units of $b_i$. Hence $a_{ij}$ has the units of $b_i/x_j$. So if $x$ is a vector of times and $b$ is a vector of positions then the entries of $A$ are velocities (or at least have those units).
I do not really understand how to interpret the idea that the matrix itself has units, however.
A: The questions of the OP seem to me based on a misconception: that physicists, chemists and engineers --for some mysterious reason-- prefer to work with matrices where the entries are not numbers (as is mathematics), but have units as well ! But that is not the case at all. 
Physicists, chemists and engineers study processes in the real world. All  quantities and observables are measured and/or defined in terms of fundamental units. Now sometimes it can be useful to construct a mathematical model to simulate these processes, and to predict the outcome of future experiments. The variables and parameters in the mathematical model are not numbers but maintain their unit (in other words: dimension). 
If a model is complicated, containing many parameters and variables, it can be useful to introduce matrices and vectors. This is done as a book-keeping device, and to keep the notation concise. The fact that the entries of the matrix have a unit is just a side-effect.
Occasionally physicists use matrices to group certain variables together, just to make the point that these variables are related to each other in an unexpected yet profound way. For example in special relativity the $4$-vector $(x, y , z, t)$ is used. It symbolizes that space and time are intimately related. The fact that time has a different dimension than a spatial coordinate is sidestepped by introducing $ct$ instead of $t$. This is a natural choice, since $c$ --the speeds of light-- takes a very special place in relativity.      
