Multiplying float16 matrices gives a float32 matrix

[ViralBShah]

This should give a Float16 result.

julia> ones(Float16,5,5)*ones(Float16,5,5,)
5x5 Float32 Array:
 5.0f0  5.0f0  5.0f0  5.0f0  5.0f0
 5.0f0  5.0f0  5.0f0  5.0f0  5.0f0
 5.0f0  5.0f0  5.0f0  5.0f0  5.0f0
 5.0f0  5.0f0  5.0f0  5.0f0  5.0f0
 5.0f0  5.0f0  5.0f0  5.0f0  5.0f0

[StefanKarpinski]

I somewhat disagree. Float16 is only for storage, not for computation. In this case it's debatable, however.

[stevengj]

Julia seems pretty consistent in that all arithmetic operations on Float16 values produce Float32 results; why should matrix multiplication be different?

[JeffBezanson]

Another example:

julia> ones(Int8,5,5)*ones(Int8,5,5,)
5x5 Int64 Array:
 5  5  5  5  5
 5  5  5  5  5
 5  5  5  5  5
 5  5  5  5  5
 5  5  5  5  5

If we want to change Float16 ops to something like

+(x::Float16, y::Float16) = float16(float32(x)+float32(y))

then matrix multiply would also change accordingly.

[ViralBShah]

The documentation also does clearly say:

Half-precision floating-point numbers are also supported (``Float16``), but
only as a storage format. In calculations they'll be converted to ``Float32``::

    julia> sizeof(float16(4.))
    2

    julia> 2*float16(4.)
    8.0f0