7Oct/11Off

Sparse 101: Calling a Sparse System Solve

by Dan

In this post, I’ll show you how to call a matrix solve using our soon to be released CULA Sparse package. For this example, I will show the main tasks to calling our Cg solver with a Block Jacobi preconditioner for a CSR matrix with double-precision data.

The first set of parameters to consider is the matrix system to be solved, Ax=b. For these inputs, you need to consider which matrix formats and which precision you are using; see this blog post for a discussion of matrix formats. The relevant parameters for this system are:

• n - size of the matrix system
• nnz - number of non-zero elements
• A, colInd, rowPtr - a CSR representation of the matrix to be solved
• x - the solution vector
• b - the right-hand-side to solve against

These parameters will be passed directly to a function with “Dcsr” in its name to denote the double-precision data (D) and CSR representation (csr), such as in the line below:

`culaDcsr{Function}(..., N, nnz, A, colInd, rowPtr, x, b, ...);`

Now that I’ve discussed the matrix system, the next parameter to consider is the configuration structure for setting options that are common to all the solvers. Among these options are the solution relative tolerance, maximum number of iterations, maximum runtime, indexing-format, and whether debugging mode has been enabled. An example configuration may look like:

```culaIterativeConfig config;
culaIterativeConfigInit(&config);
config.maxIterations = 1000;
config.tolerance = 1e-6;
config.indexing = 1;```

As you can see, setting up problem options is an easy task. Each option is clear and self-documenting. The config initialization routine will set sensible defaults for most problems, but it’s worth double checking to see if they meet up with your own needs and perhaps overriding them as we have done above.

The last set of options to consider are the solver- and preconditioner-specific options. These options are done with a set of structures that are initialized similarly to the general configuration structure. To use the Cg solver with a Blockjacobi preconditioner, you would write:

```culaCgOptions solverOptions;
culaCgOptionsInit(&solverOptions);

culaBlockjacobiOptions preOptions;
culaBlockjacobiOptionsInit(&preOptions);
preOptions.blockSize = 4;```

Above, we have default initialized both structures but then overrided the Block Jacobi preconditioner block size. Because each solver and preconditioner is initialized very similarly, this makes trying out different solvers and preconditioners an easy task.

Putting it all of the parameters together, we end up with the following line:

```culaDcsrCgBlockjacobi(&config, &solverOptions, &preOptions,
N, nnz, A, colInd, rowPtr, x, b, &result);```

That’s it! We built CULA Sparse so that it would be easy to set up and work with different options while making sure that the code is always clear about what it is doing.