Why not have a cov_exp_quad() with vector-valued rho?

I noticed that there is no cov_exp_quad() for a covariance of this form (using TeX notation):

(\alpha^2)*\exp[-0.5*\sum_{d=1}^D ((x_{i,d} - x_{j,d})/rho_d)^2]

only for the special case where rho_1 = rho_2 = … = rho. Was there a reason for handling only this particular case?

Yes. Difficulty in implementation at the C++ level. It was easy to start with a scalar rho. Now that we have scalar rho, someone needs to implement the container version. If you know some C++, would love the help.