1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{HMatrix} 18 19Hierarchical 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).\\ 20 21To manage this tree representation, three template classes are provided: 22\begin{itemize} 23 \item the {\class HMatrixNode<T,I>} class dealing with a node of the matrix tree 24 \item the {\class HMatrix<T,I>} front end class 25 \item the {\class HMatrixEntry<I>} interface class that handles some pointers to {\class HMatrix<real\_t,I>}, {\class HMatrix<complex\_t,I>}, ... . 26\end{itemize} 27These 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. 28 29\subsection{The {\classtitle HMatrixNode} class} 30 31\vspace{.1cm} 32\begin{lstlisting} 33template <typename T, typename I> 34class HMatrixNode 35{public : 36 HMatrixNode<T,I>* parent_; //pointer to its parent, if 0 root node 37 HMatrixNode<T,I>* child_; //pointer to its first child, if 0 no child 38 HMatrixNode<T,I>* next_; //pointer to its brother, if 0 no brother 39 ClusterNode<I>* rowNode_; //row cluster node to access to row numbering 40 ClusterNode<I>* colNode_; //column cluster node to access to col numbering 41 LargeMatrix<T>* mat_; //pointer to a LargeMatrix 42 ApproximateMatrix<T>* appmat_; //pointer to an ApproximateMatrix 43 bool admissible_; //block can be approximated 44 number_t depth_; //depth of node in tree, root has 0 depth 45 number_t rowSub_, colSub_; //sub block numbering (i,j) 46 bool isDiag_; //true if matrix node is located on the diagonal 47 \end{lstlisting} 48 \vspace{.1cm} 49 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\_}. \\ 50 51 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: 52\vspace{.1cm} 53\begin{lstlisting} 54enum HMatrixMethod {_standardHM, _denseHM}; 55enum HMAdmissibilityRule {_boxesRule}; 56\end{lstlisting} 57\vspace{.1cm} 58The {\class HMatrixNode} class has few constructors: 59\vspace{.1cm} 60\begin{lstlisting} 61HMatrixNode(); //default constructor 62HMatrixNode(HMatrixNode<T,I>*, HMatrixNode<T,I>*, HMatrixNode<T,I>*, number_t, 63 ClusterNode<I>*, ClusterNode<I>*, number_t, number_t, 64 LargeMatrix<T>* , ApproximateMatrix<T>* );//full constructor 65HMatrixNode(HMatrixNode<T,I>*, number_t); //node constructor 66HMatrixNode(LargeMatrix<T>*, HMatrixNode<T,I>*, 67 number_t); //LargeMatrix leaf constructor 68HMatrixNode(ApproximateMatrix<T>*, HMatrixNode<T,I>*, 69 number_t); //ApproximateMatrix leaf constructor 70HMatrixNode(const HMatrixNode<T,I>&); //copy constructor 71HMatrixNode<T,I>& operator=(const HMatrixNode<T,I>&); //assign operator 72void copy(const HMatrixNode<T,I>&); //copy function 73\end{lstlisting} 74\vspace{.1cm} 75 and some stuff related to destructor or cleaner: 76 \vspace{.1cm} 77\begin{lstlisting} 78~HMatrixNode(); //destructor 79void clear(); //clear all except ClusterNode pointers 80void clearMatrices(); //clear (deallocate) only matrices 81\end{lstlisting} 82\vspace{.1cm} 83\begin{remark} 84 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. 85\end{remark}\\ 86 87The most important function is 88\vspace{.1cm} 89\begin{lstlisting} 90void divide(number_t rmin, number_t cmin, number_t maxdepth = 0, 91 HMAdmissibilityRule=_boxesRule, bool sym=false); 92\end{lstlisting} 93\vspace{.1cm} 94that creates recursively the sub-tree of the current node according to the admissibility rule given by the enumeration 95\vspace{.1cm} 96\begin{lstlisting} 97enum HMAdmissibilityRule {_boxesRule}; 98\end{lstlisting} 99\vspace{.1cm} 100Up 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: 101$$\text{diam}(B_r) \leq 2*\eta*\text{dist}(B_r,B_c).$$ 102\begin{remark} 103 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. 104\end{remark}\\ 105 106Once the tree is built, some various member functions are available: 107\vspace{.1cm} 108\begin{lstlisting} 109number_t numberOfRows() const; //number of rows counted in T 110number_t numberOfCols() const; //number of columns counted in T 111dimPair dimValues() const; //dimensions of matrix values 112void setClusterRow(ClusterNode<I>*); //update row cluster node pointers 113void setClusterCol(ClusterNode<I>*); //update col cluster node pointers 114number_t nbNonZero() const; //number of T coefficients used 115void getLeaves(list<HMatrixNode<T,I>* >&, 116 bool = true) const; //get leaves as a list 117\end{lstlisting} 118\vspace{.1cm} 119The {\class HMatrixNode} class implements all the stuff required by algebraic operations on {\class HMatrix} class: 120\vspace{.1cm} 121\begin{lstlisting} 122vector<T>& multMatrixVectorNode(const vector<T>&, vector<T>&, 123 SymType symt=_noSymmetry) const; 124vector<T>& multMatrixVector(const vector<T>&, vector<T>&, 125 SymType symt=_noSymmetry) const; 126real_t norm2() const; //Frobenius norm 127real_t norminfty() const; //infinite norm 128\end{lstlisting} 129\vspace{.1cm} 130and by printing operations 131\vspace{.1cm} 132\begin{lstlisting} 133void printNode(ostream&) const; //print current node contents 134void print(ostream&) const; //print tree from current node 135void printStructure(ostream&, number_t, number_t, 136 bool all=false, bool shift=false) const; 137\end{lstlisting} 138\vspace{.1cm} 139 140Finally, for performance measurement purpose, the class manages also some static data and static functions: 141\vspace{.1cm} 142\begin{lstlisting} 143static bool counterOn; //flag to activate the operations counter 144static number_t counter; //operations counter 145static void initCounter(number_t n=0); //init the counter and enable it 146static void stopCounter(); //disable counter 147\end{lstlisting} 148\vspace{.1cm} 149 150\subsection{The {\classtitle HMatrix} class} 151 152The {\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: 153 \vspace{.1cm} 154 \begin{lstlisting} 155template <typename T, typename I> 156class HMatrix 157{private : 158 HMatrixNode<T,I>* root_; //root node of the tree representation 159 ClusterTree<I>* rowCT_; //row cluster tree 160 ClusterTree<I>* colCT_; //col cluster tree 161 public : 162 string_t name; //optional name, useful for doc 163 ValueType valueType_; //type of values (real, complex) 164 StrucType strucType_; //structure of values (scalar, vector, matrix) 165 HMatrixMethod method_; //method to define admissible block 166 HMAdmissibilityRule admRule_; //block admissible rule 167 real_t eta_; //ratio in the admissibility criteria 168 number_t rowmin_, colmin_;//minimum size of block matrix 169 SymType sym_; //type of symmetry 170 number_t depth; //maximal depth (info) 171 number_t nbNodes; //number of nodes (info) 172 number_t nbLeaves; //number of leaves (info) 173 number_t nbAdmissibles; //number of admissible blocks (info) 174 number_t nbAppMatrices; //number of approximate matrices (info) 175 ... 176 \end{lstlisting} 177 \vspace{.1cm} 178 The data member {\tt depth, nbNodes, nbLeaves, nbAdmissibles, nbAppMatrices} are set by the division algorithm or computed on demand.\\ 179 180 The class redefines the default constructor (initializing data members) and an explicit constructor filling the cluster parameters and building the tree structure recursively: 181 \vspace{.1cm} 182 \begin{lstlisting} 183HMatrix(); 184HMatrix(ClusterTree<I>&, ClusterTree<I>&, number_t, number_t, number_t =0, 185 const string_t& = "", SymType =_noSymmetry, 186 HMatrixMethod =_standardHM, HMAdmissibilityRule=_boxesRule, real_t =1.); 187void buildTree(); //build the tree 188void updateInfo(); //update tree info (depth, ...) 189HMatrix(const HMatrix<T,I>&); //copy constructor 190HMatrix<T,I>& operator=(const HMatrix<T,I>&); //assign operator 191~HMatrix(); //destructor 192void copy(const HMatrix<T,I>&); //copy all 193void clear(); //clear all 194void clearMatrices(); //deallocate matrix pointers 195\end{lstlisting} 196\vspace{.1cm} 197 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!\\ 198 199The {\class HMatrix} class provides some stuff to get or extract various informations about the tree structure: 200\vspace{.1cm} 201\begin{lstlisting} 202number_t numberOfRows() const; //number of rows counted in T 203number_t numberOfCols() const; //number of cols counted in T 204dimPair dimValues() const; //dimensions of values, (1,1) when scalar 205number_t nbNonZero() const; //number of coefficients used 206const ClusterTree<I>* rowTree() const; //access to row ClusterTree pointer 207const ClusterTree<I>* colTree() const; //access to col ClusterTree pointer 208void setClusterRow(ClusterTree<I>*); //change the row ClusterTree pointer 209void setClusterCol(ClusterTree<I>*); //change the col ClusterTree pointer 210pair<number_t, number_t> averageSize() const; //average size of leaves 211number_t averageRank() const; //average rank of admissible leaves 212list<HMatrixNode<T,I>* > getLeaves(bool = true) const; //get all leaves 213\end{lstlisting} 214\vspace{.1cm} 215It 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). 216\vspace{.1cm} 217\begin{lstlisting} 218void load(const LargeMatrix<T>&,HMApproximationMethod = _noHMApproximation); 219LargeMatrix<T> toLargeMatrix(StorageType st=_dense, AccessType at=_row) const; 220\end{lstlisting} 221\vspace{.1cm} 222Some algebraic operations on {\class HMatrix} are available: 223\vspace{.1cm} 224\begin{lstlisting} 225vector<T>& multMatrixVector(const vector<T>&, vector<T>&) const; 226vector<T>& multMatrixVectorOmp(const vector<T>&, vector<T>&) const; 227void addFELargeMatrix(const LargeMatrix<T>&); //add FE LargeMatrix to current 228real_t norm2() const; //Frobenius norm 229real_t norminfty() const; //infinite norm 230\end{lstlisting} 231\vspace{.1cm} 232Note 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.\\ 233 234The {\class Hmatrix} class provides the following print stuff: 235\vspace{.1cm} 236\begin{lstlisting} 237void print(std::ostream&) const; 238void printSummary(std::ostream&) const; 239void printStructure(std::ostream&,bool all=false, bool =false) const; 240void saveStructureToFile(const string_t&) const; 241ostream& operator<<(ostream& os, const HMatrix<X,J>& hm); 242\end{lstlisting} 243\vspace{.2cm} 244As 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 245\vspace{.1cm} 246\begin{lstlisting} 247Mesh meshd(Sphere(_center=Point(0.,0.,0.),_radius=1.,_nnodes=2, 248 _domain_name="Omega"),_triangle,1,_subdiv); 249Domain omega=meshd.domain("Omega"); 250Space V(omega,P0,"V",false); 251ClusterTree<FeDof> ct(V.feSpace()->dofs,_cardinalityBisection,10); 252HMatrix<Real,FeDof> hm(ct,ct,nbox,nbox); 253hm.saveStructureToFile("hmatrix.dat"); 254\end{lstlisting} 255\vspace{.1cm} 256\begin{figure}[H] 257 \centering 258 \includePict[width=15cm]{HMatrix.png} 259 \caption{ HMatrix with a sphere cluster, non admissible blocks in red} 260 \label{fig_hmatrix} 261\end{figure} 262 263 264\subsection{The {\classtitle HMatrixEntry} class} 265 266The 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>}: 267 \vspace{.1cm} 268 \begin{lstlisting} 269 class HMatrixEntry 270 { 271 public : 272 ValueType valueType_; //entries value type 273 StrucType strucType_; //entries structure type 274 HMatrix<real_t,I>* rEntries_p; //pointer to real HMatrix 275 HMatrix<complex_t,I>* cEntries_p; //pointer to complex HMatrix 276 HMatrix<Matrix<real_t>,I>* rmEntries_p; //pointer to HMatrix of real matrix 277 HMatrix<Matrix<complex_t>,I>* cmEntries_p;//pointer to HMatrix of cmplx matrix 278 dimPair nbOfComponents; //nb of rows and columns of matrix values 279 ... 280 \end{lstlisting} 281 \vspace{.1cm} 282 The following constructor/destructor stuff is available: 283 \vspace{.1cm} 284 \begin{lstlisting} 285 HMatrixEntry(ValueType, StrucType, ClusterTree<I>&, ClusterTree<I>&, 286 number_t, number_t, number_t =1, number_t=1, 287 SymType sy = _noSymmetry); 288 HMatrixEntry(const HMatrixEntry<I>&); //copy constructor 289 HMatrixEntry<I>& operator=(const HMatrixEntry<I>&);//assign operator 290 ~HMatrixEntry(){clear();} //destructor 291 void copy(const HMatrixEntry<I>&); //full copy of members 292 void clear(); //deallocate memory used 293 \end{lstlisting} 294 \vspace{.1cm} 295 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.\\ 296 297 The class provides some useful functionalities: 298 \vspace{.1cm} 299 \begin{lstlisting} 300 void setClusterRow(ClusterTree<I>*); //change the ClusterTree row pointer 301 void setClusterCol(ClusterTree<I>*); //change the ClusterTree col pointer 302 template <typename T> 303 HMatrix<T,I>& getHMatrix() const //get HMatrix object 304 real_t norm2() const; //Frobenius norm 305 real_t norminfty() const; //infinite norm 306 void print(ostream&) const; //print HMatrixEntry 307 void printSummary(ostream&) const; //print HMatrixEntry in brief 308\end{lstlisting} 309\vspace{.1cm} 310 The member function {\tt getHMatrix} has {\tt real\_t, complex\_t, Matrix<real\_t>} and {\tt Matrix<complex\_t>} specializations. 311 312 \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}} 313 314\subsection{The {\classtitle HMatrixIM} class} 315 316In 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. 317\vspace{.1cm} 318\begin{lstlisting} 319class HMatrixIM : public DoubleIM 320{public: 321 ClusterTree<FeDof>* rowCluster_; //row cluster pointer 322 ClusterTree<FeDof>* colCluster_; //col cluster pointer 323 mutable bool deletePointers_; //flag enabling pointers deallocation 324 ClusteringMethod clusterMethod; //clustering method 325 HMApproximationMethod hmAppMethod; //type of approximation of admissible blocks 326 number_t minRowSize, minColSize; //minimum row/col size of leaf matrix 327 number_t maxRank; //maximal rank of approximate matrices 328 real_t epsRank; //precision used in compression 329 IntegrationMethod* intgMethod; //real integration method 330 ... 331}; 332\end{lstlisting} 333\vspace{.1cm} 334This class provides some constructors 335\vspace{.1cm} 336\begin{lstlisting} 337HMatrixIM(ClusteringMethod clm, number_t minRow, number_t minCol, 338 HMApproximationMethod hmap, number_t maxr, IntegrationMethod& im); 339HMatrixIM(ClusteringMethod clm, number_t minRow, number_t minCol, 340 HMApproximationMethod hmap, real_t epsr, IntegrationMethod& im); 341HMatrixIM(HMApproximationMethod hmap, number_t maxr, IntegrationMethod& im, 342 ClusterTree<FeDof>& rowC, ClusterTree<FeDof>& colC); 343HMatrixIM(HMApproximationMethod hmap, real_t epsr, IntegrationMethod& im, 344 ClusterTree<FeDof>& rowC, ClusterTree<FeDof>& colC); 345\end{lstlisting} 346\vspace{.1cm} 347The 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. \\ 348 349The class has also a clear method to delete the cluster pointers and a destructor that clears the cluster pointers if {\tt deletePointers} is true: 350\vspace{.1cm} 351\begin{lstlisting} 352~HMatrixIM(); 353 void clear(); 354\end{lstlisting} 355\vspace{.1cm} 356\goodbreak 357The {\class HMatrixIM} class is used as following: 358\vspace{.1cm} 359\begin{lstlisting} 360Mesh meshd(Sphere(_center=Point(0.,0.,0.),_radius=1.,_nnodes=9, 361 _domain_name="Omega"), _triangle,1,_subdiv); 362Domain omega=meshd.domain("Omega"); 363Space W(omega,P0,"V",false); 364Unknown u(W,"u"); TestFunction v(u,"v"); 365Kernel G=Laplace3dKernel(); 366SauterSchwabIM ssim(5,5,4,3,2.,4.); 367HMatrixIM him(_cardinalityBisection, 20, 20,_acaplus,0.00001,ssim); 368BilinearForm alf=intg(omega,omega,u*Gl*v,him); 369TermMatrix A(alf,"A"); 370\end{lstlisting} 371\vspace{.1cm} 372 373\begin{remark2} 374Because 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. 375\end{remark2} 376 377\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}} 378