Cannot work with matrix data

[thudjx]

Sorry to ask question in the issue page here, but I really have struggled with the problem following for a time and cannot figure it out.

I am trying to use NLsolve to perform least-square-fitting. The kpts acts as xs in normal sense, which is of shape (n, 2) with each row be a point in 2d. The data to be fitted are the eigenvals at those points, which is obtained through dirac_2layer:

function dirac_1layer(kpts, p)
	γ0, V1 = p
	eigs = Array{Float64}(undef, 2, size(kpts, 1))
	for ik in axes(eigs, 2)
		f = f_gen(vec(kpts[ik, :]))
		fc = conj(f)
		h11 = [0 -γ0*f; -γ0*fc 0]
		hpot = diagm([V1, V1])
		eigs[:, ik] = eigvals(h11 .+ hpot)
	end
	return vec(eigs)
end

Random.seed!(1234)
# - initial values of params -
γ0 = 2.5
V10 = -0.097611
# - generate testing data -
kpts = rand(100, 2)
ys0 = dirac_1layer(kpts, [γ0, V10])
ys = ys0 .+ randn(size(ys0)) .* 0.1	
# - fit -
function residual!(dy, p)
	dy .= dirac_1layer(kpts, p) .- ys
end
nlsolve(residual!, [γ0, V10])

The error is like:

In fact, I also have tried the not in-place syntax:

...
# - fit -
function residual(p)
	return dirac_1layer(kpts, p) .- ys
end
	
nlsolve(residual, [γ0, V10])

And the error occur at the same line:

Is there something I am missing here?