DUNE PDELab (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.
Traits for type conversions and type information.
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 28, 22:46, 2025)