DUNE-FEM (unstable)

fmatrix.hh
Go to the documentation of this file.
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_FMATRIX_HH
6#define DUNE_FMATRIX_HH
7
8#include <cmath>
9#include <cstddef>
10#include <iostream>
11#include <algorithm>
12#include <initializer_list>
13
21#include <dune/common/matrixconcepts.hh>
22
23namespace Dune
24{
25
26 namespace Impl
27 {
28
29 template<class M>
30 class ColumnVectorView
31 {
32 public:
33
34 using value_type = typename M::value_type;
35 using size_type = typename M::size_type;
36
37 constexpr ColumnVectorView(M& matrix, size_type col) :
38 matrix_(matrix),
39 col_(col)
40 {}
41
42 constexpr size_type N () const {
43 return matrix_.N();
44 }
45
46 template<class M_ = M,
47 std::enable_if_t<std::is_same_v<M_,M> and not std::is_const_v<M_>, int> = 0>
48 constexpr value_type& operator[] (size_type row) {
49 return matrix_[row][col_];
50 }
51
52 constexpr const value_type& operator[] (size_type row) const {
53 return matrix_[row][col_];
54 }
55
56 protected:
57 M& matrix_;
58 const size_type col_;
59 };
60
61 }
62
63 template<typename M>
64 struct FieldTraits< Impl::ColumnVectorView<M> >
65 {
66 using field_type = typename FieldTraits<M>::field_type;
67 using real_type = typename FieldTraits<M>::real_type;
68 };
69
81 template< class K, int ROWS, int COLS = ROWS > class FieldMatrix;
82
83
84 template< class K, int ROWS, int COLS >
85 struct DenseMatVecTraits< FieldMatrix<K,ROWS,COLS> >
86 {
87 typedef FieldMatrix<K,ROWS,COLS> derived_type;
88
89 // each row is implemented by a field vector
90 typedef FieldVector<K,COLS> row_type;
91
92 typedef row_type &row_reference;
93 typedef const row_type &const_row_reference;
94
95 typedef std::array<row_type,ROWS> container_type;
96 typedef K value_type;
97 typedef typename container_type::size_type size_type;
98 };
99
100 template< class K, int ROWS, int COLS >
101 struct FieldTraits< FieldMatrix<K,ROWS,COLS> >
102 {
103 typedef typename FieldTraits<K>::field_type field_type;
104 typedef typename FieldTraits<K>::real_type real_type;
105 };
106
115 template<class K, int ROWS, int COLS>
116 class FieldMatrix : public DenseMatrix< FieldMatrix<K,ROWS,COLS> >
117 {
118 std::array< FieldVector<K,COLS>, ROWS > _data;
120 public:
121
123 constexpr static int rows = ROWS;
125 constexpr static int cols = COLS;
126
127 typedef typename Base::size_type size_type;
128 typedef typename Base::row_type row_type;
129
130 typedef typename Base::row_reference row_reference;
132
133 //===== constructors
136 constexpr FieldMatrix() = default;
137
140 constexpr FieldMatrix(std::initializer_list<Dune::FieldVector<K, cols> > const &l) {
141 assert(l.size() == rows); // Actually, this is not needed any more!
142 for(std::size_t i=0; i<std::min<std::size_t>(ROWS, l.size()); ++i)
143 _data[i] = std::data(l)[i];
144 }
145
147 FieldMatrix(const FieldMatrix&) = default;
148
150 template <class T,
151 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
152 constexpr FieldMatrix(T const& rhs)
153 : _data{}
154 {
155 *this = rhs;
156 }
157
158 using Base::operator=;
159
161 constexpr FieldMatrix& operator=(const FieldMatrix&) = default;
162
164 template<typename T>
166 {
167 _data = x._data;
168 return *this;
169 }
170
172 template <typename T, int rows, int cols>
174
177 {
179 for( int i = 0; i < ROWS; ++i )
180 for( int j = 0; j < COLS; ++j )
181 AT[j][i] = (*this)[i][j];
182 return AT;
183 }
184
186 template <class OtherScalar>
187 friend constexpr auto operator+ ( const FieldMatrix& matrixA,
189 {
191
192 for (size_type i = 0; i < ROWS; ++i)
193 for (size_type j = 0; j < COLS; ++j)
194 result[i][j] = matrixA[i][j] + matrixB[i][j];
195
196 return result;
197 }
198
200 template <class OtherScalar>
201 friend constexpr auto operator- ( const FieldMatrix& matrixA,
203 {
205
206 for (size_type i = 0; i < ROWS; ++i)
207 for (size_type j = 0; j < COLS; ++j)
208 result[i][j] = matrixA[i][j] - matrixB[i][j];
209
210 return result;
211 }
212
214 template <class Scalar,
215 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
216 friend constexpr auto operator* ( const FieldMatrix& matrix, Scalar scalar)
217 {
219
220 for (size_type i = 0; i < ROWS; ++i)
221 for (size_type j = 0; j < COLS; ++j)
222 result[i][j] = matrix[i][j] * scalar;
223
224 return result;
225 }
226
228 template <class Scalar,
229 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
230 friend constexpr auto operator* ( Scalar scalar, const FieldMatrix& matrix)
231 {
233
234 for (size_type i = 0; i < ROWS; ++i)
235 for (size_type j = 0; j < COLS; ++j)
236 result[i][j] = scalar * matrix[i][j];
237
238 return result;
239 }
240
242 template <class Scalar,
243 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
244 friend constexpr auto operator/ ( const FieldMatrix& matrix, Scalar scalar)
245 {
247
248 for (size_type i = 0; i < ROWS; ++i)
249 for (size_type j = 0; j < COLS; ++j)
250 result[i][j] = matrix[i][j] / scalar;
251
252 return result;
253 }
254
257 template <class OtherScalar, int otherCols>
258 friend constexpr auto operator* ( const FieldMatrix& matrixA,
260 {
262
263 for (size_type i = 0; i < matrixA.mat_rows(); ++i)
264 for (size_type j = 0; j < matrixB.mat_cols(); ++j)
265 {
266 result[i][j] = 0;
267 for (size_type k = 0; k < matrixA.mat_cols(); ++k)
268 result[i][j] += matrixA[i][k] * matrixB[k][j];
269 }
270
271 return result;
272 }
273
280 template <class OtherMatrix, std::enable_if_t<
281 Impl::IsStaticSizeMatrix_v<OtherMatrix>
282 and not Impl::IsFieldMatrix_v<OtherMatrix>
283 , int> = 0>
284 friend constexpr auto operator* ( const FieldMatrix& matrixA,
285 const OtherMatrix& matrixB)
286 {
287 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
289 for (std::size_t j=0; j<rows; ++j)
290 matrixB.mtv(matrixA[j], result[j]);
291 return result;
292 }
293
300 template <class OtherMatrix, std::enable_if_t<
301 Impl::IsStaticSizeMatrix_v<OtherMatrix>
302 and not Impl::IsFieldMatrix_v<OtherMatrix>
303 , int> = 0>
304 friend constexpr auto operator* ( const OtherMatrix& matrixA,
305 const FieldMatrix& matrixB)
306 {
307 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
309 for (std::size_t j=0; j<cols; ++j)
310 {
311 auto B_j = Impl::ColumnVectorView(matrixB, j);
312 auto result_j = Impl::ColumnVectorView(result, j);
313 matrixA.mv(B_j, result_j);
314 }
315 return result;
316 }
317
319 template<int l>
321 {
323
324 for (size_type i=0; i<l; i++) {
325 for (size_type j=0; j<cols; j++) {
326 C[i][j] = 0;
327 for (size_type k=0; k<rows; k++)
328 C[i][j] += M[i][k]*(*this)[k][j];
329 }
330 }
331 return C;
332 }
333
335
337 template <int r, int c>
339 {
340 static_assert(r == c, "Cannot rightmultiply with non-square matrix");
341 static_assert(r == cols, "Size mismatch");
343
344 for (size_type i=0; i<rows; i++)
345 for (size_type j=0; j<cols; j++) {
346 (*this)[i][j] = 0;
347 for (size_type k=0; k<cols; k++)
348 (*this)[i][j] += C[i][k]*M[k][j];
349 }
350 return *this;
351 }
352
354 template<int l>
356 {
358
359 for (size_type i=0; i<rows; i++) {
360 for (size_type j=0; j<l; j++) {
361 C[i][j] = 0;
362 for (size_type k=0; k<cols; k++)
363 C[i][j] += (*this)[i][k]*M[k][j];
364 }
365 }
366 return C;
367 }
368
369 // make this thing a matrix
370 static constexpr size_type mat_rows() { return ROWS; }
371 static constexpr size_type mat_cols() { return COLS; }
372
373 constexpr row_reference mat_access ( size_type i )
374 {
375 DUNE_ASSERT_BOUNDS(i < ROWS);
376 return _data[i];
377 }
378
379 constexpr const_row_reference mat_access ( size_type i ) const
380 {
381 DUNE_ASSERT_BOUNDS(i < ROWS);
382 return _data[i];
383 }
384 };
385
386#ifndef DOXYGEN // hide specialization
389 template<class K>
390 class FieldMatrix<K,1,1> : public DenseMatrix< FieldMatrix<K,1,1> >
391 {
392 FieldVector<K,1> _data;
393 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
394 public:
395 // standard constructor and everything is sufficient ...
396
397 //===== type definitions and constants
398
400 typedef typename Base::size_type size_type;
401
404 constexpr static int blocklevel = 1;
405
406 typedef typename Base::row_type row_type;
407
408 typedef typename Base::row_reference row_reference;
409 typedef typename Base::const_row_reference const_row_reference;
410
413 constexpr static int rows = 1;
416 constexpr static int cols = 1;
417
418 //===== constructors
421 constexpr FieldMatrix() = default;
422
425 FieldMatrix(std::initializer_list<Dune::FieldVector<K, 1>> const &l)
426 {
427 std::copy_n(l.begin(), std::min<std::size_t>(1, l.size()), &_data);
428 }
429
430 template <class T,
431 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
432 constexpr FieldMatrix(T const& rhs)
433 {
434 *this = rhs;
435 }
436
437 using Base::operator=;
438
440 constexpr FieldMatrix<K, 1, 1> transposed() const
441 {
442 return *this;
443 }
444
446 template <class OtherScalar>
447 friend constexpr auto operator+ ( const FieldMatrix& matrixA,
448 const FieldMatrix<OtherScalar,1,1>& matrixB)
449 {
450 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] + matrixB[0][0]};
451 }
452
454 template <class Scalar,
455 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
456 friend constexpr auto operator+ ( const FieldMatrix& matrix,
457 const Scalar& scalar)
458 {
459 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
460 }
461
463 template <class Scalar,
464 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
465 friend constexpr auto operator+ ( const Scalar& scalar,
466 const FieldMatrix& matrix)
467 {
468 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar + matrix[0][0]};
469 }
470
472 template <class OtherScalar>
473 friend constexpr auto operator- ( const FieldMatrix& matrixA,
474 const FieldMatrix<OtherScalar,1,1>& matrixB)
475 {
476 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] - matrixB[0][0]};
477 }
478
480 template <class Scalar,
481 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
482 friend constexpr auto operator- ( const FieldMatrix& matrix,
483 const Scalar& scalar)
484 {
485 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
486 }
487
489 template <class Scalar,
490 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
491 friend constexpr auto operator- ( const Scalar& scalar,
492 const FieldMatrix& matrix)
493 {
494 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar - matrix[0][0]};
495 }
496
498 template <class Scalar,
499 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
500 friend constexpr auto operator* ( const FieldMatrix& matrix, Scalar scalar)
501 {
502 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] * scalar};
503 }
504
506 template <class Scalar,
507 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
508 friend constexpr auto operator* ( Scalar scalar, const FieldMatrix& matrix)
509 {
510 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {scalar * matrix[0][0]};
511 }
512
514 template <class Scalar,
515 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
516 friend constexpr auto operator/ ( const FieldMatrix& matrix, Scalar scalar)
517 {
518 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] / scalar};
519 }
520
521 //===== solve
522
525 template <class OtherScalar, int otherCols>
526 friend constexpr auto operator* ( const FieldMatrix& matrixA,
527 const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
528 {
529 FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
530
531 for (size_type j = 0; j < matrixB.mat_cols(); ++j)
532 result[0][j] = matrixA[0][0] * matrixB[0][j];
533
534 return result;
535 }
536
543 template <class OtherMatrix, std::enable_if_t<
544 Impl::IsStaticSizeMatrix_v<OtherMatrix>
545 and not Impl::IsFieldMatrix_v<OtherMatrix>
546 and (OtherMatrix::rows==1)
547 , int> = 0>
548 friend constexpr auto operator* ( const FieldMatrix& matrixA,
549 const OtherMatrix& matrixB)
550 {
551 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
553 for (std::size_t j=0; j<rows; ++j)
554 matrixB.mtv(matrixA[j], result[j]);
555 return result;
556 }
557
564 template <class OtherMatrix, std::enable_if_t<
565 Impl::IsStaticSizeMatrix_v<OtherMatrix>
566 and not Impl::IsFieldMatrix_v<OtherMatrix>
567 and (OtherMatrix::cols==1)
568 , int> = 0>
569 friend constexpr auto operator* ( const OtherMatrix& matrixA,
570 const FieldMatrix& matrixB)
571 {
572 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
574 for (std::size_t j=0; j<cols; ++j)
575 {
576 auto B_j = Impl::ColumnVectorView(matrixB, j);
577 auto result_j = Impl::ColumnVectorView(result, j);
578 matrixA.mv(B_j, result_j);
579 }
580 return result;
581 }
582
584 template<int l>
585 constexpr FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
586 {
587 FieldMatrix<K,l,1> C;
588 for (size_type j=0; j<l; j++)
589 C[j][0] = M[j][0]*(*this)[0][0];
590 return C;
591 }
592
594 constexpr FieldMatrix& rightmultiply (const FieldMatrix& M)
595 {
596 _data[0] *= M[0][0];
597 return *this;
598 }
599
601 template<int l>
602 constexpr FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
603 {
604 FieldMatrix<K,1,l> C;
605
606 for (size_type j=0; j<l; j++)
607 C[0][j] = M[0][j]*_data[0];
608 return C;
609 }
610
611 // make this thing a matrix
612 static constexpr size_type mat_rows() { return 1; }
613 static constexpr size_type mat_cols() { return 1; }
614
615 constexpr row_reference mat_access ([[maybe_unused]] size_type i)
616 {
617 DUNE_ASSERT_BOUNDS(i == 0);
618 return _data;
619 }
620
621 constexpr const_row_reference mat_access ([[maybe_unused]] size_type i) const
622 {
623 DUNE_ASSERT_BOUNDS(i == 0);
624 return _data;
625 }
626
628 constexpr FieldMatrix& operator+= (const K& k)
629 {
630 _data[0] += k;
631 return (*this);
632 }
633
635 constexpr FieldMatrix& operator-= (const K& k)
636 {
637 _data[0] -= k;
638 return (*this);
639 }
640
642 constexpr FieldMatrix& operator*= (const K& k)
643 {
644 _data[0] *= k;
645 return (*this);
646 }
647
649 constexpr FieldMatrix& operator/= (const K& k)
650 {
651 _data[0] /= k;
652 return (*this);
653 }
654
655 //===== conversion operator
656
657 constexpr operator const K& () const { return _data[0]; }
658
659 };
660
662 template<typename K>
663 std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
664 {
665 s << a[0][0];
666 return s;
667 }
668
669#endif // DOXYGEN
670
671 namespace FMatrixHelp {
672
674 template <typename K>
675 static constexpr K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
676 {
677 using real_type = typename FieldTraits<K>::real_type;
678 inverse[0][0] = real_type(1.0)/matrix[0][0];
679 return matrix[0][0];
680 }
681
683 template <typename K>
684 static constexpr K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
685 {
686 return invertMatrix(matrix,inverse);
687 }
688
689
691 template <typename K>
692 static constexpr K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
693 {
694 using real_type = typename FieldTraits<K>::real_type;
695 // code generated by maple
696 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
697 K det_1 = real_type(1.0)/det;
698 inverse[0][0] = matrix[1][1] * det_1;
699 inverse[0][1] = - matrix[0][1] * det_1;
700 inverse[1][0] = - matrix[1][0] * det_1;
701 inverse[1][1] = matrix[0][0] * det_1;
702 return det;
703 }
704
707 template <typename K>
708 static constexpr K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
709 {
710 using real_type = typename FieldTraits<K>::real_type;
711 // code generated by maple
712 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
713 K det_1 = real_type(1.0)/det;
714 inverse[0][0] = matrix[1][1] * det_1;
715 inverse[1][0] = - matrix[0][1] * det_1;
716 inverse[0][1] = - matrix[1][0] * det_1;
717 inverse[1][1] = matrix[0][0] * det_1;
718 return det;
719 }
720
722 template <typename K>
723 static constexpr K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
724 {
725 using real_type = typename FieldTraits<K>::real_type;
726 // code generated by maple
727 K t4 = matrix[0][0] * matrix[1][1];
728 K t6 = matrix[0][0] * matrix[1][2];
729 K t8 = matrix[0][1] * matrix[1][0];
730 K t10 = matrix[0][2] * matrix[1][0];
731 K t12 = matrix[0][1] * matrix[2][0];
732 K t14 = matrix[0][2] * matrix[2][0];
733
734 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
735 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
736 K t17 = real_type(1.0)/det;
737
738 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
739 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
740 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
741 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
742 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
743 inverse[1][2] = -(t6-t10) * t17;
744 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
745 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
746 inverse[2][2] = (t4-t8) * t17;
747
748 return det;
749 }
750
752 template <typename K>
753 static constexpr K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
754 {
755 using real_type = typename FieldTraits<K>::real_type;
756 // code generated by maple
757 K t4 = matrix[0][0] * matrix[1][1];
758 K t6 = matrix[0][0] * matrix[1][2];
759 K t8 = matrix[0][1] * matrix[1][0];
760 K t10 = matrix[0][2] * matrix[1][0];
761 K t12 = matrix[0][1] * matrix[2][0];
762 K t14 = matrix[0][2] * matrix[2][0];
763
764 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
765 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
766 K t17 = real_type(1.0)/det;
767
768 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
769 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
770 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
771 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
772 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
773 inverse[2][1] = -(t6-t10) * t17;
774 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
775 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
776 inverse[2][2] = (t4-t8) * t17;
777
778 return det;
779 }
780
782 template< class K, int m, int n, int p >
783 static constexpr void multMatrix ( const FieldMatrix< K, m, n > &A,
784 const FieldMatrix< K, n, p > &B,
786 {
787 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
788
789 for( size_type i = 0; i < m; ++i )
790 {
791 for( size_type j = 0; j < p; ++j )
792 {
793 ret[ i ][ j ] = K( 0 );
794 for( size_type k = 0; k < n; ++k )
795 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
796 }
797 }
798 }
799
801 template <typename K, int rows, int cols>
803 {
804 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
805
806 for(size_type i=0; i<cols; i++)
807 for(size_type j=0; j<cols; j++)
808 {
809 ret[i][j]=0.0;
810 for(size_type k=0; k<rows; k++)
811 ret[i][j]+=matrix[k][i]*matrix[k][j];
812 }
813 }
814
815 using Dune::DenseMatrixHelp::multAssign;
816
818 template <typename K, int rows, int cols>
819 static constexpr void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
820 {
821 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
822
823 for(size_type i=0; i<cols; ++i)
824 {
825 ret[i] = 0.0;
826 for(size_type j=0; j<rows; ++j)
827 ret[i] += matrix[j][i]*x[j];
828 }
829 }
830
832 template <typename K, int rows, int cols>
833 static constexpr FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
834 {
836 multAssign(matrix,x,ret);
837 return ret;
838 }
839
841 template <typename K, int rows, int cols>
843 {
845 multAssignTransposed( matrix, x, ret );
846 return ret;
847 }
848
849 } // end namespace FMatrixHelp
850
853} // end namespace
854
855#include "fmatrixev.hh"
856#endif
Macro for wrapping boundary checks.
A dense n x m matrix.
Definition: densematrix.hh:145
constexpr derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition: densematrix.hh:294
constexpr derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:326
constexpr derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition: densematrix.hh:317
constexpr size_type M() const
number of columns
Definition: densematrix.hh:708
FieldMatrix< K, ROWS, COLS > & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:650
constexpr derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:334
constexpr derived_type operator-() const
Matrix negation.
Definition: densematrix.hh:303
constexpr void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:392
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: densematrix.hh:183
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:174
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:171
Traits::const_row_reference const_row_reference
The type used to represent a reference to a constant row (usually const row_type &)
Definition: densematrix.hh:180
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:177
A dense n x m matrix.
Definition: fmatrix.hh:117
constexpr FieldMatrix()=default
Default constructor.
constexpr FieldMatrix< K, rows, l > rightmultiplyany(const FieldMatrix< K, cols, l > &M) const
Multiplies M from the right to this matrix, this matrix is not modified.
Definition: fmatrix.hh:355
constexpr FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:338
friend constexpr auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition: fmatrix.hh:216
constexpr FieldMatrix & operator=(const FieldMatrix< T, ROWS, COLS > &x)
copy assignment from FieldMatrix over a different field
Definition: fmatrix.hh:165
constexpr FieldMatrix(T const &rhs)
copy constructor from assignable type T
Definition: fmatrix.hh:152
FieldMatrix & operator=(FieldMatrix< T, rows, cols > const &)=delete
no copy assignment from FieldMatrix of different size
constexpr FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition: fmatrix.hh:140
static constexpr int rows
The number of rows.
Definition: fmatrix.hh:123
constexpr FieldMatrix & operator=(const FieldMatrix &)=default
copy assignment operator
static constexpr int cols
The number of columns.
Definition: fmatrix.hh:125
constexpr FieldMatrix< K, COLS, ROWS > transposed() const
Return transposed of the matrix as FieldMatrix.
Definition: fmatrix.hh:176
friend constexpr auto operator/(const FieldMatrix &matrix, Scalar scalar)
vector space division by scalar
Definition: fmatrix.hh:244
friend constexpr auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition: fmatrix.hh:187
constexpr FieldMatrix< K, l, cols > leftmultiplyany(const FieldMatrix< K, l, rows > &M) const
Multiplies M from the left to this matrix, this matrix is not modified.
Definition: fmatrix.hh:320
FieldMatrix(const FieldMatrix &)=default
copy constructor
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
A few common exception classes.
static constexpr void multAssignTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x, FieldVector< K, cols > &ret)
calculates ret = matrix^T * x
Definition: fmatrix.hh:819
static constexpr FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:842
static constexpr K invertMatrix_retTransposed(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:684
static constexpr void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:802
static constexpr void multMatrix(const FieldMatrix< K, m, n > &A, const FieldMatrix< K, n, p > &B, FieldMatrix< K, m, p > &ret)
calculates ret = A * B
Definition: fmatrix.hh:783
static constexpr K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:675
static constexpr FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:833
Eigenvalue computations for the FieldMatrix class.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:30
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
Dune namespace.
Definition: alignedallocator.hh:13
Various precision settings for calculations with FieldMatrix and FieldVector.
Compute type of the result of an arithmetic operation involving two different number types.
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Feb 17, 23:30, 2025)