Two specialized drivers, run_taucs.cpp and run_mumps.cpp presented previously show a common problem when one needs to deal with a sparse solver. If you expect that in the future you may need to change a sparse solver, you must be then ready to modify the part of your code that works with the sparse solver. For example all algorithms that use the sparse solver must be modified if you decide to change a linear solver.

A solution would be to generalize function calls necessary for the application and then to create an application code in terms of such calls only. Then your code will be solver independent. In C++ this is usually achieved by using virtual functions and below I present such a general interface. Please note that there are libraries (Trilinos and PETSc) where you will find a more involved interface with more functionality. Yet on the other side they are more complex and it will take definitely more time to learn them.

Let us consider class `LinearSolver`

in solvers.h. The function there are pure virtual functions and that means that this class cannot be instantiated:

`LinearSolver sol; // compiled time error`

The goal of the class is to present a generalized interface that will be implemented in the derived classes. Let us consider it function by function

` virtual void* setMatrix(SparseMatrix &mat) = 0;`

The function takes a `SparseMatrix`

and converts it to the internal representation required by the solver. The memory of the `SparseMatrix`

object will be freed during the call. We obtain a handle to the sparse matrix in the form of the pointer to void and we can use it later on to refer to that matrix.

` virtual void setVerbose() = 0;`

The function forces the solver to write internal information to the standard output if available.

` virtual void mulMatrixByVec(void *mat, double *in, double *out) = 0;`

The matrix is multiplied by a vector `in`

and the result will be written in the vector `out`

. The memory for `out`

must be allocated by the user before the function call.

` virtual void* factor(void *mat, bool ClearMatrix, const string ¶m) = 0;`

virtual void mulInverseByVec(void *L, double *in, double *out) = 0;

The main functions for linear solve, first the factorization, second the back substitution. For the factorization function we also specify if the matrix to factor should be left intact or removed from the memory as well as an extra parameter to the solver (for example an ordering). The function `factor`

returns a handle to the factor that will be used in the back substitution function.

` virtual void* copyMatrix(void *mat1) = 0;`

virtual void sumMatrices(void *mat1, void* mat2, const double &alpha) = 0;

Functions to copy a matrix if necessary and to find a sum of two matrices, `mat1 = mat1 + alpha * mat2`

.

` virtual bool isSymmetric(void *mat) = 0;`

virtual bool supportSymmetric() = 0;

Functions to query if the matrix is symmetric and if the solver supports symmetric matrices. For example UMFPACK is working only with full matrices. If one uses the matrix in the symmetric storage, then the matrix will be considered by UMFPACK as a triangular matrix and the results will be wrong.

` virtual void clearMatrix(void *mat) = 0;`

virtual void clearFactor(void *L) = 0;

virtual void clear() {}

Functions to clear memory for the matrix `mat`

, for the factor `L`

, and finally to clear all the temporary structures if they have been created.

Now suppose that two classed derived from LinearSolver have been implemented, TAUCS and MUMPS. Then it is possible to write

`LinearSolver *one = new TAUCS;`

LinearSolver *two = new MUMPS;

...

one->factor(mat1, true, "metis");

two->factor(mat2, true, "metis");

At that point object `one`

will use function `factor`

implemented in the TAUCS class and object `two `

will employ function `factor`

implemented in the MUMPS class. That shows you the power of using virtual functions.

In order to simplify the creation of `LinearSover`

, there is a static function

` static LinearSolver* create(const string &str);`

It is however not a part of the generalized interface and its implementation will be considered later after we implement classes for TAUCS, MUMPS, UMFPACK and PARDISO.

## Previous

Introduction

Class Matrix

Class SparseMatrix

Sample code in C++ to run TAUCS

Sample code in C++ to run MUMPS

## Next

Class UMFPACK

Class TAUCS

Class MUMPS and PARDISO

Gluing all together: LinearSolver::create

Sample driver to run a sparse solver

RSS