I'm looking for sparse SPD matrices with right hand side? There is UF collection of sparse matrices, however, I'm not sure how do I search of the matrices of these kind efficiently (I'm doing a naive search which hasn't given me any results so far and it takes arbitrary long for some of the matrices).
[The comment thread is long enough to convert to an answer.]
We acknowledge that most matrix collections don't include physically reasonable right hand sides. If we want to test solver performance, we have to generate problems to solve --- either the right-hand sides $b$ or solutions $x$ from which we compute $b \gets A x$ and then solve $A y = b$ to some precision. The latter is convenient because it lets us check accuracy. Solvers like unpreconditioned GMRES are optimal (over a subspace) in the 2-norm of the residual, thus the $A^* A$-norm of the error. This is a significantly weaker norm than the 2-norm of the error, for example, thus it really is useful to compare solvers in the error norm.
Dangers of choosing a random vector
A common method of generating $x$ or $b$ is to draw from an independent random distribution. This is flawed because these random vectors will be nearly orthogonal to the eigenvectors associated with the smallest eigenvalues -- exactly the modes that slow solver convergence. For example, if you use a Gaussian distribution with mean zero for a Neumann Laplacian (has a constant null space), then the average over any significant portion of the domain will be nearly zero, meaning that no long-range interactions take place. In this case, a coarse grid becomes unnecessary since you only have to move information far enough for the random process to homogenize.
If your distribution is orthogonal to that lowest eigenvector, every random vector ($x$ or $b$ doesn't matter when comparing to an eigenvector $z$, though if $z^T x$ is close to zero, $z^T b = z^T A x$ will be even smaller) that you draw will have this property of being nearly orthogonal to the "problematic modes", thus converging artificially fast. An independent random distribution will be orthogonal to most or all these "problematic modes" (e.g., translations and rotations for elasticity).
This mistake gets published from time to time, so it's worth increasing awareness.
If you don't know anything about the problem, it's inconvenient to tune the random distribution so that its expected value has nonzero component in the direction of all the low-energy eigenvectors---you would have to solve an expensive eigenproblem to find those "bad" vectors. Instead, I would recommend the following ad-hoc procedure:
- Recursively bisect the graph.
- At each level of the bisection tree, choose a random value from a random distribution with nonzero mean.
- Define the vector $x$ to be the sum of the random value assigned to all the parts that contain it (i.e., the sum of the path from the root of the tree).
This is much like random sampling in a wavelet basis and provides correlation at all scales. It doesn't require any special knowledge about the SPD matrix other than its sparsity pattern.
Extensions and alternatives
[Collecting comments from offline discussion.]
Michael Grant asks "What about dense systems?" My reply: For an unstructured dense matrix, I would threshold since that should be some measure of locality. For a structured dense matrix (e.g., an $H$-matrix), you already have a hierarchy. The details of the bisection/partitioning isn't that important. All it needs to do is provide correlation at different scales so that the distribution isn't orthogonal to the problematic eigenvectors.
Jack Poulson suggests using a combination of "point-sources, plane waves, wave packets, and Gaussian random vectors." My reply: Plane waves and wave packets require some additional knowledge beyond the matrix itself. Point sources aren't great for heterogeneous media problems because the Green's functions might be local almost everywhere, but global in very select places (e.g., faults, wires). The chance of a random point source activating that long-range coupling is low, but it's important for the real system. The bulk of the answer above was explaining why vectors drawn from an independent random distribution are inadequate.