ENH: for one variable data, np.cov should return either a scalar or a (1, 1) array, not a ndim 0 array.

As it stands,

>>> np.cov([[1, 2, 3]])
array(1.)

I feel that this is a bit unnatural, considering that any input array of first dimension > 1 will return an ndim 2 array.

For example,

>>> np.cov([[1, 2, 3], [4, 5, 6]])
array([[1., 1.],
       [1., 1.]])

I suggest that either:

1. A scalar is returned. This would be useful for those people that prefer to use a single function for everything, and thus use cov for computing variance as well.
2. A (1, 1) array is returned. This would provide greater consistency with what is returned. Of these two, I prefer this one. I mean, np.dot doesn't squeeze, so why should cov?

While ndim 0 arrays can be used as scalars, it feels weird.

For the lazy, https://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html

Reproduce the two np.cov examples from the issue and read the linked np.cov documentation to confirm the current shape behavior. The issue leaves scalar versus (1, 1) unresolved; first get a decision on the API contract, then add coverage for the chosen one-variable result while preserving existing multi-variable behavior.