This is how standard floating point computations work.

In Mathematica floating point calculations do not work like they do in C or with IEEE floating point values. This is due to Mathematica's roots as a computer algebra system. There are many differences, including the lack of separation between +0 and -0, the fact that there isn't only NaN, +inf, -inf, but a more general DirectedInfinity (try 1/0), etc. Overflow or underflow doesn't result in infinity or zero. Instead "machine numbers" get converted to "arbitrary precision numbers" which can represent much larger or smaller values with any number of digits. The calculation can continue, however, it will be slower.

I suppose Indeterminate is represented as some kind of NaN when it is in a PackedArray although the same is not done for Matrix.

There are multiple misunderstandings here. There is no such thing as Matrix in Mathematica, only nested Lists. If the nested lists have the same structure as an n-dimensional array, then this data structure may become a "packed array". This simply means that it uses the most efficient possible internal presentation.

Packed arrays cannot contains Indeterminate. They can only contain machine numbers. Indeterminate is a symbol. Inserting a symbol unpacks the array.

They also cannot contain arbitrary precision numbers. If such a result appears, it unpacks the array.

Overall I agree with you that there are some cases when the ability to store either missing values, indeterminates, or infinities in packed arrays would be nice. R can handle missing values (NA) without trading off storage efficiency, and this is frequently used by statisticians. But in Mathematica currently it is not possible.

Will this ever be added in the future? Given what I know about how Mathematica works, I highly doubt it. It is my impression that adding this would be a huge and highly risky task, which would also generate a lot of extra work for the future (more special cases will need to be handled by each new function that is added to the language).

That is how it is used, but one should prevent it in my opinion. I never make use of NaNs (Indeterminate), and they are really not necessary. If you want to use packed array, you'd better prevent them before doing e.g. division as they can't be packed.

Developer`ToPackedArray will return the original input not a packed version if one tries to pack an indeterminate or so...

There is no NaN or equivalent that can be stored in packed arrays. Indeterminate, Infinity, DirectedInfinity, and Missing are your options.

I'm not sure why you would have these 'as a result of a computation'; with all the computations I do (which is quite a bit) I never 'needed' NaNs. I think it is just bad practice in general. divisions by zero (main cause) should be caught and be dealt with accordingly.

Edit: There is a 'NaN' but you need to load a (built-in) package:

Needs["ComputerArithmetic`"]
Developer`PackedArrayQ[Developer`ToPackedArray[{1, 2, 3, NaN}]]

But, as you can see, it can't be used inside a PackedArray. So then one can perhaps better use Indeterminate, Infinity and so on...

Also note that:

var = Array[Function[x, 3.1], 1000000];
Developer`PackedArrayQ[var]

var = ConstantArray[3.1, 1000000];
Developer`PackedArrayQ[var]

Array does not (in general) return a packed array. You should use ConstantArray in such a case, which does return a packed array.