Dune Core Modules (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 template<class,int,int> friend class FieldMatrix;
119 std::array< FieldVector<K,COLS>, ROWS > _data;
121 public:
122
124 constexpr static int rows = ROWS;
126 constexpr static int cols = COLS;
127
128 typedef typename Base::size_type size_type;
129 typedef typename Base::row_type row_type;
130
131 typedef typename Base::row_reference row_reference;
133
134 //===== constructors
137 constexpr FieldMatrix() = default;
138
141 constexpr FieldMatrix(std::initializer_list<Dune::FieldVector<K, cols> > const &l) {
142 assert(l.size() == rows); // Actually, this is not needed any more!
143 for(std::size_t i=0; i<std::min<std::size_t>(ROWS, l.size()); ++i)
144 _data[i] = std::data(l)[i];
145 }
146
148 FieldMatrix(const FieldMatrix&) = default;
149
151 template <class T,
152 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
153 constexpr FieldMatrix(T const& rhs)
154 : _data{}
155 {
156 *this = rhs;
157 }
158
159 using Base::operator=;
160
162 constexpr FieldMatrix& operator=(const FieldMatrix&) = default;
163
165 template<typename T>
167 {
168 // The copy must be done element-by-element since a std::array does not have
169 // a converting assignment operator.
170 for (std::size_t i = 0; i < _data.size(); ++i)
171 _data[i] = x._data[i];
172 return *this;
173 }
174
176 template <typename T, int rows, int cols>
178
181 {
183 for( int i = 0; i < ROWS; ++i )
184 for( int j = 0; j < COLS; ++j )
185 AT[j][i] = (*this)[i][j];
186 return AT;
187 }
188
190 template <class OtherScalar>
191 friend constexpr auto operator+ ( const FieldMatrix& matrixA,
193 {
195
196 for (size_type i = 0; i < ROWS; ++i)
197 for (size_type j = 0; j < COLS; ++j)
198 result[i][j] = matrixA[i][j] + matrixB[i][j];
199
200 return result;
201 }
202
204 template <class OtherScalar>
205 friend constexpr auto operator- ( const FieldMatrix& matrixA,
207 {
209
210 for (size_type i = 0; i < ROWS; ++i)
211 for (size_type j = 0; j < COLS; ++j)
212 result[i][j] = matrixA[i][j] - matrixB[i][j];
213
214 return result;
215 }
216
218 template <class Scalar,
219 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
220 friend constexpr auto operator* ( const FieldMatrix& matrix, Scalar scalar)
221 {
223
224 for (size_type i = 0; i < ROWS; ++i)
225 for (size_type j = 0; j < COLS; ++j)
226 result[i][j] = matrix[i][j] * scalar;
227
228 return result;
229 }
230
232 template <class Scalar,
233 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
234 friend constexpr auto operator* ( Scalar scalar, const FieldMatrix& matrix)
235 {
237
238 for (size_type i = 0; i < ROWS; ++i)
239 for (size_type j = 0; j < COLS; ++j)
240 result[i][j] = scalar * matrix[i][j];
241
242 return result;
243 }
244
246 template <class Scalar,
247 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
248 friend constexpr auto operator/ ( const FieldMatrix& matrix, Scalar scalar)
249 {
251
252 for (size_type i = 0; i < ROWS; ++i)
253 for (size_type j = 0; j < COLS; ++j)
254 result[i][j] = matrix[i][j] / scalar;
255
256 return result;
257 }
258
261 template <class OtherScalar, int otherCols>
262 friend constexpr auto operator* ( const FieldMatrix& matrixA,
264 {
266
267 for (size_type i = 0; i < matrixA.mat_rows(); ++i)
268 for (size_type j = 0; j < matrixB.mat_cols(); ++j)
269 {
270 result[i][j] = 0;
271 for (size_type k = 0; k < matrixA.mat_cols(); ++k)
272 result[i][j] += matrixA[i][k] * matrixB[k][j];
273 }
274
275 return result;
276 }
277
284 template <class OtherMatrix, std::enable_if_t<
285 Impl::IsStaticSizeMatrix_v<OtherMatrix>
286 and not Impl::IsFieldMatrix_v<OtherMatrix>
287 , int> = 0>
288 friend constexpr auto operator* ( const FieldMatrix& matrixA,
289 const OtherMatrix& matrixB)
290 {
291 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
293 for (std::size_t j=0; j<rows; ++j)
294 matrixB.mtv(matrixA[j], result[j]);
295 return result;
296 }
297
304 template <class OtherMatrix, std::enable_if_t<
305 Impl::IsStaticSizeMatrix_v<OtherMatrix>
306 and not Impl::IsFieldMatrix_v<OtherMatrix>
307 , int> = 0>
308 friend constexpr auto operator* ( const OtherMatrix& matrixA,
309 const FieldMatrix& matrixB)
310 {
311 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
313 for (std::size_t j=0; j<cols; ++j)
314 {
315 auto B_j = Impl::ColumnVectorView(matrixB, j);
316 auto result_j = Impl::ColumnVectorView(result, j);
317 matrixA.mv(B_j, result_j);
318 }
319 return result;
320 }
321
323 template<int l>
325 {
327
328 for (size_type i=0; i<l; i++) {
329 for (size_type j=0; j<cols; j++) {
330 C[i][j] = 0;
331 for (size_type k=0; k<rows; k++)
332 C[i][j] += M[i][k]*(*this)[k][j];
333 }
334 }
335 return C;
336 }
337
339
341 template <int r, int c>
343 {
344 static_assert(r == c, "Cannot rightmultiply with non-square matrix");
345 static_assert(r == cols, "Size mismatch");
347
348 for (size_type i=0; i<rows; i++)
349 for (size_type j=0; j<cols; j++) {
350 (*this)[i][j] = 0;
351 for (size_type k=0; k<cols; k++)
352 (*this)[i][j] += C[i][k]*M[k][j];
353 }
354 return *this;
355 }
356
358 template<int l>
360 {
362
363 for (size_type i=0; i<rows; i++) {
364 for (size_type j=0; j<l; j++) {
365 C[i][j] = 0;
366 for (size_type k=0; k<cols; k++)
367 C[i][j] += (*this)[i][k]*M[k][j];
368 }
369 }
370 return C;
371 }
372
373 // make this thing a matrix
374 static constexpr size_type mat_rows() { return ROWS; }
375 static constexpr size_type mat_cols() { return COLS; }
376
377 constexpr row_reference mat_access ( size_type i )
378 {
379 DUNE_ASSERT_BOUNDS(i < ROWS);
380 return _data[i];
381 }
382
383 constexpr const_row_reference mat_access ( size_type i ) const
384 {
385 DUNE_ASSERT_BOUNDS(i < ROWS);
386 return _data[i];
387 }
388 };
389
390#ifndef DOXYGEN // hide specialization
393 template<class K>
394 class FieldMatrix<K,1,1> : public DenseMatrix< FieldMatrix<K,1,1> >
395 {
396 FieldVector<K,1> _data;
397 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
398 public:
399 // standard constructor and everything is sufficient ...
400
401 //===== type definitions and constants
402
404 typedef typename Base::size_type size_type;
405
408 constexpr static int blocklevel = 1;
409
410 typedef typename Base::row_type row_type;
411
412 typedef typename Base::row_reference row_reference;
413 typedef typename Base::const_row_reference const_row_reference;
414
417 constexpr static int rows = 1;
420 constexpr static int cols = 1;
421
422 //===== constructors
425 constexpr FieldMatrix() = default;
426
429 FieldMatrix(std::initializer_list<Dune::FieldVector<K, 1>> const &l)
430 {
431 std::copy_n(l.begin(), std::min<std::size_t>(1, l.size()), &_data);
432 }
433
434 template <class T,
435 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
436 constexpr FieldMatrix(T const& rhs)
437 {
438 *this = rhs;
439 }
440
441 using Base::operator=;
442
444 constexpr FieldMatrix<K, 1, 1> transposed() const
445 {
446 return *this;
447 }
448
450 template <class OtherScalar>
451 friend constexpr auto operator+ ( const FieldMatrix& matrixA,
452 const FieldMatrix<OtherScalar,1,1>& matrixB)
453 {
454 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] + matrixB[0][0]};
455 }
456
458 template <class Scalar,
459 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
460 friend constexpr auto operator+ ( const FieldMatrix& matrix,
461 const Scalar& scalar)
462 {
463 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
464 }
465
467 template <class Scalar,
468 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
469 friend constexpr auto operator+ ( const Scalar& scalar,
470 const FieldMatrix& matrix)
471 {
472 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar + matrix[0][0]};
473 }
474
476 template <class OtherScalar>
477 friend constexpr auto operator- ( const FieldMatrix& matrixA,
478 const FieldMatrix<OtherScalar,1,1>& matrixB)
479 {
480 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] - matrixB[0][0]};
481 }
482
484 template <class Scalar,
485 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
486 friend constexpr auto operator- ( const FieldMatrix& matrix,
487 const Scalar& scalar)
488 {
489 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
490 }
491
493 template <class Scalar,
494 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
495 friend constexpr auto operator- ( const Scalar& scalar,
496 const FieldMatrix& matrix)
497 {
498 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar - matrix[0][0]};
499 }
500
502 template <class Scalar,
503 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
504 friend constexpr auto operator* ( const FieldMatrix& matrix, Scalar scalar)
505 {
506 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] * scalar};
507 }
508
510 template <class Scalar,
511 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
512 friend constexpr auto operator* ( Scalar scalar, const FieldMatrix& matrix)
513 {
514 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {scalar * matrix[0][0]};
515 }
516
518 template <class Scalar,
519 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
520 friend constexpr auto operator/ ( const FieldMatrix& matrix, Scalar scalar)
521 {
522 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] / scalar};
523 }
524
525 //===== solve
526
529 template <class OtherScalar, int otherCols>
530 friend constexpr auto operator* ( const FieldMatrix& matrixA,
531 const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
532 {
533 FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
534
535 for (size_type j = 0; j < matrixB.mat_cols(); ++j)
536 result[0][j] = matrixA[0][0] * matrixB[0][j];
537
538 return result;
539 }
540
547 template <class OtherMatrix, std::enable_if_t<
548 Impl::IsStaticSizeMatrix_v<OtherMatrix>
549 and not Impl::IsFieldMatrix_v<OtherMatrix>
550 and (OtherMatrix::rows==1)
551 , int> = 0>
552 friend constexpr auto operator* ( const FieldMatrix& matrixA,
553 const OtherMatrix& matrixB)
554 {
555 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
557 for (std::size_t j=0; j<rows; ++j)
558 matrixB.mtv(matrixA[j], result[j]);
559 return result;
560 }
561
568 template <class OtherMatrix, std::enable_if_t<
569 Impl::IsStaticSizeMatrix_v<OtherMatrix>
570 and not Impl::IsFieldMatrix_v<OtherMatrix>
571 and (OtherMatrix::cols==1)
572 , int> = 0>
573 friend constexpr auto operator* ( const OtherMatrix& matrixA,
574 const FieldMatrix& matrixB)
575 {
576 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
578 for (std::size_t j=0; j<cols; ++j)
579 {
580 auto B_j = Impl::ColumnVectorView(matrixB, j);
581 auto result_j = Impl::ColumnVectorView(result, j);
582 matrixA.mv(B_j, result_j);
583 }
584 return result;
585 }
586
588 template<int l>
589 constexpr FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
590 {
591 FieldMatrix<K,l,1> C;
592 for (size_type j=0; j<l; j++)
593 C[j][0] = M[j][0]*(*this)[0][0];
594 return C;
595 }
596
598 constexpr FieldMatrix& rightmultiply (const FieldMatrix& M)
599 {
600 _data[0] *= M[0][0];
601 return *this;
602 }
603
605 template<int l>
606 constexpr FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
607 {
608 FieldMatrix<K,1,l> C;
609
610 for (size_type j=0; j<l; j++)
611 C[0][j] = M[0][j]*_data[0];
612 return C;
613 }
614
615 // make this thing a matrix
616 static constexpr size_type mat_rows() { return 1; }
617 static constexpr size_type mat_cols() { return 1; }
618
619 constexpr row_reference mat_access ([[maybe_unused]] size_type i)
620 {
621 DUNE_ASSERT_BOUNDS(i == 0);
622 return _data;
623 }
624
625 constexpr const_row_reference mat_access ([[maybe_unused]] size_type i) const
626 {
627 DUNE_ASSERT_BOUNDS(i == 0);
628 return _data;
629 }
630
632 constexpr FieldMatrix& operator+= (const K& k)
633 {
634 _data[0] += k;
635 return (*this);
636 }
637
639 constexpr FieldMatrix& operator-= (const K& k)
640 {
641 _data[0] -= k;
642 return (*this);
643 }
644
646 constexpr FieldMatrix& operator*= (const K& k)
647 {
648 _data[0] *= k;
649 return (*this);
650 }
651
653 constexpr FieldMatrix& operator/= (const K& k)
654 {
655 _data[0] /= k;
656 return (*this);
657 }
658
659 //===== conversion operator
660
661 constexpr operator const K& () const { return _data[0]; }
662
663 };
664
666 template<typename K>
667 std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
668 {
669 s << a[0][0];
670 return s;
671 }
672
673#endif // DOXYGEN
674
675 namespace FMatrixHelp {
676
678 template <typename K>
679 static constexpr K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
680 {
681 using real_type = typename FieldTraits<K>::real_type;
682 inverse[0][0] = real_type(1.0)/matrix[0][0];
683 return matrix[0][0];
684 }
685
687 template <typename K>
688 static constexpr K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
689 {
690 return invertMatrix(matrix,inverse);
691 }
692
693
695 template <typename K>
696 static constexpr K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
697 {
698 using real_type = typename FieldTraits<K>::real_type;
699 // code generated by maple
700 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
701 K det_1 = real_type(1.0)/det;
702 inverse[0][0] = matrix[1][1] * det_1;
703 inverse[0][1] = - matrix[0][1] * det_1;
704 inverse[1][0] = - matrix[1][0] * det_1;
705 inverse[1][1] = matrix[0][0] * det_1;
706 return det;
707 }
708
711 template <typename K>
712 static constexpr K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
713 {
714 using real_type = typename FieldTraits<K>::real_type;
715 // code generated by maple
716 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
717 K det_1 = real_type(1.0)/det;
718 inverse[0][0] = matrix[1][1] * det_1;
719 inverse[1][0] = - matrix[0][1] * det_1;
720 inverse[0][1] = - matrix[1][0] * det_1;
721 inverse[1][1] = matrix[0][0] * det_1;
722 return det;
723 }
724
726 template <typename K>
727 static constexpr K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
728 {
729 using real_type = typename FieldTraits<K>::real_type;
730 // code generated by maple
731 K t4 = matrix[0][0] * matrix[1][1];
732 K t6 = matrix[0][0] * matrix[1][2];
733 K t8 = matrix[0][1] * matrix[1][0];
734 K t10 = matrix[0][2] * matrix[1][0];
735 K t12 = matrix[0][1] * matrix[2][0];
736 K t14 = matrix[0][2] * matrix[2][0];
737
738 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
739 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
740 K t17 = real_type(1.0)/det;
741
742 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
743 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
744 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
745 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
746 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
747 inverse[1][2] = -(t6-t10) * t17;
748 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
749 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
750 inverse[2][2] = (t4-t8) * t17;
751
752 return det;
753 }
754
756 template <typename K>
757 static constexpr K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
758 {
759 using real_type = typename FieldTraits<K>::real_type;
760 // code generated by maple
761 K t4 = matrix[0][0] * matrix[1][1];
762 K t6 = matrix[0][0] * matrix[1][2];
763 K t8 = matrix[0][1] * matrix[1][0];
764 K t10 = matrix[0][2] * matrix[1][0];
765 K t12 = matrix[0][1] * matrix[2][0];
766 K t14 = matrix[0][2] * matrix[2][0];
767
768 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
769 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
770 K t17 = real_type(1.0)/det;
771
772 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
773 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
774 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
775 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
776 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
777 inverse[2][1] = -(t6-t10) * t17;
778 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
779 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
780 inverse[2][2] = (t4-t8) * t17;
781
782 return det;
783 }
784
786 template< class K, int m, int n, int p >
787 static constexpr void multMatrix ( const FieldMatrix< K, m, n > &A,
788 const FieldMatrix< K, n, p > &B,
790 {
791 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
792
793 for( size_type i = 0; i < m; ++i )
794 {
795 for( size_type j = 0; j < p; ++j )
796 {
797 ret[ i ][ j ] = K( 0 );
798 for( size_type k = 0; k < n; ++k )
799 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
800 }
801 }
802 }
803
805 template <typename K, int rows, int cols>
807 {
808 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
809
810 for(size_type i=0; i<cols; i++)
811 for(size_type j=0; j<cols; j++)
812 {
813 ret[i][j]=0.0;
814 for(size_type k=0; k<rows; k++)
815 ret[i][j]+=matrix[k][i]*matrix[k][j];
816 }
817 }
818
819 using Dune::DenseMatrixHelp::multAssign;
820
822 template <typename K, int rows, int cols>
823 static constexpr void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
824 {
825 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
826
827 for(size_type i=0; i<cols; ++i)
828 {
829 ret[i] = 0.0;
830 for(size_type j=0; j<rows; ++j)
831 ret[i] += matrix[j][i]*x[j];
832 }
833 }
834
836 template <typename K, int rows, int cols>
837 static constexpr FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
838 {
840 multAssign(matrix,x,ret);
841 return ret;
842 }
843
845 template <typename K, int rows, int cols>
847 {
849 multAssignTransposed( matrix, x, ret );
850 return ret;
851 }
852
853 } // end namespace FMatrixHelp
854
857} // end namespace
858
859#include "fmatrixev.hh"
860#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:359
constexpr FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:342
friend constexpr auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition: fmatrix.hh:220
constexpr FieldMatrix & operator=(const FieldMatrix< T, ROWS, COLS > &x)
copy assignment from FieldMatrix over a different field
Definition: fmatrix.hh:166
constexpr FieldMatrix(T const &rhs)
copy constructor from assignable type T
Definition: fmatrix.hh:153
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:141
static constexpr int rows
The number of rows.
Definition: fmatrix.hh:124
constexpr FieldMatrix & operator=(const FieldMatrix &)=default
copy assignment operator
static constexpr int cols
The number of columns.
Definition: fmatrix.hh:126
constexpr FieldMatrix< K, COLS, ROWS > transposed() const
Return transposed of the matrix as FieldMatrix.
Definition: fmatrix.hh:180
friend constexpr auto operator/(const FieldMatrix &matrix, Scalar scalar)
vector space division by scalar
Definition: fmatrix.hh:248
friend constexpr auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition: fmatrix.hh:191
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:324
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:823
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:846
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:688
static constexpr void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:806
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:787
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:679
static constexpr FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:837
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 (Sep 21, 22:40, 2025)