DUNE-FEM (unstable)

diagonalmatrix.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_DIAGONAL_MATRIX_HH
6#define DUNE_DIAGONAL_MATRIX_HH
7
12#include <algorithm>
13#include <cassert>
14#include <cmath>
15#include <complex>
16#include <cstddef>
17#include <initializer_list>
18#include <iostream>
19#include <memory>
20
29#include <dune/common/matrixconcepts.hh>
31
32
33namespace Dune {
34
35 template< class K, int n > class DiagonalRowVectorConst;
36 template< class K, int n > class DiagonalRowVector;
37 template< class DiagonalMatrixType > class DiagonalMatrixWrapper;
38 template< class C, class T, class R> class ContainerWrapperIterator;
39
54 template<class K, int n>
56 {
57 template<class, int> friend class DiagonalMatrix;
58 typedef DiagonalMatrixWrapper< DiagonalMatrix<K,n> > WrapperType;
59
60 public:
61 //===== type definitions and constants
62
64 typedef K value_type;
65 typedef value_type field_type;
66
68 typedef K block_type;
69
71 typedef std::size_t size_type;
72
74 constexpr static int blocklevel = 1;
75
77 typedef DiagonalRowVector<K,n> row_type;
78 typedef row_type reference;
79 typedef row_type row_reference;
80 typedef DiagonalRowVectorConst<K,n> const_row_type;
81 typedef const_row_type const_reference;
82 typedef const_row_type const_row_reference;
83
85 constexpr static int rows = n;
87 constexpr static int cols = n;
88
89 //==== size
90
91 static constexpr size_type size ()
92 {
93 return rows;
94 }
95
96 //===== constructors
97
99 constexpr DiagonalMatrix() = default;
100
102 DiagonalMatrix (const K& k)
103 : diag_(k)
104 {}
105
108 : diag_(diag)
109 {}
110
119 DiagonalMatrix (std::initializer_list<K> const &l)
120 {
121 std::copy_n(l.begin(), std::min<std::size_t>(rows, l.size()),
122 diag_.begin());
123 }
124
126 template <class OtherK>
128 : diag_(other.diag_)
129 {}
130
133 {
134 diag_ = k;
135 return *this;
136 }
137
139 template <class OtherK>
141 {
142 diag_ = other.diag_;
143 return *this;
144 }
145
147 bool identical(const DiagonalMatrix<K,n>& other) const
148 {
149 return (this==&other);
150 }
151
154 {
155 return *this;
156 }
157
158 //===== iterator interface to rows of the matrix
167
170 {
171 return Iterator(WrapperType(this),0);
172 }
173
176 {
177 return Iterator(WrapperType(this),n);
178 }
179
183 {
184 return Iterator(WrapperType(this),n-1);
185 }
186
190 {
191 return Iterator(WrapperType(this),-1);
192 }
193
194
203
206 {
207 return ConstIterator(WrapperType(this),0);
208 }
209
212 {
213 return ConstIterator(WrapperType(this),n);
214 }
215
219 {
220 return ConstIterator(WrapperType(this),n-1);
221 }
222
226 {
227 return ConstIterator(WrapperType(this),-1);
228 }
229
230
231
232 //===== vector space arithmetic
233
236 {
237 diag_ += y.diag_;
238 return *this;
239 }
240
243 {
244 diag_ -= y.diag_;
245 return *this;
246 }
247
250 {
251 diag_ += k;
252 return *this;
253 }
254
257 {
258 diag_ -= k;
259 return *this;
260 }
261
264 {
265 diag_ *= k;
266 return *this;
267 }
268
271 {
272 diag_ /= k;
273 return *this;
274 }
275
276 //===== comparison ops
277
279 bool operator==(const DiagonalMatrix& other) const
280 {
281 return diag_==other.diagonal();
282 }
283
285 bool operator!=(const DiagonalMatrix& other) const
286 {
287 return diag_!=other.diagonal();
288 }
289
290
291 //===== linear maps
292
294 template<class X, class Y>
295 void mv (const X& x, Y& y) const
296 {
297#ifdef DUNE_FMatrix_WITH_CHECKING
298 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
299 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
300#endif
301 for (size_type i=0; i<n; ++i)
302 y[i] = diag_[i] * x[i];
303 }
304
306 template<class X, class Y>
307 void mtv (const X& x, Y& y) const
308 {
309 mv(x, y);
310 }
311
313 template<class X, class Y>
314 void umv (const X& x, Y& y) const
315 {
316#ifdef DUNE_FMatrix_WITH_CHECKING
317 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
318 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
319#endif
320 for (size_type i=0; i<n; ++i)
321 y[i] += diag_[i] * x[i];
322 }
323
325 template<class X, class Y>
326 void umtv (const X& x, Y& y) const
327 {
328#ifdef DUNE_FMatrix_WITH_CHECKING
329 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
330 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
331#endif
332 for (size_type i=0; i<n; ++i)
333 y[i] += diag_[i] * x[i];
334 }
335
337 template<class X, class Y>
338 void umhv (const X& x, Y& y) const
339 {
340#ifdef DUNE_FMatrix_WITH_CHECKING
341 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
342 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
343#endif
344 for (size_type i=0; i<n; i++)
345 y[i] += conjugateComplex(diag_[i])*x[i];
346 }
347
349 template<class X, class Y>
350 void mmv (const X& x, Y& y) const
351 {
352#ifdef DUNE_FMatrix_WITH_CHECKING
353 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
354 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
355#endif
356 for (size_type i=0; i<n; ++i)
357 y[i] -= diag_[i] * x[i];
358 }
359
361 template<class X, class Y>
362 void mmtv (const X& x, Y& y) const
363 {
364#ifdef DUNE_FMatrix_WITH_CHECKING
365 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
366 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
367#endif
368 for (size_type i=0; i<n; ++i)
369 y[i] -= diag_[i] * x[i];
370 }
371
373 template<class X, class Y>
374 void mmhv (const X& x, Y& y) const
375 {
376#ifdef DUNE_FMatrix_WITH_CHECKING
377 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
378 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
379#endif
380 for (size_type i=0; i<n; i++)
381 y[i] -= conjugateComplex(diag_[i])*x[i];
382 }
383
385 template<class X, class Y>
386 void usmv (const typename FieldTraits<Y>::field_type & alpha,
387 const X& x, Y& y) const
388 {
389#ifdef DUNE_FMatrix_WITH_CHECKING
390 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
391 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
392#endif
393 for (size_type i=0; i<n; i++)
394 y[i] += alpha * diag_[i] * x[i];
395 }
396
398 template<class X, class Y>
399 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
400 const X& x, Y& y) const
401 {
402#ifdef DUNE_FMatrix_WITH_CHECKING
403 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
404 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
405#endif
406 for (size_type i=0; i<n; i++)
407 y[i] += alpha * diag_[i] * x[i];
408 }
409
411 template<class X, class Y>
412 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
413 const X& x, Y& y) const
414 {
415#ifdef DUNE_FMatrix_WITH_CHECKING
416 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
417 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
418#endif
419 for (size_type i=0; i<n; i++)
420 y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
421 }
422
423 //===== norms
424
426 double frobenius_norm () const
427 {
428 return diag_.two_norm();
429 }
430
432 double frobenius_norm2 () const
433 {
434 return diag_.two_norm2();
435 }
436
438 double infinity_norm () const
439 {
440 return diag_.infinity_norm();
441 }
442
444 double infinity_norm_real () const
445 {
446 return diag_.infinity_norm_real();
447 }
448
449
450
451 //===== solve
452
454 template<class V>
455 void solve (V& x, const V& b) const
456 {
457 for (int i=0; i<n; i++)
458 x[i] = b[i]/diag_[i];
459 }
460
462 void invert()
463 {
464 using real_type = typename FieldTraits<K>::real_type;
465 for (int i=0; i<n; i++)
466 diag_[i] = real_type(1.0)/diag_[i];
467 }
468
470 K determinant () const
471 {
472 K det = diag_[0];
473 for (int i=1; i<n; i++)
474 det *= diag_[i];
475 return det;
476 }
477
478
479
480 //===== matrix-matrix multiplication
481
484 template <class OtherScalar>
485 friend auto operator* ( const DiagonalMatrix& matrixA,
486 const DiagonalMatrix<OtherScalar, n>& matrixB)
487 {
489 for(int i=0; i<n; ++i)
490 result.diagonal(i) = matrixA.diagonal(i)*matrixB.diagonal(i);
491 return result;
492 }
493
503 template <class OtherMatrix,
504 std::enable_if_t<(Impl::IsDenseMatrix<OtherMatrix>::value), int> = 0,
505 std::enable_if_t<(not Impl::IsFieldMatrix<OtherMatrix>::value), int> = 0>
506 friend auto operator* ( const DiagonalMatrix& matrixA,
507 const OtherMatrix& matrixB)
508 {
509 using OtherField = typename FieldTraits<OtherMatrix>::field_type;
510 using F = typename PromotionTraits<field_type, OtherField>::PromotedType;
511
512 auto result = [&]{
513 if constexpr (Impl::IsStaticSizeMatrix_v<OtherMatrix>) {
514 static_assert(n == OtherMatrix::rows);
516 } else {
517 assert(n == matrixB.N());
518 return DynamicMatrix<F>{n,matrixB.M()};
519 }
520 }();
521
522 for (int i = 0; i < result.N(); ++i)
523 for (int j = 0; j < result.M(); ++j)
524 result[i][j] = matrixA.diagonal(i) * matrixB[i][j];
525 return result;
526 }
527
528 //===== sizes
529
531 static constexpr size_type N ()
532 {
533 return n;
534 }
535
537 static constexpr size_type M ()
538 {
539 return n;
540 }
541
542
543
544 //===== query
545
547 bool exists (size_type i, size_type j) const
548 {
549 DUNE_ASSERT_BOUNDS(i >= 0 && i < n);
550 DUNE_ASSERT_BOUNDS(j >= 0 && j < n);
551 return i==j;
552 }
553
554
555
557 friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
558 {
559 for (size_type i=0; i<n; i++) {
560 for (size_type j=0; j<n; j++)
561 s << ((i==j) ? a.diag_[i] : 0) << " ";
562 s << std::endl;
563 }
564 return s;
565 }
566
569 {
570 return reference(const_cast<K*>(&diag_[i]), i);
571 }
572
574 const_reference operator[](size_type i) const
575 {
576 return const_reference(const_cast<K*>(&diag_[i]), i);
577 }
578
580 const K& diagonal(size_type i) const
581 {
582 return diag_[i];
583 }
584
587 {
588 return diag_[i];
589 }
590
593 {
594 return diag_;
595 }
596
599 {
600 return diag_;
601 }
602
603 private:
604
605 // the data, a FieldVector storing the diagonal
606 FieldVector<K,n> diag_;
607 };
608
609 template< class K, int n >
610 struct FieldTraits< DiagonalMatrix<K,n> >
611 {
612 typedef typename FieldTraits<K>::field_type field_type;
613 typedef typename FieldTraits<K>::real_type real_type;
614 };
615
616
617#ifndef DOXYGEN // hide specialization
620 template< class K >
621 class DiagonalMatrix<K, 1> : public FieldMatrix<K, 1, 1>
622 {
623 typedef FieldMatrix<K,1,1> Base;
624 public:
626 typedef typename Base::size_type size_type;
627
630 constexpr static int blocklevel = 1;
631
632 typedef typename Base::row_type row_type;
633
634 typedef typename Base::row_reference row_reference;
635 typedef typename Base::const_row_reference const_row_reference;
636
639 constexpr static int rows = 1;
642 constexpr static int cols = 1;
643
644
646 constexpr DiagonalMatrix() = default;
647
649 DiagonalMatrix(const K& scalar)
650 {
651 (*this)[0][0] = scalar;
652 }
653
655 template <class OtherK>
656 DiagonalMatrix(const DiagonalMatrix<OtherK,1>& other)
657 : Base(FieldMatrix<OtherK,1,1>(other))
658 {}
659
661 const K& diagonal(size_type) const
662 {
663 return (*this)[0][0];
664 }
665
668 {
669 return (*this)[0][0];
670 }
671
673 const FieldVector<K,1>& diagonal() const
674 {
675 return (*this)[0];
676 }
677
679 FieldVector<K,1>& diagonal()
680 {
681 return (*this)[0];
682 }
683
685 DiagonalMatrix<K, 1> transposed() const
686 {
687 return *this;
688 }
689
692 template <class OtherScalar>
693 friend auto operator* ( const DiagonalMatrix& matrixA,
694 const DiagonalMatrix<OtherScalar, 1>& matrixB)
695 {
696 return DiagonalMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType, 1>{matrixA.diagonal(0)*matrixB.diagonal(0)};
697 }
698
699 template <class OtherMatrix,
700 std::enable_if_t<(Impl::IsDenseMatrix<OtherMatrix>::value), int> = 0,
701 std::enable_if_t<(not Impl::IsFieldMatrix<OtherMatrix>::value), int> = 0>
702 friend auto operator* ( const DiagonalMatrix& matrixA,
703 const OtherMatrix& matrixB)
704 {
705 using OtherField = typename FieldTraits<OtherMatrix>::field_type;
706 using F = typename PromotionTraits<K, OtherField>::PromotedType;
707
708 auto result = [&]{
709 if constexpr (Impl::IsStaticSizeMatrix_v<OtherMatrix>) {
710 static_assert(1 == OtherMatrix::rows);
711 return FieldMatrix<F, 1, OtherMatrix::cols>{};
712 } else {
713 assert(1 == matrixB.N());
714 return DynamicMatrix<F>{1,matrixB.M()};
715 }
716 }();
717
718 for (int i = 0; i < result.N(); ++i)
719 for (int j = 0; j < result.M(); ++j)
720 result[i][j] = matrixA.diagonal(i) * matrixB[i][j];
721 return result;
722 }
723
724 };
725#endif
726
727
728 template<class DiagonalMatrixType>
729 class DiagonalMatrixWrapper
730 {
731 typedef typename DiagonalMatrixType::reference reference;
732 typedef typename DiagonalMatrixType::const_reference const_reference;
733 typedef typename DiagonalMatrixType::field_type K;
734 typedef DiagonalRowVector<K, DiagonalMatrixType::rows> row_type;
735 typedef std::size_t size_type;
736 typedef DiagonalMatrixWrapper< DiagonalMatrixType> MyType;
737
738 friend class ContainerWrapperIterator<const MyType, reference, reference>;
739 friend class ContainerWrapperIterator<const MyType, const_reference, const_reference>;
740
741 public:
742
743 DiagonalMatrixWrapper() :
744 mat_(0)
745 {}
746
747 DiagonalMatrixWrapper(const DiagonalMatrixType* mat) :
748 mat_(const_cast<DiagonalMatrixType*>(mat))
749 {}
750
751 size_type realIndex(int i) const
752 {
753 return i;
754 }
755
756 row_type* pointer(int i) const
757 {
758 row_ = row_type(&(mat_->diagonal(i)), i);
759 return &row_;
760 }
761
762 bool identical(const DiagonalMatrixWrapper& other) const
763 {
764 return mat_==other.mat_;
765 }
766
767 private:
768
769 mutable DiagonalMatrixType* mat_;
770 mutable row_type row_;
771 };
772
776 template< class K, int n >
777 class DiagonalRowVectorConst
778 {
779 template<class DiagonalMatrixType>
780 friend class DiagonalMatrixWrapper;
781 friend class ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&>;
782
783 public:
784 // remember size of vector
785 constexpr static int dimension = n;
786
787 // standard constructor and everything is sufficient ...
788
789 //===== type definitions and constants
790
792 typedef K field_type;
793
795 typedef K block_type;
796
798 typedef std::size_t size_type;
799
801 constexpr static int blocklevel = 1;
802
804 constexpr static int size = n;
805
808 p_(0),
809 row_(0)
810 {}
811
813 explicit DiagonalRowVectorConst (K* p, int col) :
814 p_(p),
815 row_(col)
816 {}
817
818 //===== access to components
819
821 const K& operator[] ([[maybe_unused]] size_type i) const
822 {
823 DUNE_ASSERT_BOUNDS(i == row_);
824 return *p_;
825 }
826
827 // check if row is identical to other row (not only identical values)
828 // since this is a proxy class we need to check equality of the stored pointer
829 bool identical(const DiagonalRowVectorConst<K,n>& other) const
830 {
831 return ((p_ == other.p_)and (row_ == other.row_));
832 }
833
838
841 {
842 return ConstIterator(*this,0);
843 }
844
847 {
848 return ConstIterator(*this,1);
849 }
850
854 {
855 return ConstIterator(*this,0);
856 }
857
861 {
862 return ConstIterator(*this,-1);
863 }
864
866 bool operator== (const DiagonalRowVectorConst& y) const
867 {
868 return ((p_==y.p_)and (row_==y.row_));
869 }
870
871 //===== sizes
872
874 size_type N () const
875 {
876 return n;
877 }
878
880 size_type dim () const
881 {
882 return n;
883 }
884
887 {
888 return row_;
889 }
890
892 const K& diagonal() const
893 {
894 return *p_;
895 }
896
897 protected:
898
899 size_type realIndex([[maybe_unused]] int i) const
900 {
901 return rowIndex();
902 }
903
904 K* pointer([[maybe_unused]] size_type i) const
905 {
906 return const_cast<K*>(p_);
907 }
908
909 DiagonalRowVectorConst* operator&()
910 {
911 return this;
912 }
913
914 // the data, very simply a pointer to the diagonal value and the row number
915 K* p_;
916 size_type row_;
917 };
918
919 template< class K, int n >
920 class DiagonalRowVector : public DiagonalRowVectorConst<K,n>
921 {
922 template<class DiagonalMatrixType>
923 friend class DiagonalMatrixWrapper;
924 friend class ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&>;
925
926 public:
927 // standard constructor and everything is sufficient ...
928
929 //===== type definitions and constants
930
932 typedef K field_type;
933
935 typedef K block_type;
936
938 typedef std::size_t size_type;
939
941 DiagonalRowVector() : DiagonalRowVectorConst<K,n>()
942 {}
943
945 explicit DiagonalRowVector (K* p, int col) : DiagonalRowVectorConst<K,n>(p, col)
946 {}
947
948 //===== assignment from scalar
950 DiagonalRowVector& operator= (const K& k)
951 {
952 *p_ = k;
953 return *this;
954 }
955
956 //===== access to components
957
959 K& operator[] ([[maybe_unused]] size_type i)
960 {
961 DUNE_ASSERT_BOUNDS(i == row_);
962 return *p_;
963 }
964
969
972 {
973 return Iterator(*this, 0);
974 }
975
978 {
979 return Iterator(*this, 1);
980 }
981
985 {
986 return Iterator(*this, 0);
987 }
988
992 {
993 return Iterator(*this, -1);
994 }
995
1000
1001 using DiagonalRowVectorConst<K,n>::identical;
1002 using DiagonalRowVectorConst<K,n>::operator[];
1003 using DiagonalRowVectorConst<K,n>::operator==;
1004 using DiagonalRowVectorConst<K,n>::begin;
1005 using DiagonalRowVectorConst<K,n>::end;
1006 using DiagonalRowVectorConst<K,n>::beforeEnd;
1007 using DiagonalRowVectorConst<K,n>::beforeBegin;
1008 using DiagonalRowVectorConst<K,n>::N;
1009 using DiagonalRowVectorConst<K,n>::dim;
1010 using DiagonalRowVectorConst<K,n>::rowIndex;
1011 using DiagonalRowVectorConst<K,n>::diagonal;
1012
1013 protected:
1014
1015 DiagonalRowVector* operator&()
1016 {
1017 return this;
1018 }
1019
1020 private:
1021
1022 using DiagonalRowVectorConst<K,n>::p_;
1023 using DiagonalRowVectorConst<K,n>::row_;
1024 };
1025
1026
1027 // implement type traits
1028 template<class K, int n>
1029 struct const_reference< DiagonalRowVector<K,n> >
1030 {
1031 typedef DiagonalRowVectorConst<K,n> type;
1032 };
1033
1034 template<class K, int n>
1035 struct const_reference< DiagonalRowVectorConst<K,n> >
1036 {
1037 typedef DiagonalRowVectorConst<K,n> type;
1038 };
1039
1040 template<class K, int n>
1041 struct mutable_reference< DiagonalRowVector<K,n> >
1042 {
1043 typedef DiagonalRowVector<K,n> type;
1044 };
1045
1046 template<class K, int n>
1047 struct mutable_reference< DiagonalRowVectorConst<K,n> >
1048 {
1049 typedef DiagonalRowVector<K,n> type;
1050 };
1051
1052
1053
1076 template<class CW, class T, class R>
1077 class ContainerWrapperIterator : public BidirectionalIteratorFacade<ContainerWrapperIterator<CW,T,R>,T, R, int>
1078 {
1079 typedef typename std::remove_const<CW>::type NonConstCW;
1080
1081 friend class ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type>;
1082 friend class ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type>;
1083
1084 typedef ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type> MyType;
1085 typedef ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type> MyConstType;
1086
1087 public:
1088
1089 // Constructors needed by the facade iterators.
1091 containerWrapper_(),
1092 position_(0)
1093 {}
1094
1095 ContainerWrapperIterator(CW containerWrapper, int position) :
1096 containerWrapper_(containerWrapper),
1097 position_(position)
1098 {}
1099
1100 template<class OtherContainerWrapperIteratorType>
1101 ContainerWrapperIterator(OtherContainerWrapperIteratorType& other) :
1102 containerWrapper_(other.containerWrapper_),
1103 position_(other.position_)
1104 {}
1105
1106 ContainerWrapperIterator(const MyType& other) :
1107 containerWrapper_(other.containerWrapper_),
1108 position_(other.position_)
1109 {}
1110
1112 containerWrapper_(other.containerWrapper_),
1113 position_(other.position_)
1114 {}
1115
1116 template<class OtherContainerWrapperIteratorType>
1117 ContainerWrapperIterator& operator=(OtherContainerWrapperIteratorType& other)
1118 {
1119 containerWrapper_ = other.containerWrapper_;
1120 position_ = other.position_;
1121 return *this;
1122 }
1123
1124 // This operator is needed since we can not get the address of the
1125 // temporary object returned by dereference
1126 T* operator->() const
1127 {
1128 return containerWrapper_.pointer(position_);
1129 }
1130
1131 // Methods needed by the forward iterator
1132 bool equals(const MyType& other) const
1133 {
1134 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1135 }
1136
1137 bool equals(const MyConstType& other) const
1138 {
1139 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1140 }
1141
1142 R dereference() const
1143 {
1144 return *containerWrapper_.pointer(position_);
1145 }
1146
1147 void increment()
1148 {
1149 ++position_;
1150 }
1151
1152 // Additional function needed by BidirectionalIterator
1153 void decrement()
1154 {
1155 --position_;
1156 }
1157
1158 // Additional function needed by RandomAccessIterator
1159 R elementAt(int i) const
1160 {
1161 return *containerWrapper_.pointer(position_+i);
1162 }
1163
1164 void advance(int n)
1165 {
1166 position_=position_+n;
1167 }
1168
1169 template<class OtherContainerWrapperIteratorType>
1170 std::ptrdiff_t distanceTo(OtherContainerWrapperIteratorType& other) const
1171 {
1172 assert(containerWrapper_.identical(other));
1173 return other.position_ - position_;
1174 }
1175
1176 std::ptrdiff_t index() const
1177 {
1178 return containerWrapper_.realIndex(position_);
1179 }
1180
1181 private:
1182 NonConstCW containerWrapper_;
1183 size_t position_;
1184 };
1185
1186 template <class DenseMatrix, class field, int N>
1188 static void apply(DenseMatrix& denseMatrix,
1189 DiagonalMatrix<field, N> const& rhs) {
1190#ifdef DUNE_CHECK_BOUNDS
1191 if (denseMatrix.M() != N || denseMatrix.N() != N)
1192 DUNE_THROW(Dune::RangeError, "Incompatible matrix dimensions in assignment.");
1193#endif // DUNE_CHECK_BOUNDS
1194
1195 denseMatrix = field(0);
1196 for (int i = 0; i < N; ++i)
1197 denseMatrix[i][i] = rhs.diagonal()[i];
1198 }
1199 };
1200 /* @} */
1201} // end namespace
1202#endif
Macro for wrapping boundary checks.
Facade class for stl conformant bidirectional iterators.
Definition: iteratorfacades.hh:275
Iterator class for sparse vector-like containers.
Definition: diagonalmatrix.hh:1078
A dense n x m matrix.
Definition: densematrix.hh:145
constexpr size_type M() const
number of columns
Definition: densematrix.hh:708
constexpr size_type N() const
number of rows
Definition: densematrix.hh:702
constexpr FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:642
constexpr FieldTraits< value_type >::real_type two_norm2() const
square of two norm (sum over squared values of entries), need for block recursion
Definition: densevector.hh:651
constexpr FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: densevector.hh:662
constexpr FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densevector.hh:678
constexpr Iterator begin()
begin iterator
Definition: densevector.hh:348
A diagonal matrix of static size.
Definition: diagonalmatrix.hh:56
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:131
A dense n x m matrix.
Definition: fmatrix.hh:117
Default exception class for range errors.
Definition: exceptions.hh:348
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
A few common exception classes.
This file implements a dense matrix with dynamic numbers of rows and columns.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Type traits to determine the type of reals (when working with complex numbers)
Implements a vector constructed from a given type representing a field and a compile-time given size.
Implements a generic iterator class for writing stl conformant iterators.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:30
ConstIterator beforeBegin() const
Definition: diagonalmatrix.hh:225
DiagonalMatrix< K, n > transposed() const
Return transposed of the matrix as DiagonalMatrix.
Definition: diagonalmatrix.hh:153
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: diagonalmatrix.hh:374
DiagonalMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: diagonalmatrix.hh:263
std::size_t size_type
The type used for the index access and size operations.
Definition: diagonalmatrix.hh:71
size_type dim() const
dimension of the vector space
Definition: diagonalmatrix.hh:880
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:200
static constexpr int rows
The number of rows.
Definition: diagonalmatrix.hh:85
static constexpr size_type M()
number of blocks in column direction
Definition: diagonalmatrix.hh:537
FieldVector< K, n > & diagonal()
Get reference to diagonal vector.
Definition: diagonalmatrix.hh:598
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: diagonalmatrix.hh:412
Iterator iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:968
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: diagonalmatrix.hh:74
K field_type
export the type representing the field
Definition: diagonalmatrix.hh:792
ConstIterator beforeEnd() const
Definition: diagonalmatrix.hh:853
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:202
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: diagonalmatrix.hh:547
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: diagonalmatrix.hh:196
DiagonalRowVector(K *p, int col)
Constructor making vector with identical coordinates.
Definition: diagonalmatrix.hh:945
K & diagonal(size_type i)
Get reference to diagonal entry.
Definition: diagonalmatrix.hh:586
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: diagonalmatrix.hh:455
static constexpr size_type N()
number of blocks in row direction
Definition: diagonalmatrix.hh:531
Iterator beforeBegin()
Definition: diagonalmatrix.hh:991
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: diagonalmatrix.hh:574
Iterator iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:162
ConstIterator begin() const
begin ConstIterator
Definition: diagonalmatrix.hh:840
ConstIterator const_iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:999
DiagonalMatrix & operator-=(const DiagonalMatrix &y)
vector space subtraction
Definition: diagonalmatrix.hh:242
DiagonalRowVectorConst(K *p, int col)
Constructor making vector with identical coordinates.
Definition: diagonalmatrix.hh:813
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: diagonalmatrix.hh:362
DiagonalMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: diagonalmatrix.hh:102
ContainerWrapperIterator< DiagonalRowVector< K, n >, K, K & > Iterator
Iterator class for sequential access.
Definition: diagonalmatrix.hh:966
Iterator beforeEnd()
Definition: diagonalmatrix.hh:984
void umtv(const X &x, Y &y) const
y += A^T x
Definition: diagonalmatrix.hh:326
ConstIterator const_iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:837
double infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: diagonalmatrix.hh:444
void umv(const X &x, Y &y) const
y += A x
Definition: diagonalmatrix.hh:314
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: diagonalmatrix.hh:160
void mv(const X &x, Y &y) const
y = A x
Definition: diagonalmatrix.hh:295
double frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: diagonalmatrix.hh:426
ConstIterator end() const
end iterator
Definition: diagonalmatrix.hh:211
static constexpr int cols
The number of columns.
Definition: diagonalmatrix.hh:87
size_type rowIndex() const
index of this row in surrounding matrix
Definition: diagonalmatrix.hh:886
bool operator!=(const DiagonalMatrix &other) const
incomparison operator
Definition: diagonalmatrix.hh:285
ConstIterator begin() const
begin iterator
Definition: diagonalmatrix.hh:205
Iterator begin()
begin iterator
Definition: diagonalmatrix.hh:971
void mmv(const X &x, Y &y) const
y -= A x
Definition: diagonalmatrix.hh:350
DiagonalMatrix & operator/=(const K &k)
vector space division by scalar
Definition: diagonalmatrix.hh:270
ConstIterator const_iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:198
bool identical(const DiagonalMatrix< K, n > &other) const
Check if matrix is the same object as the other matrix.
Definition: diagonalmatrix.hh:147
friend auto operator*(const DiagonalMatrix &matrixA, const DiagonalMatrix< OtherScalar, n > &matrixB)
Matrix-matrix multiplication.
Definition: diagonalmatrix.hh:485
Iterator end()
end iterator
Definition: diagonalmatrix.hh:175
DiagonalRowVector()
Constructor making uninitialized vector.
Definition: diagonalmatrix.hh:941
double frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: diagonalmatrix.hh:432
constexpr DiagonalMatrix()=default
Default constructor.
void mtv(const X &x, Y &y) const
y = A^T x
Definition: diagonalmatrix.hh:307
ContainerWrapperIterator< DiagonalRowVectorConst< K, n >, const K, const K & > ConstIterator
ConstIterator class for sequential access.
Definition: diagonalmatrix.hh:997
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: diagonalmatrix.hh:399
const FieldVector< K, n > & diagonal() const
Get const reference to diagonal vector.
Definition: diagonalmatrix.hh:592
void invert()
Compute inverse.
Definition: diagonalmatrix.hh:462
const K & diagonal(size_type i) const
Get const reference to diagonal entry.
Definition: diagonalmatrix.hh:580
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: diagonalmatrix.hh:874
DiagonalMatrix(std::initializer_list< K > const &l)
Construct diagonal matrix from an initializer list.
Definition: diagonalmatrix.hh:119
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:166
K field_type
export the type representing the field
Definition: diagonalmatrix.hh:932
reference operator[](size_type i)
Return reference object as row replacement.
Definition: diagonalmatrix.hh:568
Iterator end()
end iterator
Definition: diagonalmatrix.hh:977
DiagonalMatrix(const FieldVector< K, n > &diag)
Constructor initializing the diagonal with a vector.
Definition: diagonalmatrix.hh:107
std::size_t size_type
The type used for the index access and size operation.
Definition: diagonalmatrix.hh:938
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: diagonalmatrix.hh:386
double infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: diagonalmatrix.hh:438
DiagonalMatrix & operator+=(const DiagonalMatrix &y)
vector space addition
Definition: diagonalmatrix.hh:235
K block_type
export the type representing the components
Definition: diagonalmatrix.hh:68
void umhv(const X &x, Y &y) const
y += A^H x
Definition: diagonalmatrix.hh:338
Iterator beforeBegin()
Definition: diagonalmatrix.hh:189
ConstIterator end() const
end ConstIterator
Definition: diagonalmatrix.hh:846
std::size_t size_type
The type used for the index access and size operation.
Definition: diagonalmatrix.hh:798
DiagonalMatrix & operator=(const K &k)
Assignment from a scalar.
Definition: diagonalmatrix.hh:132
DiagonalRowVectorConst()
Constructor making uninitialized vector.
Definition: diagonalmatrix.hh:807
Iterator beforeEnd()
Definition: diagonalmatrix.hh:182
friend std::ostream & operator<<(std::ostream &s, const DiagonalMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: diagonalmatrix.hh:557
K value_type
export the type representing the field
Definition: diagonalmatrix.hh:64
K block_type
export the type representing the components
Definition: diagonalmatrix.hh:795
K determinant() const
calculates the determinant of this matrix
Definition: diagonalmatrix.hh:470
ConstIterator beforeBegin() const
Definition: diagonalmatrix.hh:860
ContainerWrapperIterator< DiagonalRowVectorConst< K, n >, const K, const K & > ConstIterator
ConstIterator class for sequential access.
Definition: diagonalmatrix.hh:835
Iterator begin()
begin iterator
Definition: diagonalmatrix.hh:169
bool operator==(const DiagonalMatrix &other) const
comparison operator
Definition: diagonalmatrix.hh:279
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: diagonalmatrix.hh:77
Iterator RowIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:164
K block_type
export the type representing the components
Definition: diagonalmatrix.hh:935
ConstIterator beforeEnd() const
Definition: diagonalmatrix.hh:218
DiagonalMatrix(const DiagonalMatrix< OtherK, n > &other)
Converting constructor.
Definition: diagonalmatrix.hh:127
const K & diagonal() const
the diagonal value
Definition: diagonalmatrix.hh:892
#define DUNE_THROW(E,...)
Definition: exceptions.hh:314
constexpr EnableIfInterOperable< T1, T2, bool >::type operator==(const ForwardIteratorFacade< T1, V1, R1, D > &lhs, const ForwardIteratorFacade< T2, V2, R2, D > &rhs)
Checks for equality.
Definition: iteratorfacades.hh:238
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:164
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition: densematrix.hh:61
get the 'mutable' version of a reference to a const object
Definition: genericiterator.hh:116
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 4, 22:38, 2025)