Allow matrix structure in covariance matrices to be exploited

This PR provides adds a new method `covariance_matrix_solve` to the CovStruct class that solves a system of equations whose coefficient matrix is the covariance matrix represented by the CovStruct instance.

Currently, we construct this covariance matrix explicitly, and solve these systems with a general-purpose solver.  However, for many of the dependence structures, the linear algebra can be optimized to exploit the structure of the matrix.

We provide a default implementation of `covariance_matrix_solve` that uses a general-purpose solver.  

We override this with optimized methods for the independence, autoregressive, and exchangeable cases.

The speed improvement would be most noticeable when the cluster sizes are large.

Tests for the independence and exchangeable cases match results obtained from R (as did the earlier code).

We don't have a comparator for the autoregressive case, so I added a regression test.  It agrees with the results from the previous (non-optimized) implementation.