inputs/dev/UmfpackSupport.tex

%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%

\section{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:
 \begin{itemize}
\item real double precision, \textbf{int} integers
\item complex double precision, \textbf{int} integers
\item real double precision, \textbf{SuiteSparse\_long} integers
\item complex double precision, \textbf{SuiteSparse\_long} integers
\end{itemize}
Only the interfaces to the first two versions are implemented in \xlifepp. However, more versions can be easily added in the need of user. \\
Because callable UMFPACK routines are written in C, so as to call them  correctly from the C++ environment of \xlifepp,  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.
\vspace{.1cm}
\begin{lstlisting}
#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);
    ...
\end{lstlisting}
\vspace{.1cm}
Routine names beginning with DI correspond to the version of real double precision, \textbf{int} integers, the ones with ZI imply the versions of
complex double precision, \textbf{int} integers.

\section{The {\classtitle UMFPACK} class}

The {\class UMFPACK} class plays a role of high-level wrapper of routines provided by UmfPackWrappers:
\begin{itemize}
\item *SYMBOLIC: performing a symbolic decomposition on the sparcity of a matrix
\item *NUMERIC: performing a numeric decomposition of a matrix
\item *FREESYMBOLIC: freeing symbolic object created by *SYMBOLIC
\item *FREENUMERIC: freeing numeric object created by *NUMERIC
\item *GETLUNZ: returning the number of non-zeros in lower matrix L and upper matrix U in LU decomposition
\item *GETNUMERIC: retrieving lower matrix L, upper matrix U, permutation P,Q, and scale factor R
\item *GETDET: returning determinant of original matrix
\end{itemize}
The "*" corresponds to DI-real double precision, or ZI-complex double precision.\\
All these routines are organized under templated "class-like" interface.
\vspace{.1cm}
\begin{lstlisting}
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]);
    ...
\end{lstlisting}
\vspace{.1cm}
 One advantage of this orgnization is to allow the specialization of code corresponding to the two versions: real double precision and complex double precision with \textbf{int} integer.
\vspace{.1cm}
\begin{lstlisting}
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]);
    ...
\end{lstlisting}
\vspace{.1cm}
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.
\vspace{.1cm}
\begin{lstlisting}
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]);
    ...
\end{lstlisting}
\vspace{.1cm}
\displayInfos{library=umfpackSupport, header=UmfPack.hpp, implementation=UmfPack.cpp, test={test\_UmfPackSolver.cpp}, header dep={utils.h}}

\section{The {\classtitle UmfPackSolver} class}

The class, as its name, serves for solving the sparse linear system $Ax=b$. Being similar to other solvers of \xlifepp, this solver can be provoked with a call of operator(). For example:
\vspace{.1cm}
\begin{lstlisting}
UmfPackSolver umfpackSolver;
// With A largeMatrix; b and x are std::vector
umfpackSolver(A,b,x);
\end{lstlisting}
\vspace{.1cm}
By calling the solver in this way, the storage and values of the {\class 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 \xlifepp 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
 \vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] size, colPointer}]
// 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)
\end{lstlisting}
\vspace{.1cm}
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 {\itshape if-else} statements or the switch statement. However, the {\itshape 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 intially described in Alexandrescu, mapping integral Constants to Types:
 \vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] value}]
template <int v>
struct Int2Type
{
enum { value = v };
};
\end{lstlisting}
\vspace{.1cm}
{\class Int2Type} generates a distinct type for each distinct constant integral value passed. This is because
different template instantiations are distinct types; thus, {\class Int2Type<0>} is different from {\class 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 to the scalar type of $A$ and $b$, and corresponding to each case, there is a {\class Int2Type}:
\begin{itemize}
\item $A$ real and $b$ real; {\class Int2Type<0>}
\item $A$ real and $b$ complex; {\class Int2Type<1>}
\item $A$ complex and $b$ real; {\class Int2Type<2>}
\item $A$ complex and $b$ complex; {\class Int2Type<3>}
\end{itemize}
Because each case above needs processing in a different way, we define an internal class of {\class UmfPackSolver} to take care of it.
 \vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] type}]
// 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() {}
   ...

\end{lstlisting}
\vspace{.1cm}
For example, the case of $A$ real and $b$ real.
 \vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] colPointer, symbolic, numeric}]
// 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);
}
\end{lstlisting}
\vspace{.1cm}
\displayInfos{library=umfpackSupport, header=UmfPackSolver.hpp, test={test\_UmfPackSolver.cpp}, header dep={utils.h, UmfPack.hpp}}

\section{The {\classtitle 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 {\class UmfPackLU} class is for this purpose. Not like to the "non-templated" {\class UmfPackSolver} class, the {\class UmfPackLU} depends on the type of input matrix.
 \vspace{.1cm}
\begin{lstlisting}
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);
    }
...
\end{lstlisting}
\vspace{.1cm}
Apart from solving the linear system $Ax=b$, this class provides some methods to work with LU factorization.
 \vspace{.1cm}
\begin{lstlisting}
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_);
}
\end{lstlisting}
\vspace{.1cm}
All these methods above return the result of LU factorization. $PAQ=LU$.  The matrix {\bfseries U} is returned in compressed column form (with sorted columns). The matrix {\bfseries 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 \xlifepp. These two matrix are returned in form of struct {\class LUMatrixType}
 \vspace{.1cm}
\begin{lstlisting}
struct LUMatrixType
{
  std::vector<int> outerIdx;
  std::vector<int> innerIdx;
  std::vector<ScalarType> values;
};
\end{lstlisting}
\vspace{.1cm}
The {\itshape outerIdx} corresponds to the {\itshape colPointer} of CSC and {\itshape rowPointer} of CSR, while the {\itshape innerIdx} corresponds to the {\itshape rowIndex} of CSC and {\itshape colIndex} of CSR. Remember that, not like CSC or CSR, the {\itshape values} of {\class 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 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 {\itshape scaleFactorR} defines how the scale factors Rs are to be interpreted. If it's one(1), then the scale factors {\itshape scaleFactorR[i+1]} are to be used by multiplying row i of matrix A by {\itshape scaleFactorR[i]}. Otherwise, the entries in row i are to be divided by {\itshape scaleFactorR[i+1]}. \\
Before retrieving all these above values, a matrix at least needs factorizing. It can be done in two ways, by constructor
 \vspace{.1cm}
\begin{lstlisting}
UmfPackLU(const MatrixType& matrix)
{
  init();
  compute(matrix);
}
\end{lstlisting}
\vspace{.1cm}
or by the method {\bfseries compute}
 \vspace{.1cm}
\begin{lstlisting}
/*! Computes the sparse Cholesky decomposition of \a largeMatrix
 */
void compute(const MatrixType& matrix)
{
  analyzePattern(matrix);
  factorize(matrix);
}
\end{lstlisting}
\vspace{.1cm}
In this function, a symbolic decomposition on the sparsity of a matrix is performed by {\bfseries analyzePattern} then a numeric decomposition performed by {\bfseries factorize}. \\
After the matrix is factorized, its factorization can be used to solve the linear system $Ax=b$ with :
 \vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, type}]
template<typename VectorScalarType>
void solve(const std::vector<VectorScalarType>& b, std::vector<typename Conditional<NumTraits<ScalarType>::IsComplex, ScalarType, VectorScalarType>::type >& x) const;
\end{lstlisting}
\vspace{.1cm}
One advantage of {\class UmfPackLU} over {\class 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" {\class 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!!!
\displayInfos{library=umfpackSupport, header=UmfPackLU.hpp, test={test\_UmfPackSolver.cpp}, header dep={utils.h, UmfPack.hpp}}