Hi all,

I’m trying to compute every product of two functions of a set of different functions `fv`

(e.g., a set of wavefunctions) at different `x`

to obtain a matrix `y`

for each `x`

of shape `NxN`

where `N`

is the number of functions. Here, `fv`

is just a vector of arbitrary functions.

The best I could come up with is the following:

```
using BenchmarkTools
import Base.*
*(f::Function, g::Function) = @. (x...) -> f(x...) * g(x...)
function get_function_matrix(fv::Vector{Function})
N = length(fv)
Mf = Array{Function, 2}(undef, N, N)
for j in 1:N, i in 1:N
Mf[i, j] = fv[i] * fv[j]
end
return Mf
end
function evaluate_function_matrix(Mf::Array{Function, 2}, x)
N = size(Mf, 1)
y = zeros(length(x), N, N)
for j in 1:N, i in 1:N
y_ = @view y[:, i, j]
map!(Mf[i, j], y_, x) # <- Why does cause an allocation?
end
return y
end
x = rand(100)
fv = [sin, cos, sinh, cosh, tan, tanh]
Mf = get_function_matrix(fv)
@btime y = evaluate_function_matrix(Mf, x)
```

My questions are:

1.) Why does the `map!`

function cause an allocation and is there a way to optimize this? Am I wrong to expect that this step should not cause any allocations at all?

2.) If one replaces the `map!`

function with `y[i, j] = M[i, j](x), where `

x` is just a float, the allocations increase even more which I don’t understand.

3.) Is there a more Julian way to solve this problem?

Thank you very much in advance for your help!