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