A finite element discretization of elliptic or parabolic equations leads to large sparse linear equation systems \( Ax=b \). To solve such a system, a direct solver may be a choice among others (see section Solvers). Kaskade 7 provides interfaces to the direct solver libraries UMFPACK, PARDISO, MUMPS, SUPERLU, UMFPACK3264, UMFPACK64.
| solver library | description |
|---|---|
| MUMPS | MUltifrontal Massively Parallel Solver: MUMPS is a package for solving systems of linear equations of the form Ax = b, where A is a square sparse matrix and can be either unsymmetric, symmetric positive definite, or general symmetric. MUMPS implements a direct method based on a multifrontal approach which performs a Gaussian factorization A = LU where L is a lower triangular matrix and U an upper triangular matrix. If the matrix is symmetric then the factorization A = LDLT where D is block diagonal matrix with blocks of order 1 or 2 on the diagonal is performed. |
| UMFPACK | Unsymmetric MultiFrontal Package: UMFPACK is a set of routines for solving unsymmetric sparse linear systems Ax = b using the Unsymmetric MultiFrontal method. It uses dynamic memory allocation, and has a symbolic preordering and analysis phase that also reports the upper bounds on the nonzeros in L and U, flop count, and memory usage in the numeric phase. It can be used for real and complex matrices, rectangular and square, and both non-singular and singular. It is included in SuiteSparse, a suite of sparse matrix algorithms. UMFPACK comes in three versions: The default one uses 32 bit integers for indices and thus cannot treat very large systems. The 64 bit version removes this limitation at the cost of larger memory consumption and slightly lower performance. |
| SUPERLU | Supernodal LU: SuperLU is a general purpose library for the direct solution of large, sparse, nonsymmetric systems of linear equations. The library routines performs an LU decomposition with partial pivoting and triangular system solves through forward and back substitution. The LU factorization routines can handle non-square matrices but the triangular solves are performed only for square matrices. |
| PARDISO | PARallel DIrect SOlver: The package PARDISO is a thread-safe, high-performance, robust, memory efficient and easy to use software for solving large sparse unsymmetric, structurally symmetric or symmetric, real or complex, positive definite or indefinite, hermitian linear systems of equations on shared-memory and distributed-memory multiprocessors. PARDISO implements an LU decomposition with complete pivoting and fast inversion strategies. |
Some methods in Kaskade require the specification of direct solver type and matrix properties. These are defined in Kaskade as enum in utilities/enums.hh:
enum class DirectType { ANY=-1, UMFPACK, PARDISO, MUMPS, SUPERLU, UMFPACK3264, UMFPACK64 };
enum class MatrixProperties { GENERAL, SYMMETRICSTRUCTURE, SYMMETRIC, POSITIVEDEFINITE };
Note that matrix properties may be set for the use of MUMPS and PARDISO, while this is not necessary for UMFPACK and SUPERLU.
The easiest (but not necessarily fastest in terms of computation time) way to solve a linear system is by using the function directInverseOperator. Its usage is simple and straightforward. With an operator A from the assembler, a coefficient vector x made of zeros for the unknown solution and a right hand side b, the inverse operation can be performed. The linear system can be solved by linear algebra using the function applyScaleAdd(alpha,b,x), calculating \( x + \alpha A^{-1} b \) in our case:
directInverseOperator(A,DirectType::UMFPACK,MatrixProperties::GENERAL).applyscaleadd(alpha,b,x);
This is a convenience function which constructs an inverse linear operator and relies on the Dune::LinearOperator interface, defining a linear operator \( A: X \rightarrow Y \). Two operator types for A which hold the assembled matrices are possible to use as input: a Galerkin operator AssembledGalerkinOperator which can be created directy from the assembler or a MatrixRepresentedOperator, where the user can decide which matrix type to use (BCRS matrix, matrix in triplet format, …). The use of the latter could be necessary if the assembled matrix had to be extracted from the assembler for further processing. Choose a direct solver type and set sparse matrix properties if known. These properties characterize the symmetry or positive definiteness of the sparse matrix. Default setting if nothing else is specified by the user is UMFPACK for matrix factorization and general symmetry.
Factorization of a matrix is the decomposition into a product of matrices such as LU factorization of A \( A=LU \) with L lower and U upper triangular matrix, Cholesky factorization af A \( A=LL^T \) with L lower triangular matrix and so on. Various numerical methods exist to compute such decompositions and are implemented in the libraries used by Kaskade (see above Direct solver libraries). Besides using the directInverseOperator (see The direct inverse operator), the factorization can be created explicitely and used for the solution of the linear system.
For a matrix A and a right hand side rhs, resulting from e.g. the finite element assembly of a problem, the factorization with UMFPACK and the solution of the system for x can be implemented as follows:
#include "linalg/umfpack_solve.hh"
[...]
triplet = A.template get<MatrixAsTriplet<double> >(); // get matrix A in triplet format
A.getAssembler().toSequence(0,neq,rhs.begin()); // fill rhs vector with values from assembler
Factorization<double> *matrix = 0;
matrix = new UMFFactorization<double>(size,triplet.ridx,triplet.cidx,triplet.data); // in case of using UMFPACK
matrix->solve(rhs,x);
delete matrix;
Factorization is an abstract base class used for matrix factorization. The matrix is defined by the vectors ridx, cidx and data containing the row and column index and the value of the matrix entries. Method solve then solves the system \( Ax=b \) for the given right hand side \( b \), which will be overwritten with the solution \( x \). The parameter size describes the total numbers of degrees of freedom in the variables and can be retrieved e.g. from the ansatz variables.
The same holds for factorization with MUMPS
#include "linalg/mumps_solve.hh"
[...]
matrix = new MUMPSFactorization<double>(size,triplet.ridx,triplet.cidx,triplet.data);
as well as SUPERLU
#include "linalg/superlu_solve.hh"
[...]
matrix = new SUPERLUFactorization<double>(size,triplet.ridx,triplet.cidx,triplet.data);
and PARDISO
#include "linalg/pardiso_solve.hh"
[...]
matrix = new PARDISOFactorization<double>(size,triplet.ridx,triplet.cidx,triplet.data);