DUNE PDELab (git)

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
28
29
30namespace Dune {
31
32 template< class K, int n > class DiagonalRowVectorConst;
33 template< class K, int n > class DiagonalRowVector;
34 template< class DiagonalMatrixType > class DiagonalMatrixWrapper;
35 template< class C, class T, class R> class ContainerWrapperIterator;
36
51 template<class K, int n>
53 {
54 typedef DiagonalMatrixWrapper< DiagonalMatrix<K,n> > WrapperType;
55
56 public:
57 //===== type definitions and constants
58
60 typedef K value_type;
61 typedef value_type field_type;
62
64 typedef K block_type;
65
67 typedef std::size_t size_type;
68
70 constexpr static int blocklevel = 1;
71
73 typedef DiagonalRowVector<K,n> row_type;
74 typedef row_type reference;
75 typedef row_type row_reference;
76 typedef DiagonalRowVectorConst<K,n> const_row_type;
77 typedef const_row_type const_reference;
78 typedef const_row_type const_row_reference;
79
81 constexpr static int rows = n;
83 constexpr static int cols = n;
84
85 //==== size
86
87 static constexpr size_type size ()
88 {
89 return rows;
90 }
91
92 //===== constructors
93
95 constexpr DiagonalMatrix() = default;
96
98 DiagonalMatrix (const K& k)
99 : diag_(k)
100 {}
101
104 : diag_(diag)
105 {}
106
115 DiagonalMatrix (std::initializer_list<K> const &l)
116 {
117 std::copy_n(l.begin(), std::min<std::size_t>(rows, l.size()),
118 diag_.begin());
119 }
120
123 {
124 diag_ = k;
125 return *this;
126 }
127
129 bool identical(const DiagonalMatrix<K,n>& other) const
130 {
131 return (this==&other);
132 }
133
136 {
137 return *this;
138 }
139
140 //===== iterator interface to rows of the matrix
149
152 {
153 return Iterator(WrapperType(this),0);
154 }
155
158 {
159 return Iterator(WrapperType(this),n);
160 }
161
165 {
166 return Iterator(WrapperType(this),n-1);
167 }
168
172 {
173 return Iterator(WrapperType(this),-1);
174 }
175
176
185
188 {
189 return ConstIterator(WrapperType(this),0);
190 }
191
194 {
195 return ConstIterator(WrapperType(this),n);
196 }
197
201 {
202 return ConstIterator(WrapperType(this),n-1);
203 }
204
208 {
209 return ConstIterator(WrapperType(this),-1);
210 }
211
212
213
214 //===== vector space arithmetic
215
218 {
219 diag_ += y.diag_;
220 return *this;
221 }
222
225 {
226 diag_ -= y.diag_;
227 return *this;
228 }
229
232 {
233 diag_ += k;
234 return *this;
235 }
236
239 {
240 diag_ -= k;
241 return *this;
242 }
243
246 {
247 diag_ *= k;
248 return *this;
249 }
250
253 {
254 diag_ /= k;
255 return *this;
256 }
257
258 //===== comparison ops
259
261 bool operator==(const DiagonalMatrix& other) const
262 {
263 return diag_==other.diagonal();
264 }
265
267 bool operator!=(const DiagonalMatrix& other) const
268 {
269 return diag_!=other.diagonal();
270 }
271
272
273 //===== linear maps
274
276 template<class X, class Y>
277 void mv (const X& x, Y& y) const
278 {
279#ifdef DUNE_FMatrix_WITH_CHECKING
280 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
281 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
282#endif
283 for (size_type i=0; i<n; ++i)
284 y[i] = diag_[i] * x[i];
285 }
286
288 template<class X, class Y>
289 void mtv (const X& x, Y& y) const
290 {
291 mv(x, y);
292 }
293
295 template<class X, class Y>
296 void umv (const X& x, Y& y) const
297 {
298#ifdef DUNE_FMatrix_WITH_CHECKING
299 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
300 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
301#endif
302 for (size_type i=0; i<n; ++i)
303 y[i] += diag_[i] * x[i];
304 }
305
307 template<class X, class Y>
308 void umtv (const X& x, Y& y) const
309 {
310#ifdef DUNE_FMatrix_WITH_CHECKING
311 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
312 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
313#endif
314 for (size_type i=0; i<n; ++i)
315 y[i] += diag_[i] * x[i];
316 }
317
319 template<class X, class Y>
320 void umhv (const X& x, Y& y) const
321 {
322#ifdef DUNE_FMatrix_WITH_CHECKING
323 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
324 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
325#endif
326 for (size_type i=0; i<n; i++)
327 y[i] += conjugateComplex(diag_[i])*x[i];
328 }
329
331 template<class X, class Y>
332 void mmv (const X& x, Y& y) const
333 {
334#ifdef DUNE_FMatrix_WITH_CHECKING
335 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
336 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
337#endif
338 for (size_type i=0; i<n; ++i)
339 y[i] -= diag_[i] * x[i];
340 }
341
343 template<class X, class Y>
344 void mmtv (const X& x, Y& y) const
345 {
346#ifdef DUNE_FMatrix_WITH_CHECKING
347 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
348 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
349#endif
350 for (size_type i=0; i<n; ++i)
351 y[i] -= diag_[i] * x[i];
352 }
353
355 template<class X, class Y>
356 void mmhv (const X& x, Y& y) const
357 {
358#ifdef DUNE_FMatrix_WITH_CHECKING
359 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
360 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
361#endif
362 for (size_type i=0; i<n; i++)
363 y[i] -= conjugateComplex(diag_[i])*x[i];
364 }
365
367 template<class X, class Y>
368 void usmv (const typename FieldTraits<Y>::field_type & alpha,
369 const X& x, Y& y) const
370 {
371#ifdef DUNE_FMatrix_WITH_CHECKING
372 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
373 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
374#endif
375 for (size_type i=0; i<n; i++)
376 y[i] += alpha * diag_[i] * x[i];
377 }
378
380 template<class X, class Y>
381 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
382 const X& x, Y& y) const
383 {
384#ifdef DUNE_FMatrix_WITH_CHECKING
385 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
386 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
387#endif
388 for (size_type i=0; i<n; i++)
389 y[i] += alpha * diag_[i] * x[i];
390 }
391
393 template<class X, class Y>
394 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
395 const X& x, Y& y) const
396 {
397#ifdef DUNE_FMatrix_WITH_CHECKING
398 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
399 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
400#endif
401 for (size_type i=0; i<n; i++)
402 y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
403 }
404
405 //===== norms
406
408 double frobenius_norm () const
409 {
410 return diag_.two_norm();
411 }
412
414 double frobenius_norm2 () const
415 {
416 return diag_.two_norm2();
417 }
418
420 double infinity_norm () const
421 {
422 return diag_.infinity_norm();
423 }
424
426 double infinity_norm_real () const
427 {
428 return diag_.infinity_norm_real();
429 }
430
431
432
433 //===== solve
434
436 template<class V>
437 void solve (V& x, const V& b) const
438 {
439 for (int i=0; i<n; i++)
440 x[i] = b[i]/diag_[i];
441 }
442
444 void invert()
445 {
446 using real_type = typename FieldTraits<K>::real_type;
447 for (int i=0; i<n; i++)
448 diag_[i] = real_type(1.0)/diag_[i];
449 }
450
452 K determinant () const
453 {
454 K det = diag_[0];
455 for (int i=1; i<n; i++)
456 det *= diag_[i];
457 return det;
458 }
459
460
461
462 //===== matrix-matrix multiplication
463
466 template <class OtherScalar>
467 friend auto operator* ( const DiagonalMatrix& matrixA,
468 const DiagonalMatrix<OtherScalar, n>& matrixB)
469 {
471 for(int i=0; i<n; ++i)
472 result.diagonal(i) = matrixA.diagonal(i)*matrixB.diagonal(i);
473 return result;
474 }
475
476
477
478 //===== sizes
479
481 static constexpr size_type N ()
482 {
483 return n;
484 }
485
487 static constexpr size_type M ()
488 {
489 return n;
490 }
491
492
493
494 //===== query
495
497 bool exists (size_type i, size_type j) const
498 {
499 DUNE_ASSERT_BOUNDS(i >= 0 && i < n);
500 DUNE_ASSERT_BOUNDS(j >= 0 && j < n);
501 return i==j;
502 }
503
504
505
507 friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
508 {
509 for (size_type i=0; i<n; i++) {
510 for (size_type j=0; j<n; j++)
511 s << ((i==j) ? a.diag_[i] : 0) << " ";
512 s << std::endl;
513 }
514 return s;
515 }
516
519 {
520 return reference(const_cast<K*>(&diag_[i]), i);
521 }
522
524 const_reference operator[](size_type i) const
525 {
526 return const_reference(const_cast<K*>(&diag_[i]), i);
527 }
528
530 const K& diagonal(size_type i) const
531 {
532 return diag_[i];
533 }
534
537 {
538 return diag_[i];
539 }
540
543 {
544 return diag_;
545 }
546
549 {
550 return diag_;
551 }
552
553 private:
554
555 // the data, a FieldVector storing the diagonal
556 FieldVector<K,n> diag_;
557 };
558
559 template< class K, int n >
560 struct FieldTraits< DiagonalMatrix<K,n> >
561 {
562 typedef typename FieldTraits<K>::field_type field_type;
563 typedef typename FieldTraits<K>::real_type real_type;
564 };
565
566
567#ifndef DOXYGEN // hide specialization
570 template< class K >
571 class DiagonalMatrix<K, 1> : public FieldMatrix<K, 1, 1>
572 {
573 typedef FieldMatrix<K,1,1> Base;
574 public:
576 typedef typename Base::size_type size_type;
577
580 constexpr static int blocklevel = 1;
581
582 typedef typename Base::row_type row_type;
583
584 typedef typename Base::row_reference row_reference;
585 typedef typename Base::const_row_reference const_row_reference;
586
589 constexpr static int rows = 1;
592 constexpr static int cols = 1;
593
594
596 constexpr DiagonalMatrix() = default;
597
599 DiagonalMatrix(const K& scalar)
600 {
601 (*this)[0][0] = scalar;
602 }
603
605 const K& diagonal(size_type) const
606 {
607 return (*this)[0][0];
608 }
609
612 {
613 return (*this)[0][0];
614 }
615
617 const FieldVector<K,1>& diagonal() const
618 {
619 return (*this)[0];
620 }
621
623 FieldVector<K,1>& diagonal()
624 {
625 return (*this)[0];
626 }
627
629 DiagonalMatrix<K, 1> transposed() const
630 {
631 return *this;
632 }
633
636 template <class OtherScalar>
637 friend auto operator* ( const DiagonalMatrix& matrixA,
638 const DiagonalMatrix<OtherScalar, 1>& matrixB)
639 {
640 return DiagonalMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType, 1>{matrixA.diagonal(0)*matrixB.diagonal(0)};
641 }
642
643 };
644#endif
645
646
647 template<class DiagonalMatrixType>
648 class DiagonalMatrixWrapper
649 {
650 typedef typename DiagonalMatrixType::reference reference;
651 typedef typename DiagonalMatrixType::const_reference const_reference;
652 typedef typename DiagonalMatrixType::field_type K;
653 typedef DiagonalRowVector<K, DiagonalMatrixType::rows> row_type;
654 typedef std::size_t size_type;
655 typedef DiagonalMatrixWrapper< DiagonalMatrixType> MyType;
656
657 friend class ContainerWrapperIterator<const MyType, reference, reference>;
658 friend class ContainerWrapperIterator<const MyType, const_reference, const_reference>;
659
660 public:
661
662 DiagonalMatrixWrapper() :
663 mat_(0)
664 {}
665
666 DiagonalMatrixWrapper(const DiagonalMatrixType* mat) :
667 mat_(const_cast<DiagonalMatrixType*>(mat))
668 {}
669
670 size_type realIndex(int i) const
671 {
672 return i;
673 }
674
675 row_type* pointer(int i) const
676 {
677 row_ = row_type(&(mat_->diagonal(i)), i);
678 return &row_;
679 }
680
681 bool identical(const DiagonalMatrixWrapper& other) const
682 {
683 return mat_==other.mat_;
684 }
685
686 private:
687
688 mutable DiagonalMatrixType* mat_;
689 mutable row_type row_;
690 };
691
695 template< class K, int n >
696 class DiagonalRowVectorConst
697 {
698 template<class DiagonalMatrixType>
699 friend class DiagonalMatrixWrapper;
700 friend class ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&>;
701
702 public:
703 // remember size of vector
704 constexpr static int dimension = n;
705
706 // standard constructor and everything is sufficient ...
707
708 //===== type definitions and constants
709
711 typedef K field_type;
712
714 typedef K block_type;
715
717 typedef std::size_t size_type;
718
720 constexpr static int blocklevel = 1;
721
723 constexpr static int size = n;
724
727 p_(0),
728 row_(0)
729 {}
730
732 explicit DiagonalRowVectorConst (K* p, int col) :
733 p_(p),
734 row_(col)
735 {}
736
737 //===== access to components
738
740 const K& operator[] ([[maybe_unused]] size_type i) const
741 {
742 DUNE_ASSERT_BOUNDS(i == row_);
743 return *p_;
744 }
745
746 // check if row is identical to other row (not only identical values)
747 // since this is a proxy class we need to check equality of the stored pointer
748 bool identical(const DiagonalRowVectorConst<K,n>& other) const
749 {
750 return ((p_ == other.p_)and (row_ == other.row_));
751 }
752
757
760 {
761 return ConstIterator(*this,0);
762 }
763
766 {
767 return ConstIterator(*this,1);
768 }
769
773 {
774 return ConstIterator(*this,0);
775 }
776
780 {
781 return ConstIterator(*this,-1);
782 }
783
785 bool operator== (const DiagonalRowVectorConst& y) const
786 {
787 return ((p_==y.p_)and (row_==y.row_));
788 }
789
790 //===== sizes
791
793 size_type N () const
794 {
795 return n;
796 }
797
799 size_type dim () const
800 {
801 return n;
802 }
803
806 {
807 return row_;
808 }
809
811 const K& diagonal() const
812 {
813 return *p_;
814 }
815
816 protected:
817
818 size_type realIndex([[maybe_unused]] int i) const
819 {
820 return rowIndex();
821 }
822
823 K* pointer([[maybe_unused]] size_type i) const
824 {
825 return const_cast<K*>(p_);
826 }
827
828 DiagonalRowVectorConst* operator&()
829 {
830 return this;
831 }
832
833 // the data, very simply a pointer to the diagonal value and the row number
834 K* p_;
835 size_type row_;
836 };
837
838 template< class K, int n >
839 class DiagonalRowVector : public DiagonalRowVectorConst<K,n>
840 {
841 template<class DiagonalMatrixType>
842 friend class DiagonalMatrixWrapper;
843 friend class ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&>;
844
845 public:
846 // standard constructor and everything is sufficient ...
847
848 //===== type definitions and constants
849
851 typedef K field_type;
852
854 typedef K block_type;
855
857 typedef std::size_t size_type;
858
860 DiagonalRowVector() : DiagonalRowVectorConst<K,n>()
861 {}
862
864 explicit DiagonalRowVector (K* p, int col) : DiagonalRowVectorConst<K,n>(p, col)
865 {}
866
867 //===== assignment from scalar
869 DiagonalRowVector& operator= (const K& k)
870 {
871 *p_ = k;
872 return *this;
873 }
874
875 //===== access to components
876
878 K& operator[] ([[maybe_unused]] size_type i)
879 {
880 DUNE_ASSERT_BOUNDS(i == row_);
881 return *p_;
882 }
883
888
891 {
892 return Iterator(*this, 0);
893 }
894
897 {
898 return Iterator(*this, 1);
899 }
900
904 {
905 return Iterator(*this, 0);
906 }
907
911 {
912 return Iterator(*this, -1);
913 }
914
919
920 using DiagonalRowVectorConst<K,n>::identical;
921 using DiagonalRowVectorConst<K,n>::operator[];
922 using DiagonalRowVectorConst<K,n>::operator==;
923 using DiagonalRowVectorConst<K,n>::begin;
924 using DiagonalRowVectorConst<K,n>::end;
925 using DiagonalRowVectorConst<K,n>::beforeEnd;
926 using DiagonalRowVectorConst<K,n>::beforeBegin;
927 using DiagonalRowVectorConst<K,n>::N;
928 using DiagonalRowVectorConst<K,n>::dim;
929 using DiagonalRowVectorConst<K,n>::rowIndex;
930 using DiagonalRowVectorConst<K,n>::diagonal;
931
932 protected:
933
934 DiagonalRowVector* operator&()
935 {
936 return this;
937 }
938
939 private:
940
941 using DiagonalRowVectorConst<K,n>::p_;
942 using DiagonalRowVectorConst<K,n>::row_;
943 };
944
945
946 // implement type traits
947 template<class K, int n>
948 struct const_reference< DiagonalRowVector<K,n> >
949 {
950 typedef DiagonalRowVectorConst<K,n> type;
951 };
952
953 template<class K, int n>
954 struct const_reference< DiagonalRowVectorConst<K,n> >
955 {
956 typedef DiagonalRowVectorConst<K,n> type;
957 };
958
959 template<class K, int n>
960 struct mutable_reference< DiagonalRowVector<K,n> >
961 {
962 typedef DiagonalRowVector<K,n> type;
963 };
964
965 template<class K, int n>
966 struct mutable_reference< DiagonalRowVectorConst<K,n> >
967 {
968 typedef DiagonalRowVector<K,n> type;
969 };
970
971
972
995 template<class CW, class T, class R>
996 class ContainerWrapperIterator : public BidirectionalIteratorFacade<ContainerWrapperIterator<CW,T,R>,T, R, int>
997 {
998 typedef typename std::remove_const<CW>::type NonConstCW;
999
1000 friend class ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type>;
1001 friend class ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type>;
1002
1003 typedef ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type> MyType;
1004 typedef ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type> MyConstType;
1005
1006 public:
1007
1008 // Constructors needed by the facade iterators.
1010 containerWrapper_(),
1011 position_(0)
1012 {}
1013
1014 ContainerWrapperIterator(CW containerWrapper, int position) :
1015 containerWrapper_(containerWrapper),
1016 position_(position)
1017 {}
1018
1019 template<class OtherContainerWrapperIteratorType>
1020 ContainerWrapperIterator(OtherContainerWrapperIteratorType& other) :
1021 containerWrapper_(other.containerWrapper_),
1022 position_(other.position_)
1023 {}
1024
1025 ContainerWrapperIterator(const MyType& other) :
1026 containerWrapper_(other.containerWrapper_),
1027 position_(other.position_)
1028 {}
1029
1031 containerWrapper_(other.containerWrapper_),
1032 position_(other.position_)
1033 {}
1034
1035 template<class OtherContainerWrapperIteratorType>
1036 ContainerWrapperIterator& operator=(OtherContainerWrapperIteratorType& other)
1037 {
1038 containerWrapper_ = other.containerWrapper_;
1039 position_ = other.position_;
1040 return *this;
1041 }
1042
1043 // This operator is needed since we can not get the address of the
1044 // temporary object returned by dereference
1045 T* operator->() const
1046 {
1047 return containerWrapper_.pointer(position_);
1048 }
1049
1050 // Methods needed by the forward iterator
1051 bool equals(const MyType& other) const
1052 {
1053 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1054 }
1055
1056 bool equals(const MyConstType& other) const
1057 {
1058 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1059 }
1060
1061 R dereference() const
1062 {
1063 return *containerWrapper_.pointer(position_);
1064 }
1065
1066 void increment()
1067 {
1068 ++position_;
1069 }
1070
1071 // Additional function needed by BidirectionalIterator
1072 void decrement()
1073 {
1074 --position_;
1075 }
1076
1077 // Additional function needed by RandomAccessIterator
1078 R elementAt(int i) const
1079 {
1080 return *containerWrapper_.pointer(position_+i);
1081 }
1082
1083 void advance(int n)
1084 {
1085 position_=position_+n;
1086 }
1087
1088 template<class OtherContainerWrapperIteratorType>
1089 std::ptrdiff_t distanceTo(OtherContainerWrapperIteratorType& other) const
1090 {
1091 assert(containerWrapper_.identical(other));
1092 return other.position_ - position_;
1093 }
1094
1095 std::ptrdiff_t index() const
1096 {
1097 return containerWrapper_.realIndex(position_);
1098 }
1099
1100 private:
1101 NonConstCW containerWrapper_;
1102 size_t position_;
1103 };
1104
1105 template <class DenseMatrix, class field, int N>
1107 static void apply(DenseMatrix& denseMatrix,
1108 DiagonalMatrix<field, N> const& rhs) {
1109 DUNE_ASSERT_BOUNDS(denseMatrix.M() == N);
1110 DUNE_ASSERT_BOUNDS(denseMatrix.N() == N);
1111 denseMatrix = field(0);
1112 for (int i = 0; i < N; ++i)
1113 denseMatrix[i][i] = rhs.diagonal()[i];
1114 }
1115 };
1116 /* @} */
1117} // end namespace
1118#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:997
A dense n x m matrix.
Definition: densematrix.hh:140
constexpr size_type M() const
number of columns
Definition: densematrix.hh:703
constexpr size_type N() const
number of rows
Definition: densematrix.hh:697
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:650
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
Iterator begin()
begin iterator
Definition: densevector.hh:347
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: densevector.hh:661
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densevector.hh:677
A diagonal matrix of static size.
Definition: diagonalmatrix.hh:53
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:126
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.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
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:207
DiagonalMatrix< K, n > transposed() const
Return transposed of the matrix as DiagonalMatrix.
Definition: diagonalmatrix.hh:135
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: diagonalmatrix.hh:356
DiagonalMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: diagonalmatrix.hh:245
std::size_t size_type
The type used for the index access and size operations.
Definition: diagonalmatrix.hh:67
size_type dim() const
dimension of the vector space
Definition: diagonalmatrix.hh:799
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:182
static constexpr int rows
The number of rows.
Definition: diagonalmatrix.hh:81
static constexpr size_type M()
number of blocks in column direction
Definition: diagonalmatrix.hh:487
FieldVector< K, n > & diagonal()
Get reference to diagonal vector.
Definition: diagonalmatrix.hh:548
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: diagonalmatrix.hh:394
Iterator iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:887
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: diagonalmatrix.hh:70
K field_type
export the type representing the field
Definition: diagonalmatrix.hh:711
ConstIterator beforeEnd() const
Definition: diagonalmatrix.hh:772
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:184
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: diagonalmatrix.hh:497
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: diagonalmatrix.hh:178
DiagonalRowVector(K *p, int col)
Constructor making vector with identical coordinates.
Definition: diagonalmatrix.hh:864
K & diagonal(size_type i)
Get reference to diagonal entry.
Definition: diagonalmatrix.hh:536
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: diagonalmatrix.hh:437
static constexpr size_type N()
number of blocks in row direction
Definition: diagonalmatrix.hh:481
Iterator beforeBegin()
Definition: diagonalmatrix.hh:910
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: diagonalmatrix.hh:524
Iterator iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:144
ConstIterator begin() const
begin ConstIterator
Definition: diagonalmatrix.hh:759
ConstIterator const_iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:918
DiagonalMatrix & operator-=(const DiagonalMatrix &y)
vector space subtraction
Definition: diagonalmatrix.hh:224
DiagonalRowVectorConst(K *p, int col)
Constructor making vector with identical coordinates.
Definition: diagonalmatrix.hh:732
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: diagonalmatrix.hh:344
DiagonalMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: diagonalmatrix.hh:98
ContainerWrapperIterator< DiagonalRowVector< K, n >, K, K & > Iterator
Iterator class for sequential access.
Definition: diagonalmatrix.hh:885
Iterator beforeEnd()
Definition: diagonalmatrix.hh:903
void umtv(const X &x, Y &y) const
y += A^T x
Definition: diagonalmatrix.hh:308
ConstIterator const_iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:756
double infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: diagonalmatrix.hh:426
void umv(const X &x, Y &y) const
y += A x
Definition: diagonalmatrix.hh:296
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: diagonalmatrix.hh:142
void mv(const X &x, Y &y) const
y = A x
Definition: diagonalmatrix.hh:277
double frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: diagonalmatrix.hh:408
ConstIterator end() const
end iterator
Definition: diagonalmatrix.hh:193
static constexpr int cols
The number of columns.
Definition: diagonalmatrix.hh:83
size_type rowIndex() const
index of this row in surrounding matrix
Definition: diagonalmatrix.hh:805
bool operator!=(const DiagonalMatrix &other) const
incomparison operator
Definition: diagonalmatrix.hh:267
ConstIterator begin() const
begin iterator
Definition: diagonalmatrix.hh:187
Iterator begin()
begin iterator
Definition: diagonalmatrix.hh:890
void mmv(const X &x, Y &y) const
y -= A x
Definition: diagonalmatrix.hh:332
DiagonalMatrix & operator/=(const K &k)
vector space division by scalar
Definition: diagonalmatrix.hh:252
ConstIterator const_iterator
typedef for stl compliant access
Definition: diagonalmatrix.hh:180
bool identical(const DiagonalMatrix< K, n > &other) const
Check if matrix is the same object as the other matrix.
Definition: diagonalmatrix.hh:129
friend auto operator*(const DiagonalMatrix &matrixA, const DiagonalMatrix< OtherScalar, n > &matrixB)
Matrix-matrix multiplication.
Definition: diagonalmatrix.hh:467
Iterator end()
end iterator
Definition: diagonalmatrix.hh:157
DiagonalRowVector()
Constructor making uninitialized vector.
Definition: diagonalmatrix.hh:860
double frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: diagonalmatrix.hh:414
constexpr DiagonalMatrix()=default
Default constructor.
void mtv(const X &x, Y &y) const
y = A^T x
Definition: diagonalmatrix.hh:289
ContainerWrapperIterator< DiagonalRowVectorConst< K, n >, const K, const K & > ConstIterator
ConstIterator class for sequential access.
Definition: diagonalmatrix.hh:916
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: diagonalmatrix.hh:381
const FieldVector< K, n > & diagonal() const
Get const reference to diagonal vector.
Definition: diagonalmatrix.hh:542
void invert()
Compute inverse.
Definition: diagonalmatrix.hh:444
const K & diagonal(size_type i) const
Get const reference to diagonal entry.
Definition: diagonalmatrix.hh:530
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: diagonalmatrix.hh:793
DiagonalMatrix(std::initializer_list< K > const &l)
Construct diagonal matrix from an initializer list.
Definition: diagonalmatrix.hh:115
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:148
K field_type
export the type representing the field
Definition: diagonalmatrix.hh:851
reference operator[](size_type i)
Return reference object as row replacement.
Definition: diagonalmatrix.hh:518
Iterator end()
end iterator
Definition: diagonalmatrix.hh:896
DiagonalMatrix(const FieldVector< K, n > &diag)
Constructor initializing the diagonal with a vector.
Definition: diagonalmatrix.hh:103
std::size_t size_type
The type used for the index access and size operation.
Definition: diagonalmatrix.hh:857
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: diagonalmatrix.hh:368
double infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: diagonalmatrix.hh:420
DiagonalMatrix & operator+=(const DiagonalMatrix &y)
vector space addition
Definition: diagonalmatrix.hh:217
K block_type
export the type representing the components
Definition: diagonalmatrix.hh:64
void umhv(const X &x, Y &y) const
y += A^H x
Definition: diagonalmatrix.hh:320
Iterator beforeBegin()
Definition: diagonalmatrix.hh:171
ConstIterator end() const
end ConstIterator
Definition: diagonalmatrix.hh:765
std::size_t size_type
The type used for the index access and size operation.
Definition: diagonalmatrix.hh:717
DiagonalMatrix & operator=(const K &k)
Assignment from a scalar.
Definition: diagonalmatrix.hh:122
DiagonalRowVectorConst()
Constructor making uninitialized vector.
Definition: diagonalmatrix.hh:726
Iterator beforeEnd()
Definition: diagonalmatrix.hh:164
friend std::ostream & operator<<(std::ostream &s, const DiagonalMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: diagonalmatrix.hh:507
K value_type
export the type representing the field
Definition: diagonalmatrix.hh:60
K block_type
export the type representing the components
Definition: diagonalmatrix.hh:714
K determinant() const
calculates the determinant of this matrix
Definition: diagonalmatrix.hh:452
ConstIterator beforeBegin() const
Definition: diagonalmatrix.hh:779
ContainerWrapperIterator< DiagonalRowVectorConst< K, n >, const K, const K & > ConstIterator
ConstIterator class for sequential access.
Definition: diagonalmatrix.hh:754
Iterator begin()
begin iterator
Definition: diagonalmatrix.hh:151
bool operator==(const DiagonalMatrix &other) const
comparison operator
Definition: diagonalmatrix.hh:261
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: diagonalmatrix.hh:73
Iterator RowIterator
rename the iterators for easier access
Definition: diagonalmatrix.hh:146
K block_type
export the type representing the components
Definition: diagonalmatrix.hh:854
ConstIterator beforeEnd() const
Definition: diagonalmatrix.hh:200
const K & diagonal() const
the diagonal value
Definition: diagonalmatrix.hh:811
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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
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:59
get the 'mutable' version of a reference to a const object
Definition: genericiterator.hh:116
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)