I see this pattern a lot, but it implicitly relies on `T`

being closed under `f`

. This is a reasonable assumption, but can quickly break down when relying eg on multiple AD packages at the same time.

The fundamental problem is that in general, it is very hard to preallocate for an operation

```
d = f(a, b, c)
```

as

```
f!(d, a, b, c)
```

when `a`

, `b`

and `c`

can be sufficiently generic, because predicting the type of `d`

can be tricky.

My current approach to this problem is to

- initially avoid it, not think about preallocation,
- the first line of optimizations is to use immutables, eg
`StaticArrays`

, and hope for the best,
- if I decide I really need it, use some kind of a heuristic to determine the element type.

For the last one, something like

```
T = typeof(g(one(eltype(a)), one(eltybe(b)), one(eltype(c)))
d = Matrix{T}(undef, I_know, the_dimensions)
```

where `g`

somehow reflects the elementary operations for what is going on in `f`

. But this can be brittle as the `<: Real`

interface is not formally specified and packages can deviate from expectations.