Any book on the solution of nonsingular systems of equations is bound to start with Ax= J, but here, A is assumed to be symmetric. These systems arise frequently in scientific computing, for example, from the discretization by finite differences or by finite elements of partial differential equations. Usually, the resulting coefficient matrix A is large, but sparse. In many cases, the need to store the matrix factors rules out the application of direct solvers, such as Gaussian elimination in which case the only alternative is to use iterative methods. A natural way to exploit the sparsity...