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
146 template <class T,
147 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
148 FieldMatrix(T const& rhs)
149 {
150 *this = rhs;
151 }
152
153 using Base::operator=;
154
157
159 template<typename T>
161 {
162 _data = x._data;
163 return *this;
164 }
165
167 template <typename T, int rows, int cols>
169
172 {
174 for( int i = 0; i < ROWS; ++i )
175 for( int j = 0; j < COLS; ++j )
176 AT[j][i] = (*this)[i][j];
177 return AT;
178 }
179
181 template <class OtherScalar>
182 friend auto operator+ ( const FieldMatrix& matrixA,
184 {
186
187 for (size_type i = 0; i < ROWS; ++i)
188 for (size_type j = 0; j < COLS; ++j)
189 result[i][j] = matrixA[i][j] + matrixB[i][j];
190
191 return result;
192 }
193
195 template <class OtherScalar>
196 friend auto operator- ( const FieldMatrix& matrixA,
198 {
200
201 for (size_type i = 0; i < ROWS; ++i)
202 for (size_type j = 0; j < COLS; ++j)
203 result[i][j] = matrixA[i][j] - matrixB[i][j];
204
205 return result;
206 }
207
209 template <class Scalar,
210 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
211 friend auto operator* ( const FieldMatrix& matrix, Scalar scalar)
212 {
214
215 for (size_type i = 0; i < ROWS; ++i)
216 for (size_type j = 0; j < COLS; ++j)
217 result[i][j] = matrix[i][j] * scalar;
218
219 return result;
220 }
221
223 template <class Scalar,
224 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
225 friend auto operator* ( Scalar scalar, const FieldMatrix& matrix)
226 {
228
229 for (size_type i = 0; i < ROWS; ++i)
230 for (size_type j = 0; j < COLS; ++j)
231 result[i][j] = scalar * matrix[i][j];
232
233 return result;
234 }
235
237 template <class Scalar,
238 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
239 friend auto operator/ ( const FieldMatrix& matrix, Scalar scalar)
240 {
242
243 for (size_type i = 0; i < ROWS; ++i)
244 for (size_type j = 0; j < COLS; ++j)
245 result[i][j] = matrix[i][j] / scalar;
246
247 return result;
248 }
249
252 template <class OtherScalar, int otherCols>
253 friend auto operator* ( const FieldMatrix& matrixA,
255 {
257
258 for (size_type i = 0; i < matrixA.mat_rows(); ++i)
259 for (size_type j = 0; j < matrixB.mat_cols(); ++j)
260 {
261 result[i][j] = 0;
262 for (size_type k = 0; k < matrixA.mat_cols(); ++k)
263 result[i][j] += matrixA[i][k] * matrixB[k][j];
264 }
265
266 return result;
267 }
268
275 template <class OtherMatrix, std::enable_if_t<
276 Impl::IsStaticSizeMatrix_v<OtherMatrix>
277 and not Impl::IsFieldMatrix_v<OtherMatrix>
278 , int> = 0>
279 friend auto operator* ( const FieldMatrix& matrixA,
280 const OtherMatrix& matrixB)
281 {
282 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
284 for (std::size_t j=0; j<rows; ++j)
285 matrixB.mtv(matrixA[j], result[j]);
286 return result;
287 }
288
295 template <class OtherMatrix, std::enable_if_t<
296 Impl::IsStaticSizeMatrix_v<OtherMatrix>
297 and not Impl::IsFieldMatrix_v<OtherMatrix>
298 , int> = 0>
299 friend auto operator* ( const OtherMatrix& matrixA,
300 const FieldMatrix& matrixB)
301 {
302 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
304 for (std::size_t j=0; j<cols; ++j)
305 {
306 auto B_j = Impl::ColumnVectorView(matrixB, j);
307 auto result_j = Impl::ColumnVectorView(result, j);
308 matrixA.mv(B_j, result_j);
309 }
310 return result;
311 }
312
314 template<int l>
316 {
318
319 for (size_type i=0; i<l; i++) {
320 for (size_type j=0; j<cols; j++) {
321 C[i][j] = 0;
322 for (size_type k=0; k<rows; k++)
323 C[i][j] += M[i][k]*(*this)[k][j];
324 }
325 }
326 return C;
327 }
328
330
332 template <int r, int c>
334 {
335 static_assert(r == c, "Cannot rightmultiply with non-square matrix");
336 static_assert(r == cols, "Size mismatch");
338
339 for (size_type i=0; i<rows; i++)
340 for (size_type j=0; j<cols; j++) {
341 (*this)[i][j] = 0;
342 for (size_type k=0; k<cols; k++)
343 (*this)[i][j] += C[i][k]*M[k][j];
344 }
345 return *this;
346 }
347
349 template<int l>
351 {
353
354 for (size_type i=0; i<rows; i++) {
355 for (size_type j=0; j<l; j++) {
356 C[i][j] = 0;
357 for (size_type k=0; k<cols; k++)
358 C[i][j] += (*this)[i][k]*M[k][j];
359 }
360 }
361 return C;
362 }
363
364 // make this thing a matrix
365 static constexpr size_type mat_rows() { return ROWS; }
366 static constexpr size_type mat_cols() { return COLS; }
367
368 row_reference mat_access ( size_type i )
369 {
370 DUNE_ASSERT_BOUNDS(i < ROWS);
371 return _data[i];
372 }
373
374 const_row_reference mat_access ( size_type i ) const
375 {
376 DUNE_ASSERT_BOUNDS(i < ROWS);
377 return _data[i];
378 }
379 };
380
381#ifndef DOXYGEN // hide specialization
384 template<class K>
385 class FieldMatrix<K,1,1> : public DenseMatrix< FieldMatrix<K,1,1> >
386 {
387 FieldVector<K,1> _data;
388 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
389 public:
390 // standard constructor and everything is sufficient ...
391
392 //===== type definitions and constants
393
395 typedef typename Base::size_type size_type;
396
399 constexpr static int blocklevel = 1;
400
401 typedef typename Base::row_type row_type;
402
403 typedef typename Base::row_reference row_reference;
404 typedef typename Base::const_row_reference const_row_reference;
405
408 constexpr static int rows = 1;
411 constexpr static int cols = 1;
412
413 //===== constructors
416 constexpr FieldMatrix() = default;
417
420 FieldMatrix(std::initializer_list<Dune::FieldVector<K, 1>> const &l)
421 {
422 std::copy_n(l.begin(), std::min<std::size_t>(1, l.size()), &_data);
423 }
424
425 template <class T,
426 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
427 FieldMatrix(T const& rhs)
428 {
429 *this = rhs;
430 }
431
432 using Base::operator=;
433
435 FieldMatrix<K, 1, 1> transposed() const
436 {
437 return *this;
438 }
439
441 template <class OtherScalar>
442 friend auto operator+ ( const FieldMatrix& matrixA,
443 const FieldMatrix<OtherScalar,1,1>& matrixB)
444 {
445 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] + matrixB[0][0]};
446 }
447
449 template <class Scalar,
450 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
451 friend auto operator+ ( const FieldMatrix& matrix,
452 const Scalar& scalar)
453 {
454 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
455 }
456
458 template <class Scalar,
459 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
460 friend auto operator+ ( const Scalar& scalar,
461 const FieldMatrix& matrix)
462 {
463 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar + matrix[0][0]};
464 }
465
467 template <class OtherScalar>
468 friend auto operator- ( const FieldMatrix& matrixA,
469 const FieldMatrix<OtherScalar,1,1>& matrixB)
470 {
471 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] - matrixB[0][0]};
472 }
473
475 template <class Scalar,
476 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
477 friend auto operator- ( const FieldMatrix& matrix,
478 const Scalar& scalar)
479 {
480 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
481 }
482
484 template <class Scalar,
485 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
486 friend auto operator- ( const Scalar& scalar,
487 const FieldMatrix& matrix)
488 {
489 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar - matrix[0][0]};
490 }
491
493 template <class Scalar,
494 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
495 friend auto operator* ( const FieldMatrix& matrix, Scalar scalar)
496 {
497 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] * scalar};
498 }
499
501 template <class Scalar,
502 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
503 friend auto operator* ( Scalar scalar, const FieldMatrix& matrix)
504 {
505 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {scalar * matrix[0][0]};
506 }
507
509 template <class Scalar,
510 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
511 friend auto operator/ ( const FieldMatrix& matrix, Scalar scalar)
512 {
513 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] / scalar};
514 }
515
516 //===== solve
517
520 template <class OtherScalar, int otherCols>
521 friend auto operator* ( const FieldMatrix& matrixA,
522 const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
523 {
524 FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
525
526 for (size_type j = 0; j < matrixB.mat_cols(); ++j)
527 result[0][j] = matrixA[0][0] * matrixB[0][j];
528
529 return result;
530 }
531
538 template <class OtherMatrix, std::enable_if_t<
539 Impl::IsStaticSizeMatrix_v<OtherMatrix>
540 and not Impl::IsFieldMatrix_v<OtherMatrix>
541 and (OtherMatrix::rows==1)
542 , int> = 0>
543 friend auto operator* ( const FieldMatrix& matrixA,
544 const OtherMatrix& matrixB)
545 {
546 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
548 for (std::size_t j=0; j<rows; ++j)
549 matrixB.mtv(matrixA[j], result[j]);
550 return result;
551 }
552
559 template <class OtherMatrix, std::enable_if_t<
560 Impl::IsStaticSizeMatrix_v<OtherMatrix>
561 and not Impl::IsFieldMatrix_v<OtherMatrix>
562 and (OtherMatrix::cols==1)
563 , int> = 0>
564 friend auto operator* ( const OtherMatrix& matrixA,
565 const FieldMatrix& matrixB)
566 {
567 using Field = typename PromotionTraits<K, typename OtherMatrix::field_type>::PromotedType;
569 for (std::size_t j=0; j<cols; ++j)
570 {
571 auto B_j = Impl::ColumnVectorView(matrixB, j);
572 auto result_j = Impl::ColumnVectorView(result, j);
573 matrixA.mv(B_j, result_j);
574 }
575 return result;
576 }
577
579 template<int l>
580 FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
581 {
582 FieldMatrix<K,l,1> C;
583 for (size_type j=0; j<l; j++)
584 C[j][0] = M[j][0]*(*this)[0][0];
585 return C;
586 }
587
590 {
591 _data[0] *= M[0][0];
592 return *this;
593 }
594
596 template<int l>
597 FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
598 {
599 FieldMatrix<K,1,l> C;
600
601 for (size_type j=0; j<l; j++)
602 C[0][j] = M[0][j]*_data[0];
603 return C;
604 }
605
606 // make this thing a matrix
607 static constexpr size_type mat_rows() { return 1; }
608 static constexpr size_type mat_cols() { return 1; }
609
610 row_reference mat_access ([[maybe_unused]] size_type i)
611 {
612 DUNE_ASSERT_BOUNDS(i == 0);
613 return _data;
614 }
615
616 const_row_reference mat_access ([[maybe_unused]] size_type i) const
617 {
618 DUNE_ASSERT_BOUNDS(i == 0);
619 return _data;
620 }
621
623 FieldMatrix& operator+= (const K& k)
624 {
625 _data[0] += k;
626 return (*this);
627 }
628
630 FieldMatrix& operator-= (const K& k)
631 {
632 _data[0] -= k;
633 return (*this);
634 }
635
637 FieldMatrix& operator*= (const K& k)
638 {
639 _data[0] *= k;
640 return (*this);
641 }
642
644 FieldMatrix& operator/= (const K& k)
645 {
646 _data[0] /= k;
647 return (*this);
648 }
649
650 //===== conversion operator
651
652 operator const K& () const { return _data[0]; }
653
654 };
655
657 template<typename K>
658 std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
659 {
660 s << a[0][0];
661 return s;
662 }
663
664#endif // DOXYGEN
665
666 namespace FMatrixHelp {
667
669 template <typename K>
670 static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
671 {
672 using real_type = typename FieldTraits<K>::real_type;
673 inverse[0][0] = real_type(1.0)/matrix[0][0];
674 return matrix[0][0];
675 }
676
678 template <typename K>
679 static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
680 {
681 return invertMatrix(matrix,inverse);
682 }
683
684
686 template <typename K>
687 static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
688 {
689 using real_type = typename FieldTraits<K>::real_type;
690 // code generated by maple
691 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
692 K det_1 = real_type(1.0)/det;
693 inverse[0][0] = matrix[1][1] * det_1;
694 inverse[0][1] = - matrix[0][1] * det_1;
695 inverse[1][0] = - matrix[1][0] * det_1;
696 inverse[1][1] = matrix[0][0] * det_1;
697 return det;
698 }
699
702 template <typename K>
703 static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
704 {
705 using real_type = typename FieldTraits<K>::real_type;
706 // code generated by maple
707 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
708 K det_1 = real_type(1.0)/det;
709 inverse[0][0] = matrix[1][1] * det_1;
710 inverse[1][0] = - matrix[0][1] * det_1;
711 inverse[0][1] = - matrix[1][0] * det_1;
712 inverse[1][1] = matrix[0][0] * det_1;
713 return det;
714 }
715
717 template <typename K>
718 static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
719 {
720 using real_type = typename FieldTraits<K>::real_type;
721 // code generated by maple
722 K t4 = matrix[0][0] * matrix[1][1];
723 K t6 = matrix[0][0] * matrix[1][2];
724 K t8 = matrix[0][1] * matrix[1][0];
725 K t10 = matrix[0][2] * matrix[1][0];
726 K t12 = matrix[0][1] * matrix[2][0];
727 K t14 = matrix[0][2] * matrix[2][0];
728
729 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
730 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
731 K t17 = real_type(1.0)/det;
732
733 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
734 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
735 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
736 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
737 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
738 inverse[1][2] = -(t6-t10) * t17;
739 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
740 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
741 inverse[2][2] = (t4-t8) * t17;
742
743 return det;
744 }
745
747 template <typename K>
748 static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
749 {
750 using real_type = typename FieldTraits<K>::real_type;
751 // code generated by maple
752 K t4 = matrix[0][0] * matrix[1][1];
753 K t6 = matrix[0][0] * matrix[1][2];
754 K t8 = matrix[0][1] * matrix[1][0];
755 K t10 = matrix[0][2] * matrix[1][0];
756 K t12 = matrix[0][1] * matrix[2][0];
757 K t14 = matrix[0][2] * matrix[2][0];
758
759 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
760 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
761 K t17 = real_type(1.0)/det;
762
763 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
764 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
765 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
766 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
767 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
768 inverse[2][1] = -(t6-t10) * t17;
769 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
770 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
771 inverse[2][2] = (t4-t8) * t17;
772
773 return det;
774 }
775
777 template< class K, int m, int n, int p >
778 static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
779 const FieldMatrix< K, n, p > &B,
781 {
782 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
783
784 for( size_type i = 0; i < m; ++i )
785 {
786 for( size_type j = 0; j < p; ++j )
787 {
788 ret[ i ][ j ] = K( 0 );
789 for( size_type k = 0; k < n; ++k )
790 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
791 }
792 }
793 }
794
796 template <typename K, int rows, int cols>
798 {
799 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
800
801 for(size_type i=0; i<cols; i++)
802 for(size_type j=0; j<cols; j++)
803 {
804 ret[i][j]=0.0;
805 for(size_type k=0; k<rows; k++)
806 ret[i][j]+=matrix[k][i]*matrix[k][j];
807 }
808 }
809
810 using Dune::DenseMatrixHelp::multAssign;
811
813 template <typename K, int rows, int cols>
815 {
816 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
817
818 for(size_type i=0; i<cols; ++i)
819 {
820 ret[i] = 0.0;
821 for(size_type j=0; j<rows; ++j)
822 ret[i] += matrix[j][i]*x[j];
823 }
824 }
825
827 template <typename K, int rows, int cols>
829 {
831 multAssign(matrix,x,ret);
832 return ret;
833 }
834
836 template <typename K, int rows, int cols>
838 {
840 multAssignTransposed( matrix, x, ret );
841 return ret;
842 }
843
844 } // end namespace FMatrixHelp
845
848} // end namespace
849
850#include "fmatrixev.hh"
851#endif
Macro for wrapping boundary checks.
A dense n x m matrix.
Definition: densematrix.hh:140
derived_type operator-() const
Matrix negation.
Definition: densematrix.hh:298
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:387
constexpr size_type M() const
number of columns
Definition: densematrix.hh:703
FieldMatrix< K, ROWS, COLS > & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:645
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:329
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:321
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition: densematrix.hh:312
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: densematrix.hh:178
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:169
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:166
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:175
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:172
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition: densematrix.hh:289
A dense n x m matrix.
Definition: fmatrix.hh:117
constexpr FieldMatrix()=default
Default constructor.
FieldMatrix & operator=(const FieldMatrix< T, ROWS, COLS > &x)
copy assignment from FieldMatrix over a different field
Definition: fmatrix.hh:160
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:350
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:315
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:333
friend auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition: fmatrix.hh:211
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
FieldMatrix< K, COLS, ROWS > transposed() const
Return transposed of the matrix as FieldMatrix.
Definition: fmatrix.hh:171
static constexpr int cols
The number of columns.
Definition: fmatrix.hh:125
friend auto operator/(const FieldMatrix &matrix, Scalar scalar)
vector space division by scalar
Definition: fmatrix.hh:239
friend auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition: fmatrix.hh:182
FieldMatrix & operator=(const FieldMatrix &)=default
copy assignment operator
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
A few common exception classes.
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:837
static K invertMatrix_retTransposed(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:679
static 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:778
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:670
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:828
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:797
static 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:814
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  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)