inputs/dev/HMatrix.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{HMatrix}

Hierarchical matrix is a tree representation of a matrix, each node representing a sub-matrix being either split in sub-matrices (tree nodes) or not if the node is a leaf of the tree. The hierarchical division in sub-matrices is based on hierarchical division of row and column matrix indices or dofs if it is a FE matrix ({\class ClusterTree} objects).\\

To manage this tree representation, three template classes are provided:
\begin{itemize}
	\item the {\class HMatrixNode<T,I>} class dealing with a node of the matrix tree
	\item the {\class HMatrix<T,I>}  front end class
	\item the {\class HMatrixEntry<I>} interface class that handles some pointers to {\class HMatrix<real\_t,I>}, {\class HMatrix<complex\_t,I>}, ... .
\end{itemize}
These classes are templated by the type of matrix coefficients ({\tt T=real\_t, complex\_t}, ...) and the type of row and column clustering ({\tt I=Point, FeDof, Element}). Let us start by describing the {\class HMatrixNode} class supporting the main material.

\subsection{The {\classtitle HMatrixNode} class}

\vspace{.1cm}
\begin{lstlisting}
template <typename T, typename I>
class HMatrixNode
{public :
 HMatrixNode<T,I>* parent_;     //pointer to its parent, if 0 root node
 HMatrixNode<T,I>* child_;      //pointer to its first child, if 0 no child
 HMatrixNode<T,I>* next_;       //pointer to its brother, if 0 no brother
 ClusterNode<I>* rowNode_;      //row cluster node to access to row numbering
 ClusterNode<I>* colNode_;      //column cluster node to access to col numbering
 LargeMatrix<T>* mat_;          //pointer to a LargeMatrix
 ApproximateMatrix<T>* appmat_; //pointer to an ApproximateMatrix
 bool admissible_;              //block can be approximated
 number_t depth_;               //depth of node in tree, root has 0 depth
 number_t rowSub_, colSub_;     //sub block numbering (i,j)
 bool isDiag_;                  //true if matrix node is located on the diagonal
 \end{lstlisting}
 \vspace{.1cm}
 Each node refers to its parent node, its first child node, its first brother node and access to its row and column numbering using some pointers to {\class ClusterNode} object. If node is a leaf (no child), it stores the node sub-matrix either as a {\class LargeMatrix} (pointer {\tt mat\_}) or as a {\class ApproximateMatrix} (pointer {\tt appmat\_}). If node is not a leaf, both the pointers {\tt mat\_} and {\tt appmat\_} are 0. Besides, {\class HMatrixNode} manages some additional useful informations that are set during construction of the tree such as {\tt depth\_, rowSub\_, colSub\_, isDiag\_}. \\

 It manages also the particular data {\tt admissible\_} that tells from some geometrical rule if the sub-matrix attached to the node can be approximated, generally using compression methods. This property is mainly relevant in the context of BEM matrix where all interactions between distant dof sets can be reduced to few interactions. The choice of admissibility rule is governed by the {\class HMatrix} front end class using the enumerations:
\vspace{.1cm}
\begin{lstlisting}
enum HMatrixMethod {_standardHM, _denseHM};
enum HMAdmissibilityRule {_boxesRule};
\end{lstlisting}
\vspace{.1cm}
The {\class HMatrixNode} class has few constructors:
\vspace{.1cm}
\begin{lstlisting}
HMatrixNode();                                        //default constructor
HMatrixNode(HMatrixNode<T,I>*, HMatrixNode<T,I>*, HMatrixNode<T,I>*, number_t,
            ClusterNode<I>*, ClusterNode<I>*,  number_t, number_t,
            LargeMatrix<T>* , ApproximateMatrix<T>* );//full constructor
HMatrixNode(HMatrixNode<T,I>*, number_t);             //node constructor
HMatrixNode(LargeMatrix<T>*, HMatrixNode<T,I>*,
            number_t);                   //LargeMatrix leaf constructor
HMatrixNode(ApproximateMatrix<T>*, HMatrixNode<T,I>*,
            number_t);                   //ApproximateMatrix leaf constructor
HMatrixNode(const HMatrixNode<T,I>&);                 //copy constructor
HMatrixNode<T,I>& operator=(const HMatrixNode<T,I>&); //assign operator
void copy(const HMatrixNode<T,I>&);                   //copy function
\end{lstlisting}
\vspace{.1cm}
 and some stuff related to destructor or cleaner:
 \vspace{.1cm}
\begin{lstlisting}
~HMatrixNode();               //destructor
void clear();                 //clear all except ClusterNode pointers
void clearMatrices();         //clear (deallocate) only matrices
\end{lstlisting}
\vspace{.1cm}
\begin{remark}
	Be cautious when copying {\class HMatrixNode} object, if allocated, the large matrix or the approximate matrix are copied, but not the {\class ClusterNode} objects (shared pointers!). In the same way, when a {\class HMatrixNode} object is cleared or deleted, matrix or approximate matrix are deleted if they have been allocated but the {\class ClusterNode} objects are not deallocated.
\end{remark}\\

The most important function is
\vspace{.1cm}
\begin{lstlisting}
void divide(number_t rmin, number_t cmin, number_t maxdepth = 0,
            HMAdmissibilityRule=_boxesRule, bool sym=false);
\end{lstlisting}
\vspace{.1cm}
that creates recursively the sub-tree of the current node according to the admissibility rule given by the enumeration
\vspace{.1cm}
\begin{lstlisting}
enum HMAdmissibilityRule {_boxesRule};
\end{lstlisting}
\vspace{.1cm}
Up to now, only one rule is available, the standard boxes rule used in BEM computation telling that a {\class HMatrixNode} is admissible if the bounding box $B_r$ of the row cluster node and the bounding box $B_c$ of the column cluster node satisfy:
$$\text{diam}(B_r) \leq 2*\eta*\text{dist}(B_r,B_c).$$
\begin{remark}
	Note that in the case of a matrix having symmetry property (symmetric, skew-symmetric, self-adjoint or skew-adjoint), the division process may take it into account (argument {\tt sym=true}). In that case, it does not generate sub-matrix nodes that are located above the diagonal, so saving memory.
\end{remark}\\

Once the tree is built, some various member functions are available:
\vspace{.1cm}
\begin{lstlisting}
number_t numberOfRows() const;       //number of rows counted in T
number_t numberOfCols() const;       //number of columns counted in T
dimPair dimValues() const;           //dimensions of matrix values
void setClusterRow(ClusterNode<I>*); //update row cluster node pointers
void setClusterCol(ClusterNode<I>*); //update col cluster node pointers
number_t nbNonZero() const;          //number of T coefficients used
void getLeaves(list<HMatrixNode<T,I>* >&,
               bool = true) const;   //get leaves as a list
\end{lstlisting}
\vspace{.1cm}
The {\class HMatrixNode} class implements all the stuff required by algebraic operations on {\class HMatrix} class:
\vspace{.1cm}
\begin{lstlisting}
vector<T>& multMatrixVectorNode(const vector<T>&, vector<T>&,
                                SymType symt=_noSymmetry) const;
vector<T>& multMatrixVector(const vector<T>&, vector<T>&,
                            SymType symt=_noSymmetry) const;
real_t norm2() const;                             //Frobenius norm
real_t norminfty() const;                         //infinite norm
\end{lstlisting}
\vspace{.1cm}
and by printing operations
\vspace{.1cm}
\begin{lstlisting}
void printNode(ostream&) const;          //print current node contents
void print(ostream&) const;              //print tree from current node
void printStructure(ostream&, number_t, number_t,
                    bool all=false, bool shift=false) const;
\end{lstlisting}
\vspace{.1cm}

Finally, for performance measurement purpose, the class manages also some static data and static functions:
\vspace{.1cm}
\begin{lstlisting}
static bool counterOn;         //flag to activate the operations counter
static number_t counter;       //operations counter
static void initCounter(number_t n=0); //init the counter and enable it
static void stopCounter();             //disable counter
\end{lstlisting}
\vspace{.1cm}

\subsection{The {\classtitle HMatrix} class}

The {\class HMatrix} class is the front end class of hierarchical matrices. It is templated by the type of matrix coefficients and the type of cluster objects. It manages the parameters of the matrix clustering , some general informations and handles the root node of the tree:
 \vspace{.1cm}
 \begin{lstlisting}
template <typename T, typename I>
class HMatrix
{private :
 HMatrixNode<T,I>* root_;  //root node of the tree representation
 ClusterTree<I>* rowCT_;   //row cluster tree
 ClusterTree<I>* colCT_;   //col cluster tree
 public :
 string_t name;            //optional name, useful for doc
 ValueType valueType_;     //type of values (real, complex)
 StrucType strucType_;     //structure of values (scalar, vector, matrix)
 HMatrixMethod method_;    //method to define admissible block
 HMAdmissibilityRule admRule_; //block admissible rule
 real_t eta_;              //ratio in the admissibility criteria
 number_t rowmin_, colmin_;//minimum size of block matrix
 SymType sym_;             //type of symmetry
 number_t depth;           //maximal depth (info)
 number_t nbNodes;         //number of nodes (info)
 number_t nbLeaves;        //number of leaves (info)
 number_t nbAdmissibles;   //number of admissible blocks (info)
 number_t nbAppMatrices;   //number of approximate matrices (info)
 ...
 \end{lstlisting}
 \vspace{.1cm}
 The data member {\tt depth, nbNodes, nbLeaves, nbAdmissibles, nbAppMatrices}  are set by the division algorithm or computed on demand.\\

 The class redefines the default constructor (initializing data members) and an explicit constructor filling the cluster parameters and building the tree structure recursively:
 \vspace{.1cm}
 \begin{lstlisting}
HMatrix();
HMatrix(ClusterTree<I>&, ClusterTree<I>&,  number_t, number_t, number_t =0,
        const string_t& = "", SymType =_noSymmetry,
        HMatrixMethod =_standardHM, HMAdmissibilityRule=_boxesRule, real_t =1.);
void buildTree();                              //build the tree
void updateInfo();                             //update tree info (depth, ...)
HMatrix(const HMatrix<T,I>&);                  //copy constructor
HMatrix<T,I>& operator=(const HMatrix<T,I>&);  //assign operator
~HMatrix();                                    //destructor
void copy(const HMatrix<T,I>&);                //copy all
void clear();                                  //clear all
void clearMatrices();                          //deallocate matrix pointers
\end{lstlisting}
\vspace{.1cm}
 To preserve integrity of {\class HMatrix} object, the copy constructor and the assign operator do a hard copy of the tree and of the matrices attached to  {\class HMatrixNode} objects but the row and col cluster trees are not copied ! When a {\class HMatrix} object is deleted or cleared, all the nodes of the tree and node matrices are deleted but not the row and col cluster trees!\\

The {\class HMatrix} class provides some stuff to get or extract various informations about the tree structure:
\vspace{.1cm}
\begin{lstlisting}
number_t numberOfRows() const;         //number of rows counted in T
number_t numberOfCols() const;         //number of cols counted in T
dimPair dimValues() const;             //dimensions of values, (1,1) when scalar
number_t nbNonZero() const;            //number of coefficients used
const ClusterTree<I>* rowTree() const; //access to row ClusterTree pointer
const ClusterTree<I>* colTree() const; //access to col ClusterTree pointer
void setClusterRow(ClusterTree<I>*);   //change the row ClusterTree pointer
void setClusterCol(ClusterTree<I>*);   //change the col ClusterTree pointer
pair<number_t, number_t> averageSize() const; //average size of leaves
number_t averageRank() const;          //average rank of admissible leaves
list<HMatrixNode<T,I>* > getLeaves(bool = true) const;   //get all leaves
\end{lstlisting}
\vspace{.1cm}
It is possible to import a {\class LargeMatrix} in a {\class Hmatrix} only if the tree structure already exists and to export a {\class Hmatrix} to a {\class LargeMatrix} (dense storage).
\vspace{.1cm}
\begin{lstlisting}
void load(const LargeMatrix<T>&,HMApproximationMethod = _noHMApproximation);
LargeMatrix<T> toLargeMatrix(StorageType st=_dense, AccessType at=_row) const;
\end{lstlisting}
\vspace{.1cm}
Some algebraic operations on {\class HMatrix} are available:
\vspace{.1cm}
\begin{lstlisting}
vector<T>& multMatrixVector(const vector<T>&, vector<T>&) const;
vector<T>& multMatrixVectorOmp(const vector<T>&, vector<T>&) const;
void addFELargeMatrix(const LargeMatrix<T>&); //add FE LargeMatrix to current
real_t norm2() const;                         //Frobenius norm
real_t norminfty() const;                       //infinite norm
\end{lstlisting}
\vspace{.1cm}
Note that the {\tt multMatrixVector} function is recursive so it is not safe in multithread computation. This is the reason why there exists a non recursive function {\tt multMatrixVectorOmp} that implements omp pragma to parallelize the computation.\\

The {\class Hmatrix} class provides the following print stuff:
\vspace{.1cm}
\begin{lstlisting}
void print(std::ostream&) const;
void printSummary(std::ostream&) const;
void printStructure(std::ostream&,bool all=false, bool =false) const;
void saveStructureToFile(const string_t&) const;
ostream& operator<<(ostream& os, const HMatrix<X,J>& hm);
\end{lstlisting}
\vspace{.2cm}
As an illustration, we show on figure \ref{fig_hmatrix} the structure of a HMatrix based on the clustering of a mesh of a sphere, get by the following \xlifepp code
\vspace{.1cm}
\begin{lstlisting}
Mesh meshd(Sphere(_center=Point(0.,0.,0.),_radius=1.,_nnodes=2,
                  _domain_name="Omega"),_triangle,1,_subdiv);
Domain omega=meshd.domain("Omega");
Space V(omega,P0,"V",false);
ClusterTree<FeDof> ct(V.feSpace()->dofs,_cardinalityBisection,10);
HMatrix<Real,FeDof> hm(ct,ct,nbox,nbox);
hm.saveStructureToFile("hmatrix.dat");
\end{lstlisting}
\vspace{.1cm}
\begin{figure}[H]
	\centering
	\includePict[width=15cm]{HMatrix.png}
	\caption{ HMatrix with a sphere cluster, non admissible blocks in red}
	\label{fig_hmatrix}
\end{figure}


\subsection{The {\classtitle HMatrixEntry} class}

The template {\class HMatrixEntry<I>} class is similar to the {\class MatrixEntry} class. In order to shadow the type of the HMatrix coefficients it handles pointers to {\class HMatrix<real\_t,I>},  {\class HMatrix<complex\_t,I>},  {\class HMatrix<Matrix<real\_t>,I>} and {\class HMatrix<Matrix<complex\_t>,I>}:
 \vspace{.1cm}
 \begin{lstlisting}
 class HMatrixEntry
 {
 public :
 ValueType valueType_;                     //entries value type
 StrucType strucType_;                     //entries structure type
 HMatrix<real_t,I>* rEntries_p;            //pointer to real HMatrix
 HMatrix<complex_t,I>* cEntries_p;         //pointer to complex HMatrix
 HMatrix<Matrix<real_t>,I>* rmEntries_p;   //pointer to HMatrix of real matrix
 HMatrix<Matrix<complex_t>,I>* cmEntries_p;//pointer to HMatrix of cmplx matrix
 dimPair nbOfComponents;             //nb of rows and columns of matrix values
 ...
 \end{lstlisting}
 \vspace{.1cm}
 The following constructor/destructor stuff is available:
 \vspace{.1cm}
 \begin{lstlisting}
 HMatrixEntry(ValueType, StrucType, ClusterTree<I>&, ClusterTree<I>&,
              number_t, number_t, number_t =1, number_t=1,
              SymType sy = _noSymmetry);
 HMatrixEntry(const HMatrixEntry<I>&);              //copy constructor
 HMatrixEntry<I>& operator=(const HMatrixEntry<I>&);//assign operator
 ~HMatrixEntry(){clear();}                          //destructor
 void copy(const HMatrixEntry<I>&);                 //full copy of members
 void clear();                                      //deallocate memory used
 \end{lstlisting}
 \vspace{.1cm}
 To preserve integrity of {\class HMatrixEntry} class, the copy constructor and assignment operator do a hard copy of {\class HMatrix} objects pointed as soon as the pointer is not null.\\

 The class provides some useful functionalities:
 \vspace{.1cm}
 \begin{lstlisting}
 void setClusterRow(ClusterTree<I>*);  //change the ClusterTree row pointer
 void setClusterCol(ClusterTree<I>*);  //change the ClusterTree col pointer
 template <typename T>
 HMatrix<T,I>& getHMatrix() const      //get HMatrix object
 real_t norm2() const;                 //Frobenius norm
 real_t norminfty() const;             //infinite norm
 void print(ostream&) const;           //print HMatrixEntry
 void printSummary(ostream&) const;    //print HMatrixEntry in brief
\end{lstlisting}
\vspace{.1cm}
 The member function {\tt getHMatrix} has {\tt real\_t, complex\_t, Matrix<real\_t>} and {\tt Matrix<complex\_t>} specializations.

 \displayInfos{library=hierarchicalMatrix, header=HMatrix.hpp, implementation=HMatrix.cpp, test={test\_HMatrix.cpp}, header dep={ApproximateMatrix.hpp, ClusterTree.hpp, largeMatrix.h ,config.h, utils.h}}

\subsection{The {\classtitle HMatrixIM} class}

In order to inform the code that HMatrix has to be used, a special class of integration method has been developed : {\class HMatrixIM }. It inherits from  the {\class DoubleIM} class inheriting itself  from the {\class IntegrationMethod} class. Attaching such integration method to a BEM bilinear form, HMatrix computation algorithm will be involved.
\vspace{.1cm}
\begin{lstlisting}
class HMatrixIM : public DoubleIM
{public:
 ClusterTree<FeDof>* rowCluster_;   //row cluster pointer
 ClusterTree<FeDof>* colCluster_;   //col cluster pointer
 mutable bool deletePointers_;      //flag enabling pointers deallocation
 ClusteringMethod clusterMethod;    //clustering method
 HMApproximationMethod hmAppMethod; //type of approximation of admissible blocks
 number_t minRowSize, minColSize;   //minimum row/col size of leaf matrix
 number_t maxRank;                  //maximal rank of approximate matrices
 real_t epsRank;                    //precision used in compression
 IntegrationMethod* intgMethod;     //real integration method
 ...
};
\end{lstlisting}
\vspace{.1cm}
This class provides some constructors
\vspace{.1cm}
\begin{lstlisting}
HMatrixIM(ClusteringMethod clm, number_t minRow, number_t minCol,
          HMApproximationMethod hmap, number_t maxr, IntegrationMethod& im);
HMatrixIM(ClusteringMethod clm, number_t minRow, number_t minCol,
          HMApproximationMethod hmap, real_t epsr, IntegrationMethod& im);
HMatrixIM(HMApproximationMethod hmap,  number_t maxr, IntegrationMethod& im,
          ClusterTree<FeDof>& rowC, ClusterTree<FeDof>& colC);
HMatrixIM(HMApproximationMethod hmap,  real_t epsr, IntegrationMethod& im,
          ClusterTree<FeDof>& rowC, ClusterTree<FeDof>& colC);
\end{lstlisting}
\vspace{.1cm}
The two last constructors allow to load some row/column clusters whereas the two first ones provide the parameters to build them. When building {\class HMatrixIM} object from clustering parameters, the {\tt deletePointers} flag is true, so destroying the {\class HMatrixIM} object induces the deallocation of the cluster pointers. When building {\class HMatrixIM} object from cluster pointers, these pointers will not be deallocated by destructor.  \\

The class has also a clear method to delete the cluster pointers and a destructor that clears the cluster pointers if {\tt deletePointers} is true:
\vspace{.1cm}
\begin{lstlisting}
~HMatrixIM();
 void clear();
\end{lstlisting}
\vspace{.1cm}
\goodbreak
The {\class HMatrixIM} class is used as following:
\vspace{.1cm}
\begin{lstlisting}
Mesh meshd(Sphere(_center=Point(0.,0.,0.),_radius=1.,_nnodes=9,
           _domain_name="Omega"), _triangle,1,_subdiv);
Domain omega=meshd.domain("Omega");
Space W(omega,P0,"V",false);
Unknown u(W,"u"); TestFunction v(u,"v");
Kernel G=Laplace3dKernel();
SauterSchwabIM ssim(5,5,4,3,2.,4.);
HMatrixIM him(_cardinalityBisection, 20, 20,_acaplus,0.00001,ssim);
BilinearForm alf=intg(omega,omega,u*Gl*v,him);
TermMatrix A(alf,"A");
\end{lstlisting}
\vspace{.1cm}

\begin{remark2}
Because the compression methods do not work yet for matrix of matrices, {\class HMatrix} is not fully operational when dealing with some problems with vector unknown.
\end{remark2}

\displayInfos{library=hierarchicalMatrix, header=HMatrix.hpp, implementation=HMatrix.cpp, test={test\_HMatrix.cpp}, header dep={ApproximateMatrix.hpp, ClusterTree.hpp, largeMatrix.h ,config.h, utils.h}}