GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/H2Matrix/H2Matrix.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 99 103 96.1%
Functions: 20 20 100.0%
Branches: 59 89 66.3%
Line Branch Exec Source
1 // This file is part of Bembel, the higher order C++ boundary element library.
2 //
3 // Copyright (C) 2022 see <http://www.bembel.eu>
4 //
5 // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
6 // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
7 // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
8 // source code is subject to the GNU General Public License version 3 and
9 // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
10 // information.
11 #ifndef BEMBEL_SRC_H2MATRIX_H2MATRIX_HPP_
12 #define BEMBEL_SRC_H2MATRIX_H2MATRIX_HPP_
13
14 /**
15 * \class H2Matrix
16 * \brief Hierarchical Matrix class, which extends the EigenBase class.
17 *
18 * The idea is to provide an easy to use interface to the H2-matrix
19 * from the fast boundary element method. At the moment, we inherit the
20 * traits of an Eigen::SparseMatrix, since this seems to be the minimum
21 * properties for a derived object to ensure that the matrix-vector
22 * multiplication can be specialised for H2Matrix.
23 * In particular, this allows for the use of the Eigen iterative solvers
24 * with a Hierarchical matrix.
25 **/
26 namespace Eigen {
27 /// forward definition of the H2Matrix Class in order to define traits
28 template <typename ScalarT>
29 class H2Matrix;
30 /// small modification of internal traits of SparseMatrix
31 namespace internal {
32 template <typename ScalarT>
33 struct traits<H2Matrix<ScalarT>> {
34 typedef ScalarT Scalar;
35 typedef int StorageIndex;
36 typedef H2 StorageKind;
37 typedef MatrixXpr XprKind;
38 enum {
39 RowsAtCompileTime = Dynamic,
40 ColsAtCompileTime = Dynamic,
41 MaxRowsAtCompileTime = Dynamic,
42 MaxColsAtCompileTime = Dynamic,
43 Flags = NestByRefBit
44 };
45 };
46
47 // this struct is necessary for compatibility with iterative solvers
48 template <typename ScalarT>
49 struct is_ref_compatible<H2Matrix<ScalarT>> {
50 enum { value = false };
51 };
52 template <typename ScalarT>
53 struct is_ref_compatible<const H2Matrix<ScalarT>>
54 : is_ref_compatible<H2Matrix<ScalarT>> {};
55
56 } // namespace internal
57
58 // actual definition of the class
59 template <typename ScalarT>
60 class H2Matrix : public H2MatrixBase<H2Matrix<ScalarT>> {
61 public:
62 //////////////////////////////////////////////////////////////////////////////
63 /// Eigen related things
64 //////////////////////////////////////////////////////////////////////////////
65 // Required typedefs, constants and so on.
66 typedef ScalarT Scalar;
67 typedef typename NumTraits<ScalarT>::Real RealScalar;
68 typedef Index StorageIndex;
69 enum {
70 ColsAtCompileTime = Dynamic,
71 MaxColsAtCompileTime = Dynamic,
72 IsRowMajor = false,
73 Flags = NestByRefBit
74 };
75 // Minimum specialisation of EigenBase methods
76 826 Index rows() const { return transformation_matrix_.cols(); }
77 500 Index cols() const { return transformation_matrix_.cols(); }
78 // Definition of the matrix multiplication
79 template <typename Rhs>
80 228 Product<H2Matrix, Rhs, AliasFreeProduct> operator*(
81 const MatrixBase<Rhs>& x) const {
82 228 return Product<H2Matrix, Rhs, AliasFreeProduct>(*this, x.derived());
83 }
84 //////////////////////////////////////////////////////////////////////////////
85 /// constructors
86 //////////////////////////////////////////////////////////////////////////////
87
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 H2Matrix() {}
88 /**
89 * \brief Assemble H2-Matrix for linear operator linOp and AnsatzSpace
90 * ansatz_space with number_of_points interpolation points in one direction of
91 * the unit square (standard is number_of_points=9)
92 */
93 template <typename Derived>
94 5 void init_H2Matrix(const Derived& linOp,
95 const Bembel::AnsatzSpace<Derived>& ansatz_space,
96 int number_of_points = 9) {
97 // get transformation matrix from ansatz space
98
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 transformation_matrix_ = ansatz_space.get_transformation_matrix();
99 /**
100 * \todo Juergen Discuss with Michael where to initialize the parameters:
101 * min_cluster_level depends on number_of_points, but this is not
102 * implemented yet, also, what do the parameter mean?
103 */
104 // initialize block-cluster-trees and get parameters
105 5 const int vector_dimension = Bembel::getFunctionSpaceVectorDimension<
106 Bembel::LinearOperatorTraits<Derived>::Form>();
107
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 block_cluster_tree_.resize(vector_dimension, vector_dimension);
108 {
109
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Bembel::BlockClusterTree<Scalar> bt(linOp, ansatz_space);
110
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 5 times.
 11 for (int i = 0; i < vector_dimension; ++i)
111
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 6 times.
 14 for (int j = 0; j < vector_dimension; ++j)
112 8 block_cluster_tree_(i, j) =
113 bt; // .init_BlockClusterTree(linOp, ansatz_space);
114 5 }
115
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 auto parameters = block_cluster_tree_(0, 0).get_parameters();
116 // compute transfer and moment matrices for fmm
117 5 int cluster_level = parameters.max_level_ - parameters.min_cluster_level_;
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (cluster_level < 0) cluster_level = 0;
119 5 int cluster_refinement = parameters.min_cluster_level_;
120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 5 if (cluster_refinement > parameters.max_level_)
121 ✗ cluster_refinement = parameters.max_level_;
122
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 fmm_transfer_matrices_ = Bembel::H2Multipole::computeTransferMatrices<
123 Bembel::H2Multipole::ChebychevRoots>(number_of_points);
124
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 fmm_moment_matrix_ = Bembel::H2Multipole::
125 Moment2D<Bembel::H2Multipole::ChebychevRoots, Derived>::compute2DMoment(
126 ansatz_space.get_superspace(), cluster_level, cluster_refinement,
127 number_of_points);
128 // compute interpolation points
129 5 Eigen::VectorXd interpolation_points1D =
130
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Bembel::H2Multipole::ChebychevRoots(number_of_points).points_;
131
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Eigen::MatrixXd interpolation_points2D =
132 Bembel::H2Multipole::interpolationPoints2D(interpolation_points1D);
133 // compute content of tree leafs
134 5 int polynomial_degree = ansatz_space.get_polynomial_degree();
135 5 int polynomial_degree_plus_one_squared =
136 5 (polynomial_degree + 1) * (polynomial_degree + 1);
137
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Bembel::GaussSquare<Bembel::Constants::maximum_quadrature_degree> GS;
138 5 auto super_space = ansatz_space.get_superspace();
139 5 auto ffield_deg = linOp.get_FarfieldQuadratureDegree(polynomial_degree);
140
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 std::vector<ElementSurfacePoints> ffield_qnodes =
141 Bembel::DuffyTrick::computeFfieldQnodes(super_space, GS[ffield_deg]);
142 5 const int NumberOfFMMComponents =
143 Bembel::LinearOperatorTraits<Derived>::NumberOfFMMComponents;
144 #pragma omp parallel
145 {
146 #pragma omp single
147 {
148 // build vector of iterators
149 std::vector<
150 typename std::vector<Bembel::BlockClusterTree<Scalar>*>::iterator>
151 5 leafs;
152
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 leafs.resize(vector_dimension * vector_dimension);
153
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 5 times.
 11 for (int i = 0; i < vector_dimension; ++i)
154
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 6 times.
 14 for (int j = 0; j < vector_dimension; ++j)
155 8 leafs[i * vector_dimension + j] =
156 16 block_cluster_tree_(j, i).lbegin();
157 // iterate over leafs
158
2/2
✓ Branch 4 taken 2880 times.
✓ Branch 5 taken 5 times.
 2885 for (; leafs[0] != block_cluster_tree_(0, 0).lend();) {
159 #pragma omp task firstprivate(leafs)
160 {
161
2/3
✓ Branch 3 taken 1824 times.
✓ Branch 4 taken 1056 times.
✗ Branch 5 not taken.
 2880 switch ((*(leafs[0]))->get_cc()) {
162 // assemble dense matrix blocks
163 1824 case Bembel::BlockClusterAdmissibility::Dense: {
164 const Bembel::ElementTreeNode* cluster1 =
165 1824 (*(leafs[0]))->get_cluster1();
166 const Bembel::ElementTreeNode* cluster2 =
167 1824 (*(leafs[0]))->get_cluster2();
168 1824 int block_size =
169
1/2
✓ Branch 3 taken 1824 times.
✗ Branch 4 not taken.
 1824 std::distance(cluster2->begin(), cluster2->end()) *
170 polynomial_degree_plus_one_squared;
171 std::vector<
172 Eigen::Matrix<ScalarT, Eigen::Dynamic, Eigen::Dynamic>>
173 1824 F;
174
2/2
✓ Branch 0 taken 2976 times.
✓ Branch 1 taken 1824 times.
 4800 for (int i = 0; i < vector_dimension * vector_dimension; ++i)
175
1/2
✓ Branch 1 taken 2976 times.
✗ Branch 2 not taken.
 2976 F.push_back(
176
1/2
✓ Branch 1 taken 2976 times.
✗ Branch 2 not taken.
 5952 Eigen::Matrix<ScalarT, Eigen::Dynamic, Eigen::Dynamic>(
177 block_size, block_size));
178 // iterate over elements in dense matrix block
179 1824 unsigned int cl1index = 0;
180 1824 unsigned int cl2index = 0;
181
2/2
✓ Branch 4 taken 7296 times.
✓ Branch 5 taken 1824 times.
 9120 for (const auto& element1 : *cluster1) {
182 7296 cl2index = 0;
183
2/2
✓ Branch 5 taken 29184 times.
✓ Branch 6 taken 7296 times.
 36480 for (const auto& element2 : *cluster2) {
184 Eigen::Matrix<ScalarT, Eigen::Dynamic, Eigen::Dynamic>
185 ✗ intval(vector_dimension *
186 polynomial_degree_plus_one_squared,
187
1/2
✓ Branch 1 taken 29184 times.
✗ Branch 2 not taken.
 29184 vector_dimension *
188 polynomial_degree_plus_one_squared);
189 // do integration
190
1/2
✓ Branch 1 taken 29184 times.
✗ Branch 2 not taken.
 29184 Bembel::DuffyTrick::evaluateBilinearForm(
191 linOp, super_space, element1, element2, GS,
192 29184 ffield_qnodes[element1.id_],
193 29184 ffield_qnodes[element2.id_], &intval);
194 // insert into dense matrices of all block cluster trees
195
2/2
✓ Branch 0 taken 35328 times.
✓ Branch 1 taken 29184 times.
 64512 for (int i = 0; i < vector_dimension; ++i)
196
2/2
✓ Branch 0 taken 47616 times.
✓ Branch 1 taken 35328 times.
 82944 for (int j = 0; j < vector_dimension; ++j)
197 47616 F[i * vector_dimension + j].block(
198 47616 polynomial_degree_plus_one_squared * cl1index,
199
1/2
✓ Branch 1 taken 47616 times.
✗ Branch 2 not taken.
 47616 polynomial_degree_plus_one_squared * cl2index,
200 polynomial_degree_plus_one_squared,
201
1/2
✓ Branch 1 taken 47616 times.
✗ Branch 2 not taken.
 47616 polynomial_degree_plus_one_squared) =
202 47616 intval.block(j * polynomial_degree_plus_one_squared,
203
1/2
✓ Branch 1 taken 47616 times.
✗ Branch 2 not taken.
 47616 i * polynomial_degree_plus_one_squared,
204 polynomial_degree_plus_one_squared,
205 polynomial_degree_plus_one_squared);
206 29184 ++cl2index;
207 }
208 7296 ++cl1index;
209 }
210
2/2
✓ Branch 0 taken 2976 times.
✓ Branch 1 taken 1824 times.
 4800 for (int i = 0; i < vector_dimension * vector_dimension; ++i)
211
1/2
✓ Branch 5 taken 2976 times.
✗ Branch 6 not taken.
 2976 (*(leafs[i]))->get_leaf().set_F(F[i]);
212 1824 } break;
213 // interpolation for low-rank blocks
214 1056 case Bembel::BlockClusterAdmissibility::LowRank: {
215
1/2
✓ Branch 1 taken 1056 times.
✗ Branch 2 not taken.
 1056 auto F = Bembel::H2Multipole::interpolateKernel<Derived>(
216 linOp, super_space, interpolation_points2D,
217 1056 *((*(leafs[0]))->get_cluster1()),
218 1056 *((*(leafs[0]))->get_cluster2()));
219
2/2
✓ Branch 0 taken 1248 times.
✓ Branch 1 taken 1056 times.
 2304 for (int i = 0; i < vector_dimension; ++i) {
220
2/2
✓ Branch 0 taken 1632 times.
✓ Branch 1 taken 1248 times.
 2880 for (int j = 0; j < vector_dimension; ++j) {
221 1632 (*(leafs[i * vector_dimension + j]))
222 1632 ->get_leaf()
223
4/8
✓ Branch 1 taken 768 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1632 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1632 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 864 times.
✗ Branch 11 not taken.
 2496 .set_F(F.block(j * interpolation_points2D.rows() *
224 NumberOfFMMComponents,
225 1632 i * interpolation_points2D.rows() *
226 NumberOfFMMComponents,
227 768 interpolation_points2D.rows() *
228 NumberOfFMMComponents,
229 768 interpolation_points2D.rows() *
230 NumberOfFMMComponents));
231 1632 (*(leafs[i * vector_dimension + j]))
232 1632 ->get_leaf()
233 1632 .set_low_rank_flag(true);
234 }
235 }
236 1056 } break;
237 // this leaf is not a low-rank block and not a dense block,
238 // thus it is not a leaf -> error
239 ✗ default:
240 ✗ assert(0 && "This should never happen");
241 break;
242 }
243 }
244 // increment leafs of all block-cluster trees
245
2/2
✓ Branch 4 taken 4608 times.
✓ Branch 5 taken 2880 times.
 7488 for (auto leafsit = leafs.begin(); leafsit != leafs.end(); ++leafsit)
246 4608 ++(*leafsit);
247 }
248 5 }
249 }
250 5 }
251 //////////////////////////////////////////////////////////////////////////////
252 /// methods
253 //////////////////////////////////////////////////////////////////////////////
254 Eigen::Matrix<ScalarT, Eigen::Dynamic, Eigen::Dynamic> get_dense() const {
255 Eigen::Matrix<ScalarT, Eigen::Dynamic, Eigen::Dynamic> dense(rows(),
256 cols());
257 for (int i = 0; i < cols(); ++i) {
258 Eigen::VectorXd unit(cols());
259 unit.setZero();
260 unit(i) = 1.;
261 dense.col(i) = (*this) * unit;
262 }
263 return dense;
264 }
265 //////////////////////////////////////////////////////////////////////////////
266 /// getter
267 //////////////////////////////////////////////////////////////////////////////
268 606 const Eigen::SparseMatrix<double> get_transformation_matrix() const {
269 606 return transformation_matrix_;
270 }
271 303 const Eigen::MatrixXd get_fmm_transfer_matrices() const {
272 303 return fmm_transfer_matrices_;
273 }
274 303 const std::vector<Eigen::MatrixXd> get_fmm_moment_matrix() const {
275 303 return fmm_moment_matrix_;
276 }
277 const Bembel::GenericMatrix<Bembel::BlockClusterTree<ScalarT>>&
278 748074 get_block_cluster_tree() const {
279 748074 return block_cluster_tree_;
280 }
281 //////////////////////////////////////////////////////////////////////////////
282 /// private members
283 //////////////////////////////////////////////////////////////////////////////
284 private:
285 // we declare functionality which has not been implemented (yet)
286 // to be private
287 H2Matrix(const H2Matrix<ScalarT>& H);
288 H2Matrix(H2Matrix<ScalarT>&& H);
289 H2Matrix& operator=(const H2Matrix<ScalarT>& H);
290 H2Matrix& operator=(H2Matrix<ScalarT>&& H);
291
292 Eigen::SparseMatrix<double> transformation_matrix_;
293 Bembel::GenericMatrix<Bembel::BlockClusterTree<ScalarT>> block_cluster_tree_;
294 Eigen::MatrixXd fmm_transfer_matrices_;
295 std::vector<Eigen::MatrixXd> fmm_moment_matrix_;
296 };
297
298 namespace internal {
299
300 // adaption from SparseMatrix from Eigen
301 template <typename Scalar>
302 struct evaluator<H2Matrix<Scalar>> : evaluator<H2MatrixBase<H2Matrix<Scalar>>> {
303 typedef evaluator<H2MatrixBase<H2Matrix<Scalar>>> Base;
304 typedef H2Matrix<Scalar> H2MatrixType;
305 evaluator() : Base() {}
306 explicit evaluator(const H2MatrixType& mat) : Base(mat) {}
307 };
308
309 } // namespace internal
310 } // namespace Eigen
311
312 #endif // BEMBEL_SRC_H2MATRIX_H2MATRIX_HPP_
313