fixed_discrete_covariance: get diagonal indices faster

Profiling repeated calls to model(), way too much time was being spent on `diagonal = diag(K)`. This way doesn't copy the diagonal.