DUNE PDELab (git)

checkmatrixinterface.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_COMMON_CHECK_MATRIX_INTERFACE_HH
6#define DUNE_COMMON_CHECK_MATRIX_INTERFACE_HH
7
8#include <algorithm>
9#include <limits>
10
11#include <dune/common/dynvector.hh>
12#include <dune/common/exceptions.hh>
13#include <dune/common/ftraits.hh>
14#include <dune/common/fvector.hh>
15
16
17/*
18 * @file
19 * @brief This file provides an interface check for dense matrices.
20 * @author Christoph Gersbacher
21 */
22
23
24namespace Dune
25{
26
27 // External forward declarations for namespace Dune
28 // ------------------------------------------------
29
30 template< class, int, int > class FieldMatrix;
31 template< class, int > class DiagonalMatrix;
32
33} // namespace Dune
34
35
36
37namespace CheckMatrixInterface
38{
39
40 namespace Capabilities
41 {
42
43 // hasStaticSizes
44 // --------------
45
46 template< class Matrix >
47 struct hasStaticSizes
48 {
49 static const bool v = false;
50 static const int rows = ~0;
51 static const int cols = ~0;
52 };
53
54 template< class Matrix >
55 struct hasStaticSizes< const Matrix >
56 {
57 static const bool v = hasStaticSizes< Matrix >::v;
58 static const int rows = hasStaticSizes< Matrix >::rows;
59 static const int cols = hasStaticSizes< Matrix >::cols;
60 };
61
62
63
64 // isSquare
65 // --------
66
67 template< class Matrix >
68 struct isSquare
69 {
70 static const bool v = false;
71 };
72
73 template< class Matrix >
74 struct isSquare< const Matrix >
75 {
76 static const bool v = isSquare< Matrix >::v;
77 };
78
79
80
81 // Template specializations for Dune::FieldMatrix
82 // ----------------------------------------------
83
84 template< class K, int r, int c >
85 struct hasStaticSizes< Dune::FieldMatrix< K, r, c > >
86 {
87 static const bool v = true;
88 static const int rows = r;
89 static const int cols = c;
90 };
91
92 template< class K, int rows, int cols >
93 struct isSquare< Dune::FieldMatrix< K, rows, cols > >
94 {
95 static const bool v = ( rows == cols );
96 };
97
98
99
100 // Template specializations for Dune::DiagonalMatrix
101 // -------------------------------------------------
102
103 template< class K, int n >
104 struct hasStaticSizes< Dune::DiagonalMatrix<K,n> >
105 {
106 static const bool v = true;
107 static const int rows = n;
108 static const int cols = n;
109 };
110
111 template< class K, int n >
112 struct isSquare< Dune::DiagonalMatrix<K,n> >
113 {
114 static const bool v = true;
115 };
116
117 } // namespace Capabilities
118
119
120
121 // UseDynamicVector
122 // ----------------
123
124 template< class Matrix >
125 struct UseDynamicVector
126 {
127 typedef typename Matrix::value_type value_type;
128
129 typedef Dune::DynamicVector< value_type > domain_type;
130 typedef domain_type range_type;
131
132 static domain_type domain ( const Matrix &matrix, value_type v = value_type() )
133 {
134 return domain_type( matrix.M(), v );
135 }
136
137 static range_type range ( const Matrix &matrix, value_type v = value_type() )
138 {
139 return range_type( matrix.N(), v );
140 }
141 };
142
143
144
145 // UseFieldVector
146 // --------------
147
148 template< class K, int rows, int cols >
149 struct UseFieldVector
150 {
151 typedef K value_type;
152
153 typedef Dune::FieldVector< K, cols > domain_type;
154 typedef Dune::FieldVector< K, rows > range_type;
155
156 template< class Matrix >
157 static domain_type domain ( const Matrix &, value_type v = value_type() )
158 {
159 return domain_type( v );
160 }
161
162 template< class Matrix >
163 static range_type range ( const Matrix &, value_type v = value_type() )
164 {
165 return range_type( v );
166 }
167 };
168
169
170
171 // MatrixSizeHelper
172 // ----------------
173
174 template< class Matrix, bool hasStaticSizes = Capabilities::hasStaticSizes< Matrix >::v >
175 struct MatrixSizeHelper;
176
177 template< class Matrix >
178 struct MatrixSizeHelper< Matrix, false >
179 {
180 typedef typename Matrix::size_type size_type;
181 static const size_type rows ( const Matrix &matrix ) { return matrix.rows(); }
182 static const size_type cols ( const Matrix &matrix ) { return matrix.cols(); }
183 };
184
185 template< class Matrix >
186 struct MatrixSizeHelper< Matrix, true >
187 {
188 typedef typename Matrix::size_type size_type;
189 static const size_type rows ( const Matrix & ) { return Matrix::rows; }
190 static const size_type cols ( const Matrix & ) { return Matrix::cols; }
191 };
192
193
194
195 // CheckIfSquareMatrix
196 // -------------------
197
198 template< class Matrix, class Traits, bool isSquare = Capabilities::isSquare< Matrix >::v >
199 struct CheckIfSquareMatrix;
200
201 template< class Matrix, class Traits >
202 struct CheckIfSquareMatrix< Matrix, Traits, false >
203 {
204 static void apply ( const Matrix &) {}
205
206 static void apply ( Matrix &) {}
207 };
208
209 template< class Matrix, class Traits >
210 struct CheckIfSquareMatrix< Matrix, Traits, true >
211 {
212 typedef typename Matrix::value_type value_type;
213
214 static value_type tolerance ()
215 {
216 return value_type( 16 ) * std::numeric_limits< value_type >::epsilon();
217 }
218
219 static void apply ( const Matrix &matrix )
220 {
221 const value_type determinant = matrix.determinant();
222 if( determinant > tolerance() )
223 {
224 typename Traits::domain_type x = Traits::domain( matrix );
225 const typename Traits::range_type b = Traits::range( matrix );
226 matrix.solve( x, b );
227 }
228 }
229
230 static void apply ( Matrix &matrix )
231 {
232 apply( const_cast< const Matrix & >( matrix ) );
233 if( matrix.determinant() > tolerance() )
234 matrix.invert();
235 }
236 };
237
238
239
240 // CheckConstMatrix
241 // ----------------
242
243 template< class Matrix, class Traits >
244 struct CheckConstMatrix
245 {
246 // required type definitions
247 typedef typename Matrix::size_type size_type;
248
249 typedef typename Matrix::value_type value_type;
250 typedef typename Matrix::field_type field_type;
251 typedef typename Matrix::block_type block_type;
252
253 typedef typename Matrix::const_row_reference const_row_reference;
254
255 typedef typename Matrix::ConstIterator ConstIterator;
256
257 static void apply ( const Matrix &matrix )
258 {
259 checkDataAccess ( matrix );
260 checkIterators ( matrix );
261 checkLinearAlgebra ( matrix );
262 checkNorms ( matrix );
263 checkSizes ( matrix );
264 CheckIfSquareMatrix< Matrix, Traits >::apply( matrix );
265
266 // TODO: check comparison
267 // bool operator == ( const Matrix &other );
268 // bool operator != ( const Matrix &other );
269 }
270
271 // check size methods
272 static void checkSizes ( const Matrix &matrix )
273 {
274 [[maybe_unused]] const size_type size = matrix.size();
275 const size_type rows = MatrixSizeHelper< Matrix >::rows( matrix );
276 const size_type cols = MatrixSizeHelper< Matrix >::cols( matrix );
277 const size_type N = matrix.N();
278 const size_type M = matrix.M();
279
280 if( N != rows || M != cols )
281 DUNE_THROW( Dune::RangeError, "Returned inconsistent sizes." );
282 }
283
284 // check read-only access to data
285 static void checkDataAccess ( const Matrix &matrix )
286 {
287 const size_type size = matrix.size();
288 for( size_type i = size_type( 0 ); i < size; ++i )
289 [[maybe_unused]] const_row_reference row = matrix[ i ];
290
291 const size_type rows = MatrixSizeHelper< Matrix >::rows( matrix );
292 const size_type cols = MatrixSizeHelper< Matrix >::cols( matrix );
293 for( size_type i = size_type( 0 ); i < rows; ++i )
294 {
295 for( size_type j = size_type( 0 ); j < cols; ++j )
296 {
297 [[maybe_unused]] bool exists = matrix.exists( i, j );
298 [[maybe_unused]] const value_type &value = matrix[ i ][ j ];
299 }
300 }
301 }
302
303 // check norms
304 static void checkNorms ( const Matrix &matrix )
305 {
306 typedef typename Dune::FieldTraits< value_type >::real_type real_type;
307 real_type frobenius_norm = matrix.frobenius_norm();
308 real_type frobenius_norm2 = matrix.frobenius_norm2();
309 real_type infinity_norm = matrix.infinity_norm() ;
310 real_type infinity_norm_real = matrix.infinity_norm_real();
311
312 if( std::min( std::min( frobenius_norm, frobenius_norm2 ),
313 std::min( infinity_norm, infinity_norm_real ) ) < real_type( 0 ) )
314 DUNE_THROW( Dune::InvalidStateException, "Norms must return non-negative value." );
315 }
316
317 // check basic linear algebra methods
318 static void checkLinearAlgebra ( const Matrix &matrix )
319 {
320 typename Traits::domain_type domain = Traits::domain( matrix );
321 typename Traits::range_type range = Traits::range( matrix );
322 typename Traits::value_type alpha( 1 );
323
324 matrix.mv( domain, range );
325 matrix.mtv( range, domain );
326 matrix.umv( domain, range );
327 matrix.umtv( range, domain );
328 matrix.umhv( range, domain );
329 matrix.mmv( domain, range );
330 matrix.mmtv( range, domain );
331 matrix.mmhv( range, domain );
332 matrix.usmv( alpha, domain, range );
333 matrix.usmtv( alpha, range, domain );
334 matrix.usmhv( alpha, range, domain );
335 }
336
337 // check iterator methods
338 static void checkIterators ( const Matrix &matrix )
339 {
340 const ConstIterator end = matrix.end();
341 for( ConstIterator it = matrix.begin(); it != end; ++it )
342 [[maybe_unused]] const_row_reference row = *it;
343 }
344 };
345
346
347
348 // CheckNonConstMatrix
349 // -------------------
350
351 template< class Matrix, class Traits >
352 struct CheckNonConstMatrix
353 {
354 // required type definitions
355 typedef typename Matrix::size_type size_type;
356 typedef typename Matrix::value_type value_type;
357 typedef typename Matrix::row_reference row_reference;
358 typedef typename Matrix::row_type row_type;
359 typedef typename Matrix::Iterator Iterator;
360
361 static void apply ( Matrix &matrix )
362 {
363 checkIterators( matrix );
364 checkAssignment( matrix );
365
366 CheckIfSquareMatrix< Matrix, Traits >::apply( matrix );
367
368 // TODO: check scalar/matrix and matrix/matrix operations
369 // Matrix &operator+= ( const Matrix &other );
370 // Matrix &operator-= ( const Matrix &other );
371 // Matrix &operator*= ( const value_type &v );
372 // Matrix &operator/= ( const value_type &v );
373 // Matrix &axpy += ( const value_type &v, const Matrix &other );
374 // Matrix &axpy += ( const value_type &v, const Matrix &other );
375 // Matrix &leftmultiply ( const Matrix &other );
376 // Matrix &rightmultiply ( const Matrix &other );
377 }
378
379 // check assignment
380 static void checkAssignment ( Matrix &matrix )
381 {
382 matrix = value_type( 1 );
383
384 const size_type size = matrix.size();
385 for( size_type i = size_type( 0 ); i < size; ++i )
386 {
387 row_reference row = matrix[ i ];
388 row = row_type( value_type( 0 ) );
389 }
390
391 const size_type rows = MatrixSizeHelper< Matrix >::rows( matrix );
392 const size_type cols = MatrixSizeHelper< Matrix >::cols( matrix );
393 for( size_type i = size_type( 0 ); i < rows; ++i )
394 {
395 for( size_type j = size_type( 0 ); j < cols; ++j )
396 matrix[ i ][ j ] = ( i == j ? value_type( 1 ) : value_type( 0 ) );
397 }
398 }
399
400 // check iterator methods
401 static void checkIterators ( Matrix &matrix )
402 {
403 const Iterator end = matrix.end();
404 for( Iterator it = matrix.begin(); it != end; ++it )
405 {
406 row_reference row = *it;
407 row = row_type( value_type( 0 ) );
408 }
409 }
410 };
411
412} // namespace CheckMatrixInterface
413
414
415
416namespace Dune
417{
418
419 // checkMatrixInterface
420 // --------------------
421
422 template< class Matrix, class Traits = CheckMatrixInterface::UseDynamicVector< Matrix > >
423 void checkMatrixInterface ( const Matrix &matrix )
424 {
425 CheckMatrixInterface::CheckConstMatrix< Matrix, Traits >::apply( matrix );
426 }
427
428 template< class Matrix, class Traits = CheckMatrixInterface::UseDynamicVector< Matrix > >
429 void checkMatrixInterface ( Matrix &matrix )
430 {
431 checkMatrixInterface( const_cast< const Matrix & >( matrix ) );
432 CheckMatrixInterface::CheckNonConstMatrix< Matrix, Traits >::apply( matrix );
433 }
434
435} // namespace Dune
436
437#endif // #ifndef DUNE_COMMON_CHECK_MATRIX_INTERFACE_HH
Dune::DynamicVector
Construct a vector with a dynamic size.
Definition: dynvector.hh:59
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Dune::InvalidStateException
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:254
exceptions.hh
A few common exception classes.
dynvector.hh
This file implements a dense vector with a dynamic size.
ftraits.hh
Type traits to determine the type of reals (when working with complex numbers)
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune::Hybrid::min
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)