Using Umfpack#

Besides some built-in direct solvers for sparse linear system, XLiFE++ provides an interface to Umfpack, a set of routines for solving non-symmetric sparse linear systems, \(Ax=b\) , using the Unsymmetric MultiFrontal method. More information can be found at UmfPack.

The aim of umfpackSupport is to make ease of the Umfpack for XLiFE++’s users, allowing them to solve the linear system \(Ax=b\) with a simple call. Instead of invoking “multi-argument” Umfpack solver, user can make call of them as any built-in iterative solver in the context of XLiFE++.

The Umfpack wrappers#

Written in ANSI/ISO C, UMFPACK library consists of 32 user-callable routines. Nearly all these routines come in four versions, with different sizes of integers and for real or complex floating-point numbers:

real double precision, int integers
complex double precision, int integers
real double precision, SuiteSparse_long integers
complex double precision, SuiteSparse_long integers

Only the interfaces to the first two versions are implemented in XLiFE++. However, more versions can be easily added in the need of user.

Because callable UMFPACK routines are written in C, to call them correctly from the C++ environment of XLiFE++, we must make these routines “extern C”. And it’s the role of UmfPackWrappers. Only some primary routines of UMFPACK are made interface; however, it’s not difficult to include others.

#ifdef __cplusplus
extern "C" {
#endif

int DISYMBOLIC(int n_row, int n_col, const int Ap[], const int Ai[], const double Ax[], void** Symbolic, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
int DINUMERIC(const int Ap[], const int Ai[], const double Ax[], void* Symbolic, void** Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
int DISOLVE(int sys, const int Ap[], const int Ai[], const double Ax[], double X[], const double B[], void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
void DIFREESYMBOLIC(void** Symbolic);
void DIFREENUMERIC(void** Numeric);
  ...
int ZISYMBOLIC(int n_row, int n_col, const int Ap[], const int Ai[], const double Ax[], const double Az[], void** Symbolic, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
int ZINUMERIC(const int Ap[], const int Ai[], const double Ax[], const double Az[], void* Symbolic, void** Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
int ZISOLVE(int sys, const int Ap[], const int Ai[], const double Ax[], const double Az[], double Xx[], double Xz[], const double Bx[], const double Bz[], void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
void ZIFREESYMBOLIC(void** Symbolic);
void ZIFREENUMERIC(void** Numeric);
  ...

Routine names beginning with DI correspond to the version of real double precision, int integers, the ones with ZI imply the versions of complex double precision, int integers.

The `Umfpack` class#

The UMFPACK class plays a role of high-level wrapper of routines provided by UmfPackWrappers:

SYMBOLIC: performing a symbolic decomposition on the sparsity of a matrix
NUMERIC: performing a numeric decomposition of a matrix
FREESYMBOLIC: freeing symbolic object created by SYMBOLIC
FREENUMERIC: freeing numeric object created by NUMERIC
GETLUNZ: returning the number of nonzeros in lower matrix L and upper matrix U in LU decomposition
GETNUMERIC: retrieving lower matrix L, upper matrix U, permutation P, Q, and scale factor R
GETDET: returning determinant of original matrix

The “*” corresponds to DI-real double precision, or ZI-complex double precision.

All these routines are organized under templated “class-like” interface.

template<typename OrdinalType, typename ScalarType>
class UMFPACK
{
  public:
    typedef typename NumTraits<ScalarType>::magnitudeType MagnitudeType;

    //! Default Constructor.
    inline UMFPACK(void) {}

    //! Destructor.
    inline ~UMFPACK(void) {}

    OrdinalType symbolic(OrdinalType nRow, OrdinalType nCol, const OrdinalType Ap[], const OrdinalType Ai[], const ScalarType Ax[], void** Symbolic, const MagnitudeType Control[UMFPACK_CONTROL], MagnitudeType Info[UMFPACK_INFO]);

    OrdinalType numeric(const OrdinalType Ap[], const OrdinalType Ai[], const ScalarType Ax[], void* Symbolic, void** Numeric, const MagnitudeType Control[UMFPACK_CONTROL], MagnitudeType Info[UMFPACK_INFO]);

    void freeSymbolic(void** Symbolic);
    void freeNumeric(void** Numeric);

    OrdinalType solve(OrdinalType sys, const OrdinalType Ap[], const OrdinalType Ai[], const ScalarType Ax[], ScalarType X[], const ScalarType B[], void* Numeric, const MagnitudeType Control[UMFPACK_CONTROL], MagnitudeType Info[UMFPACK_INFO]);

    OrdinalType getLunz(OrdinalType* lnz, OrdinalType* unz, OrdinalType* nRow, OrdinalType* nCol, OrdinalType* nzUdiag, void* Numeric);

    OrdinalType getNumeric(OrdinalType Lp[], OrdinalType Lj[], ScalarType Lx[], OrdinalType Up[], OrdinalType Ui[], ScalarType Ux[], OrdinalType P[], OrdinalType Q[], ScalarType Dx[], OrdinalType* doRecip, ScalarType Rs[], void* Numeric);

    OrdinalType getDeterminant(ScalarType* Mx, ScalarType* Ex, void* NumericHandle, ScalarType UserInfo[UMFPACK_INFO]);
   ...

One advantage of this organization is to allow the specialization of code corresponding to the two versions: real double precision and complex double precision with int integer.

template<>
class UMFPACK<int, double>
{
  public:
    inline UMFPACK(void) {}
    inline ~UMFPACK(void) {}

    int symbolic(int nRow, int nCol, const int Ap[], const int Ai[], const double Ax[], void** Symbolic, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);

    int numeric(const int Ap[], const int Ai[], const double Ax[], void* Symbolic, void** Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);

    void freeSymbolic(void** Symbolic);

    void freeNumeric(void** Numeric);

    int solve(int sys, const int Ap[], const int Ai[], const double Ax[], double X[], const double B[], void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
   ...
}

The specialization of complex version comes with an extra method for the linear system \(Ax=b\) with \(A\) complex double precision and \(b\) real double precision.

template<>
class UMFPACK<int, std::complex<double> >
{
  public:
    inline UMFPACK(void) {}
    inline ~UMFPACK(void) {}

    int symbolic(int nRow, int nCol, const int Ap[], const int Ai[], const std::complex<double> Ax[], void** Symbolic, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);

    int numeric(const int Ap[], const int Ai[], const std::complex<double> Ax[], void* Symbolic, void** Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);

    void freeSymbolic(void** Symbolic);

    void freeNumeric(void** Numeric);

    int solve(int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], std::complex<double> X[], const std::complex<double> B[], void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);

    // Solve in case AX=B with A complex, B real
    int solveCR(int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], std::complex<double> X[], const double Bx[], const double Bz[], void* Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]);
    ...

The `UmfPackSolver` class#

The class, as its name, serves for solving the sparse linear system \(Ax=b\) . Being similar to other solvers of XLiFE++, this solver can be provoked with a call of operator(). For example:

UmfPackSolver umfpackSolver;
// With A largeMatrix; b and x are std::vector
umfpackSolver(A,b,x);

By calling the solver in this way, the storage and values of the largeMatrix \(A\) are extracted and these extractions are converted to the format of Umfpack including: column pointer, row index and values. These three are inputs to underlying UMFPACK routines. Because extraction and conversion between largeMatrix of XLiFE++ and Umfpack format are expensive, in some cases, it’s better to directly make use of the conversion and it can be easily done by calling operator() with more inputs

// Solve AX=B with UmfPack
// \param[in] colPointer column pointer of CSC-like Umfpack format (after calling function extract2UmfPack of largeMatrix)
// \param[in] rowIdx row index of CSC-like Umfpack format (after calling function extract2UmfPack of largeMatrix)
// \param[in] values values of CSC-like Umfpack format (after calling function extract2UmfPack of largeMatrix)
// \param[in] vecB vector B (vector B should be std::vector)
// \param[in] vecX vector X (vector X should be std::vector)
template<class ScalarTypeMat, class VecB, class VecX>
void operator()(int size, const std::vector<int>& colPointer, const std::vector<int>& rowIdx, const std::vector<ScalarTypeMat>& values, const VecB& vecB, VecX& vecX, UmfPackComputationMode sys=_A)

One of a technical issue on working with different types of scalar value, in this case, there are Real and Complex values, is the static dispatching. For dispatching at runtime, we can use simple if-else statements or the switch statement. However, the if-else statement requires both branches to compile successfully, even when the condition tested by if is known at compile time. To overcome this problem, we use a simple technique that is initially described in Alexandrescu, mapping integral Constants to Types:

template <int v>
struct Int2Type
{
  enum { value = v };
};

Int2Type generates a distinct type for each distinct constant integral value passed. This is because different template instantiations are distinct types; thus, Int2Type<0> is different from Int2Type<1>, and so on. In addition, the value that generates the type is “saved” in the enum member value.

For a linear system \(Ax=b\), we can divide into four cases depending on the scalar type of \(A\) and \(b\) , and corresponding to each case, there is a Int2Type:

\(A\) real and \(b\) real; Int2Type<0>
\(A\) real and \(b\) complex; Int2Type<1>
\(A\) complex and \(b\) real; Int2Type<2>
\(A\) complex and \(b\) complex; Int2Type<3>

Because each case above needs processing differently, we define an internal class of UmfPackSolver to take care of it.

// Auxiliary internal class to process 4 cases of equation AX = B:
//      + A:real, B: real
//      + A:real, B:complex
//      + A:complex, B: real
//      + A:complex, B:complex
template<typename ScalarTypeMat, typename VecB, typename VecX>
class SolverIntern
{
  public:
    typedef typename NumTraits<ScalarTypeMat>::magnitudeType MagType;
    typedef typename Conditional<NumTraits<ScalarTypeMat>::IsComplex, ScalarTypeMat, typename VectorTrait<VecB>::Type>::type ScalarType;

    SolverIntern() {}
   ...

For example, the case of \(A\) real and \(b\) real.

// Solve AX=B in case A Real, B Real
void solve(int sys, int nRows, int nCols, const std::vector<int>& colPointer, const std::vector<int>& rowIdx, const std::vector<ScalarTypeMat>& values, const VecB& vecB, VecX& vecX, Int2Type<0>)
{
  void* symbolic;
  void* numeric;
  double* null = (double*) NULL;

  UMFPACK<int, ScalarType> umfPack;

  (void)umfPack.symbolic(nRows, nCols, &(colPointer[0]), &(rowIdx[0]),&(values[0]),&symbolic,null,null);
  (void)umfPack.numeric(&(colPointer[0]), &(rowIdx[0]), &(values[0]), symbolic, &numeric, null,null);
  (void)umfPack.solve((int)sys, &(colPointer[0]), &(rowIdx[0]),&(values[0]), &vecX[0], &vecB[0], numeric, null, null);
  umfPack.freeSymbolic(&symbolic);
  umfPack.freeNumeric(&numeric);
}

The `UmfPackLU` class#

UMFPACK library allows not only to solve the linear system \(Ax=b\) but also to retrieve more information. There are cases where users may wish to do more with the LU factorization of a matrix rather than solve a linear system. And the UmfPackLU class is for this purpose. Not like to the “non-templated” UmfPackSolver class, the UmfPackLU depends on the type of input matrix.

template<typename MatrixType>
class UmfPackLU
{
  public:
    typedef typename MatrixType::ScalarType ScalarType;

    struct LUMatrixType {
      std::vector<int> outerIdx;
      std::vector<int> innerIdx;
      std::vector<ScalarType> values;
     };

  public:
    UmfPackLU() { init();}
    UmfPackLU(const MatrixType& matrix)
    {
      init();
      compute(matrix);
    }
 ...

Apart from solving the linear system \(Ax=b\) , this class provides some methods to work with LU factorization.

inline const LUMatrixType& matrixL() const
{
  if (extractedDataAreDirty_) extractData();
  return (lMatrix_);
}

inline const LUMatrixType& matrixU() const
{
  if (extractedDataAreDirty_) extractData();
  return (uMatrix_);
}

inline const std::vector<int>& permutationP() const
{
  if (extractedDataAreDirty_) extractData();
  return (pPerm_);
}

inline const std::vector<int>& permutationQ() const
{
  if (extractedDataAreDirty_) extractData();
  return (qPerm_);
}

inline const std::vector<Real>& scaleFactorR() const
{
  if (extractedDataAreDirty_) extractData();
  return (rScaleFact_);
}

All these methods above return the result of LU factorization. \(PAQ=LU\). The matrix \(U\) is returned in compressed column form (with sorted columns). The matrix \(L\) is returned in compressed row form (with sorted rows). One remarkable thing is that the compressed column (row) form is similar to CSC (CSR) of XLiFE++. These two matrix are returned in form of struct LUMatrixType

struct LUMatrixType
{
  std::vector<int> outerIdx;
  std::vector<int> innerIdx;
  std::vector<ScalarType> values;
};

The outerIdx corresponds to the colPointer of CSC and rowPointer of CSR, while the innerIdx corresponds to the rowIndex of CSC and colIndex of CSR. Remember that, not like CSC or CSR, the values of LUMatrixType struct doesn’t contain the “first-position” zero (0).

The permutations P and Q are represented as permutation vectors, where \(P[k] = i\) means that row i of the original matrix is the k-th row of \(PAQ\) , and where \(Q[k] = j\) means that column \(j\) of the original matrix is the k-th column of PAQ.

The vector R is the scale factor in the LU factorization. The first value of scaleFactorR defines how the scale factors Rs are to be interpreted. If it’s one(1), then the scale factors scaleFactorR[i+1] are to be used by multiplying row i of matrix A by scaleFactorR[i]. Otherwise, the entries in row i are to be divided by scaleFactorR[i+1].

Before retrieving all these above values, a matrix at least needs factorizing. It can be done in two ways, by constructor

UmfPackLU(const MatrixType& matrix)
{
  init();
  compute(matrix);
}

or by the method compute

/*! Computes the sparse Cholesky decomposition of \a largeMatrix */
void compute(const MatrixType& matrix)
{
  analyzePattern(matrix);
  factorize(matrix);
}

In this function, a symbolic decomposition on the sparsity of a matrix is performed by analyzePattern then a numeric decomposition performed by factorize. After the matrix is factorized, its factorization can be used to solve the linear system \(Ax=b\) with :

template<typename VectorScalarType>
void solve(const std::vector<VectorScalarType>& b, std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, VectorScalarType>::type >& x) const;

One advantage of UmfPackLU over UmfPackSolver is that for a same matrix A, it only needs to do factorization once then all these factorizations can be used over and over again to solve the linear system \(Ax=b\) . Of course, it doesn’t have the flexibility of “non-templated” UmfPackSolver, it needs template arguments in order to define type of matrix. Depending on a specific demand, one can be a better choice than another!!!

Using Umfpack#

The Umfpack wrappers#

The Umfpack class#

The UmfPackSolver class#

The UmfPackLU class#

The `Umfpack` class#

The `UmfPackSolver` class#

The `UmfPackLU` class#