Dune Core Modules (2.8.0)

matrix.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#ifndef DUNE_ISTL_MATRIX_HH
4#define DUNE_ISTL_MATRIX_HH
5
10#include <cmath>
11#include <memory>
12
17
18#include <dune/istl/bvector.hh>
19#include <dune/istl/istlexception.hh>
21
22namespace Dune {
23
24namespace MatrixImp
25{
37 template<class B, class A=std::allocator<B> >
38 class DenseMatrixBase : public Imp::block_vector_unmanaged<B,A>
39 // this derivation gives us all the blas level 1 and norms
40 // on the large array. However, access operators have to be
41 // overwritten.
42 {
43 public:
44
45 //===== type definitions and constants
46
48 using field_type = typename Imp::BlockTraits<B>::field_type;
49
51 typedef A allocator_type;
52
54 typedef typename A::size_type size_type;
55
62
66
67 // just a shorthand
68 typedef Imp::BlockVectorWindow<B,A> window_type;
69
70 typedef window_type reference;
71
72 typedef const window_type const_reference;
73
74
75 //===== constructors and such
76
80 DenseMatrixBase () : Imp::block_vector_unmanaged<B,A>()
81 {
82 // nothing is known ...
83 rows_ = 0;
84 columns_ = 0;
85 }
86
93 DenseMatrixBase (size_type rows, size_type columns) : Imp::block_vector_unmanaged<B,A>()
94 {
95 // and we can allocate the big array in the base class
96 this->n = rows*columns;
97 columns_ = columns;
98 if (this->n>0)
99 {
100 this->p = allocator_.allocate(this->n);
101 new (this->p)B[this->n];
102 }
103 else
104 {
105 this->n = 0;
106 this->p = 0;
107 }
108
109 // we can allocate the windows now
110 rows_ = rows;
111 }
112
115 {
116 // allocate the big array in the base class
117 this->n = a.n;
118 columns_ = a.columns_;
119 if (this->n>0)
120 {
121 // allocate and construct objects
122 this->p = allocator_.allocate(this->n);
123 new (this->p)B[this->n];
124
125 // copy data
126 for (size_type i=0; i<this->n; i++)
127 this->p[i]=a.p[i];
128 }
129 else
130 {
131 this->n = 0;
132 this->p = nullptr;
133 }
134
135 // we can allocate the windows now
136 rows_ = a.rows_;
137 }
138
141 {
142 if (this->n>0) {
143 size_type i=this->n;
144 while (i)
145 this->p[--i].~B();
146 allocator_.deallocate(this->p,this->n);
147 }
148 }
149
151 void resize (size_type rows, size_type columns)
152 {
153 // deconstruct objects and deallocate memory if necessary
154 if (this->n>0) {
155 size_type i=this->n;
156 while (i)
157 this->p[--i].~B();
158 allocator_.deallocate(this->p,this->n);
159 }
160
161 // and we can allocate the big array in the base class
162 this->n = rows*columns;
163 if (this->n>0)
164 {
165 this->p = allocator_.allocate(this->n);
166 new (this->p)B[this->n];
167 }
168 else
169 {
170 this->n = 0;
171 this->p = nullptr;
172 }
173
174 // we can allocate the windows now
175 rows_ = rows;
176 columns_ = columns;
177 }
178
181 {
182 if (&a!=this) // check if this and a are different objects
183 {
184 columns_ = a.columns_;
185 // reallocate arrays if necessary
186 // Note: still the block sizes may vary !
187 if (this->n!=a.n || rows_!=a.rows_)
188 {
189 // deconstruct objects and deallocate memory if necessary
190 if (this->n>0) {
191 size_type i=this->n;
192 while (i)
193 this->p[--i].~B();
194 allocator_.deallocate(this->p,this->n);
195 }
196
197 // allocate the big array in the base class
198 this->n = a.n;
199 if (this->n>0)
200 {
201 // allocate and construct objects
202 this->p = allocator_.allocate(this->n);
203 new (this->p)B[this->n];
204 }
205 else
206 {
207 this->n = 0;
208 this->p = nullptr;
209 }
210
211 // Copy number of rows
212 rows_ = a.rows_;
213 }
214
215 // and copy the data
216 for (size_type i=0; i<this->n; i++)
217 this->p[i]=a.p[i];
218 }
219
220 return *this;
221 }
222
223
224 //===== assignment from scalar
225
228 {
229 (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
230 return *this;
231 }
232
233
234 //===== access to components
235 // has to be overwritten from base class because it must
236 // return access to the windows
237
240 {
241#ifdef DUNE_ISTL_WITH_CHECKING
242 if (i>=rows_) DUNE_THROW(ISTLError,"index out of range");
243#endif
244 return window_type(this->p + i*columns_, columns_);
245 }
246
248 const_reference operator[] (size_type i) const
249 {
250#ifdef DUNE_ISTL_WITH_CHECKING
251 if (i<0 || i>=rows_) DUNE_THROW(ISTLError,"index out of range");
252#endif
253 return window_type(this->p + i*columns_, columns_);
254 }
255
256 // forward declaration
257 class ConstIterator;
258
261 {
262 public:
265 : window_(nullptr,0)
266 {
267 i = 0;
268 }
269
270 Iterator (Iterator& other) = default;
271 Iterator (Iterator&& other) = default;
272
274 Iterator (B* data, size_type columns, size_type _i)
275 : i(_i),
276 window_(data + _i*columns, columns)
277 {}
278
281 {
282 i = other.i;
283 // Do NOT use window_.operator=, because that copies the window content, not just the window!
284 window_.set(other.window_.getsize(),other.window_.getptr());
285 return *this;
286 }
287
290 {
291 i = other.i;
292 // Do NOT use window_.operator=, because that copies the window content, not just the window!
293 window_.set(other.window_.getsize(),other.window_.getptr());
294 return *this;
295 }
296
299 {
300 ++i;
301 window_.setptr(window_.getptr()+window_.getsize());
302 return *this;
303 }
304
307 {
308 --i;
309 window_.setptr(window_.getptr()-window_.getsize());
310 return *this;
311 }
312
314 bool operator== (const Iterator& it) const
315 {
316 return window_.getptr() == it.window_.getptr();
317 }
318
320 bool operator!= (const Iterator& it) const
321 {
322 return window_.getptr() != it.window_.getptr();
323 }
324
326 bool operator== (const ConstIterator& it) const
327 {
328 return window_.getptr() == it.window_.getptr();
329 }
330
332 bool operator!= (const ConstIterator& it) const
333 {
334 return window_.getptr() != it.window_.getptr();
335 }
336
338 window_type& operator* () const
339 {
340 return window_;
341 }
342
344 window_type* operator-> () const
345 {
346 return &window_;
347 }
348
349 // return index corresponding to pointer
350 size_type index () const
351 {
352 return i;
353 }
354
355 friend class ConstIterator;
356
357 private:
358 size_type i;
359 mutable window_type window_;
360 };
361
364 {
365 return Iterator(this->p, columns_, 0);
366 }
367
370 {
371 return Iterator(this->p, columns_, rows_);
372 }
373
377 {
378 return Iterator(this->p, columns_, rows_-1);
379 }
380
384 {
385 return Iterator(this->p, columns_, -1);
386 }
387
390 {
391 return Iterator(this->p, columns_, std::min(i,rows_));
392 }
393
396 {
397 return ConstIterator(this->p, columns_, std::min(i,rows_));
398 }
399
402 {
403 public:
406 : window_(nullptr,0)
407 {
408 i = 0;
409 }
410
412 ConstIterator (const B* data, size_type columns, size_type _i)
413 : i(_i),
414 window_(const_cast<B*>(data + _i * columns), columns)
415 {}
416
419 : i(it.i), window_(it.window_.getptr(),it.window_.getsize())
420 {}
421
422 ConstIterator& operator=(Iterator&& other)
423 {
424 i = other.i;
425 // Do NOT use window_.operator=, because that copies the window content, not just the window!
426 window_.set(other.window_.getsize(),other.window_.getptr());
427 return *this;
428 }
429
430 ConstIterator& operator=(Iterator& other)
431 {
432 i = other.i;
433 // Do NOT use window_.operator=, because that copies the window content, not just the window!
434 window_.set(other.window_.getsize(),other.window_.getptr());
435 return *this;
436 }
437
440 {
441 ++i;
442 window_.setptr(window_.getptr()+window_.getsize());
443 return *this;
444 }
445
448 {
449 --i;
450 window_.setptr(window_.getptr()-window_.getsize());
451 return *this;
452 }
453
455 bool operator== (const ConstIterator& it) const
456 {
457 return window_.getptr() == it.window_.getptr();
458 }
459
461 bool operator!= (const ConstIterator& it) const
462 {
463 return window_.getptr() != it.window_.getptr();
464 }
465
467 bool operator== (const Iterator& it) const
468 {
469 return window_.getptr() == it.window_.getptr();
470 }
471
473 bool operator!= (const Iterator& it) const
474 {
475 return window_.getptr() != it.window_.getptr();
476 }
477
479 const window_type& operator* () const
480 {
481 return window_;
482 }
483
485 const window_type* operator-> () const
486 {
487 return &window_;
488 }
489
490 // return index corresponding to pointer
491 size_type index () const
492 {
493 return i;
494 }
495
496 friend class Iterator;
497
498 private:
499 size_type i;
500 mutable window_type window_;
501 };
502
505
508
511 {
512 return ConstIterator(this->p, columns_, 0);
513 }
514
517 {
518 return ConstIterator(this->p, columns_, rows_);
519 }
520
524 {
525 return ConstIterator(this->p, columns_, rows_-1);
526 }
527
530 {
531 return ConstIterator(this->p, columns_, -1);
532 }
533
534 //===== sizes
535
537 size_type N () const
538 {
539 return rows_;
540 }
541
542
543 private:
544 size_type rows_; // number of matrix rows
545 size_type columns_; // number of matrix columns
546
547 A allocator_;
548 };
549
550} // namespace MatrixImp
551
557 template<class T, class A=std::allocator<T> >
558 class Matrix
559 {
560 public:
561
563 using field_type = typename Imp::BlockTraits<T>::field_type;
564
566 typedef T block_type;
567
569 typedef A allocator_type;
570
572 typedef typename MatrixImp::DenseMatrixBase<T,A>::window_type row_type;
573
575 typedef typename A::size_type size_type;
576
579
581 typedef typename row_type::iterator ColIterator;
582
585
587 typedef typename row_type::const_iterator ConstColIterator;
588
590 [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
591 static constexpr auto blocklevel = blockLevel<T>()+1;
592
594 Matrix() : data_(0,0), cols_(0)
595 {}
596
599 Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
600 {}
601
606 void setSize(size_type rows, size_type cols) {
607 data_.resize(rows,cols);
608 cols_ = cols;
609 }
610
613 {
614 return data_.begin();
615 }
616
619 {
620 return data_.end();
621 }
622
626 {
627 return data_.beforeEnd();
628 }
629
633 {
634 return data_.beforeBegin();
635 }
636
639 {
640 return data_.begin();
641 }
642
645 {
646 return data_.end();
647 }
648
652 {
653 return data_.beforeEnd();
654 }
655
659 {
660 return data_.beforeBegin();
661 }
662
665 {
666 data_ = t;
667 return *this;
668 }
669
672#ifdef DUNE_ISTL_WITH_CHECKING
673 if (row<0)
674 DUNE_THROW(ISTLError, "Can't access negative rows!");
675 if (row>=N())
676 DUNE_THROW(ISTLError, "Row index out of range!");
677#endif
678 return data_[row];
679 }
680
682 const row_type operator[](size_type row) const {
683#ifdef DUNE_ISTL_WITH_CHECKING
684 if (row<0)
685 DUNE_THROW(ISTLError, "Can't access negative rows!");
686 if (row>=N())
687 DUNE_THROW(ISTLError, "Row index out of range!");
688#endif
689 return data_[row];
690 }
691
693 size_type N() const {
694 return data_.N();
695 }
696
698 size_type M() const {
699 return cols_;
700 }
701
704 data_ *= scalar;
705 return (*this);
706 }
707
710 data_ /= scalar;
711 return (*this);
712 }
713
720#ifdef DUNE_ISTL_WITH_CHECKING
721 if(N()!=b.N() || M() != b.M())
722 DUNE_THROW(RangeError, "Matrix sizes do not match!");
723#endif
724 data_ += b.data_;
725 return (*this);
726 }
727
734#ifdef DUNE_ISTL_WITH_CHECKING
735 if(N()!=b.N() || M() != b.M())
736 DUNE_THROW(RangeError, "Matrix sizes do not match!");
737#endif
738 data_ -= b.data_;
739 return (*this);
740 }
741
744 Matrix out(M(), N());
745 for (size_type i=0; i<N(); i++)
746 for (size_type j=0; j<M(); j++)
747 out[j][i] = (*this)[i][j];
748
749 return out;
750 }
751
753 friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
754 Matrix<T> out(m1.N(), m2.M());
755 out = 0;
756
757 for (size_type i=0; i<out.N(); i++ ) {
758 for ( size_type j=0; j<out.M(); j++ )
759 for (size_type k=0; k<m1.M(); k++)
760 out[i][j] += m1[i][k]*m2[k][j];
761 }
762
763 return out;
764 }
765
767 template <class X, class Y>
768 friend Y operator*(const Matrix<T>& m, const X& vec) {
769#ifdef DUNE_ISTL_WITH_CHECKING
770 if (m.M()!=vec.size())
771 DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix columns!");
772#endif
773 Y out(m.N());
774 out = 0;
775
776 for (size_type i=0; i<out.size(); i++ ) {
777 for ( size_type j=0; j<vec.size(); j++ )
778 out[i] += m[i][j]*vec[j];
779 }
780
781 return out;
782 }
783
785 template <class X, class Y>
786 void mv(const X& x, Y& y) const
787 {
788#ifdef DUNE_ISTL_WITH_CHECKING
789 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
790 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
791#endif
792 for (size_type i=0; i<data_.N(); i++) {
793 y[i]=0;
794 for (size_type j=0; j<cols_; j++)
795 {
796 auto&& xj = Impl::asVector(x[j]);
797 auto&& yi = Impl::asVector(y[i]);
798 Impl::asMatrix((*this)[i][j]).umv(xj, yi);
799 }
800 }
801 }
802
804 template<class X, class Y>
805 void mtv (const X& x, Y& y) const
806 {
807#ifdef DUNE_ISTL_WITH_CHECKING
808 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
809 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
810#endif
811 for(size_type i=0; i<y.N(); ++i)
812 y[i]=0;
813 umtv(x,y);
814 }
815
817 template <class X, class Y>
818 void umv(const X& x, Y& y) const
819 {
820#ifdef DUNE_ISTL_WITH_CHECKING
821 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
822 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
823#endif
824 for (size_type i=0; i<data_.N(); i++)
825 for (size_type j=0; j<cols_; j++)
826 {
827 auto&& xj = Impl::asVector(x[j]);
828 auto&& yi = Impl::asVector(y[i]);
829 Impl::asMatrix((*this)[i][j]).umv(xj, yi);
830 }
831 }
832
834 template<class X, class Y>
835 void mmv (const X& x, Y& y) const
836 {
837#ifdef DUNE_ISTL_WITH_CHECKING
838 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
839 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
840#endif
841 for (size_type i=0; i<data_.N(); i++)
842 for (size_type j=0; j<cols_; j++)
843 {
844 auto&& xj = Impl::asVector(x[j]);
845 auto&& yi = Impl::asVector(y[i]);
846 Impl::asMatrix((*this)[i][j]).mmv(xj, yi);
847 }
848 }
849
851 template <class X, class Y>
852 void usmv(const field_type& alpha, const X& x, Y& y) const
853 {
854#ifdef DUNE_ISTL_WITH_CHECKING
855 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
856 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
857#endif
858 for (size_type i=0; i<data_.N(); i++)
859 for (size_type j=0; j<cols_; j++)
860 {
861 auto&& xj = Impl::asVector(x[j]);
862 auto&& yi = Impl::asVector(y[i]);
863 Impl::asMatrix((*this)[i][j]).usmv(alpha, xj, yi);
864 }
865 }
866
868 template<class X, class Y>
869 void umtv (const X& x, Y& y) const
870 {
871#ifdef DUNE_ISTL_WITH_CHECKING
872 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
873 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
874#endif
875 for (size_type i=0; i<data_.N(); i++)
876 for (size_type j=0; j<cols_; j++)
877 {
878 auto&& xi = Impl::asVector(x[i]);
879 auto&& yj = Impl::asVector(y[j]);
880 Impl::asMatrix((*this)[i][j]).umtv(xi, yj);
881 }
882 }
883
885 template<class X, class Y>
886 void mmtv (const X& x, Y& y) const
887 {
888#ifdef DUNE_ISTL_WITH_CHECKING
889 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
890 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
891#endif
892 for (size_type i=0; i<data_.N(); i++)
893 for (size_type j=0; j<cols_; j++)
894 {
895 auto&& xi = Impl::asVector(x[i]);
896 auto&& yj = Impl::asVector(y[j]);
897 Impl::asMatrix((*this)[i][j]).mmtv(xi, yj);
898 }
899 }
900
902 template<class X, class Y>
903 void usmtv (const field_type& alpha, const X& x, Y& y) const
904 {
905#ifdef DUNE_ISTL_WITH_CHECKING
906 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
907 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
908#endif
909 for (size_type i=0; i<data_.N(); i++)
910 for (size_type j=0; j<cols_; j++)
911 {
912 auto&& xi = Impl::asVector(x[i]);
913 auto&& yj = Impl::asVector(y[j]);
914 Impl::asMatrix((*this)[i][j]).usmtv(alpha, xi, yj);
915 }
916 }
917
919 template<class X, class Y>
920 void umhv (const X& x, Y& y) const
921 {
922#ifdef DUNE_ISTL_WITH_CHECKING
923 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
924 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
925#endif
926 for (size_type i=0; i<data_.N(); i++)
927 for (size_type j=0; j<cols_; j++)
928 {
929 auto&& xi = Impl::asVector(x[i]);
930 auto&& yj = Impl::asVector(y[j]);
931 Impl::asMatrix((*this)[i][j]).umhv(xi,yj);
932 }
933 }
934
936 template<class X, class Y>
937 void mmhv (const X& x, Y& y) const
938 {
939#ifdef DUNE_ISTL_WITH_CHECKING
940 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
941 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
942#endif
943 for (size_type i=0; i<data_.N(); i++)
944 for (size_type j=0; j<cols_; j++)
945 {
946 auto&& xi = Impl::asVector(x[i]);
947 auto&& yj = Impl::asVector(y[j]);
948 Impl::asMatrix((*this)[i][j]).mmhv(xi,yj);
949 }
950 }
951
953 template<class X, class Y>
954 void usmhv (const field_type& alpha, const X& x, Y& y) const
955 {
956#ifdef DUNE_ISTL_WITH_CHECKING
957 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
958 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
959#endif
960 for (size_type i=0; i<data_.N(); i++)
961 for (size_type j=0; j<cols_; j++)
962 {
963 auto&& xi = Impl::asVector(x[i]);
964 auto&& yj = Impl::asVector(y[j]);
965 Impl::asMatrix((*this)[i][j]).usmhv(alpha,xi,yj);
966 }
967 }
968
969 //===== norms
970
972 typename FieldTraits<field_type>::real_type frobenius_norm () const
973 {
974 return std::sqrt(frobenius_norm2());
975 }
976
978 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
979 {
980 typename FieldTraits<field_type>::real_type sum=0;
981 for (size_type i=0; i<this->N(); i++)
982 for (size_type j=0; j<this->M(); j++)
983 sum += Impl::asMatrix(data_[i][j]).frobenius_norm2();
984 return sum;
985 }
986
988 template <typename ft = field_type,
989 typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
990 typename FieldTraits<ft>::real_type infinity_norm() const {
991 using real_type = typename FieldTraits<ft>::real_type;
992 using std::max;
993
994 real_type norm = 0;
995 for (auto const &x : *this) {
996 real_type sum = 0;
997 for (auto const &y : x)
998 sum += Impl::asMatrix(y).infinity_norm();
999 norm = max(sum, norm);
1000 isNaN += sum;
1001 }
1002
1003 return norm;
1004 }
1005
1007 template <typename ft = field_type,
1008 typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
1009 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1010 using real_type = typename FieldTraits<ft>::real_type;
1011 using std::max;
1012
1013 real_type norm = 0;
1014 for (auto const &x : *this) {
1015 real_type sum = 0;
1016 for (auto const &y : x)
1017 sum += Impl::asMatrix(y).infinity_norm_real();
1018 norm = max(sum, norm);
1019 }
1020 return norm;
1021 }
1022
1024 template <typename ft = field_type,
1025 typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1026 typename FieldTraits<ft>::real_type infinity_norm() const {
1027 using real_type = typename FieldTraits<ft>::real_type;
1028 using std::max;
1029
1030 real_type norm = 0;
1031 real_type isNaN = 1;
1032 for (auto const &x : *this) {
1033 real_type sum = 0;
1034 for (auto const &y : x)
1035 sum += Impl::asMatrix(y).infinity_norm();
1036 norm = max(sum, norm);
1037 isNaN += sum;
1038 }
1039
1040 return norm * (isNaN / isNaN);
1041 }
1042
1044 template <typename ft = field_type,
1045 typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1046 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1047 using real_type = typename FieldTraits<ft>::real_type;
1048 using std::max;
1049
1050 real_type norm = 0;
1051 real_type isNaN = 1;
1052 for (auto const &x : *this) {
1053 real_type sum = 0;
1054 for (auto const &y : x)
1055 sum += Impl::asMatrix(y).infinity_norm_real();
1056 norm = max(sum, norm);
1057 isNaN += sum;
1058 }
1059
1060 return norm * (isNaN / isNaN);
1061 }
1062
1063 //===== query
1064
1066 bool exists ([[maybe_unused]] size_type i, [[maybe_unused]] size_type j) const
1067 {
1068#ifdef DUNE_ISTL_WITH_CHECKING
1069 if (i<0 || i>=N()) DUNE_THROW(ISTLError,"row index out of range");
1070 if (j<0 || i>=M()) DUNE_THROW(ISTLError,"column index out of range");
1071#endif
1072 return true;
1073 }
1074
1075 protected:
1076
1080
1087 };
1088
1090} // end namespace Dune
1091
1092#endif
Helper functions for determining the vector/matrix block level.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A vector of blocks with memory management.
Definition: bvector.hh:393
derive error class from the base class in common
Definition: istlexception.hh:17
ConstIterator class for sequential access.
Definition: matrix.hh:402
const window_type * operator->() const
arrow
Definition: matrix.hh:485
const window_type & operator*() const
dereferencing
Definition: matrix.hh:479
ConstIterator & operator++()
prefix increment
Definition: matrix.hh:439
ConstIterator(const B *data, size_type columns, size_type _i)
constructor from pointer
Definition: matrix.hh:412
ConstIterator & operator--()
prefix decrement
Definition: matrix.hh:447
ConstIterator(const Iterator &it)
constructor from non_const iterator
Definition: matrix.hh:418
bool operator!=(const ConstIterator &it) const
inequality
Definition: matrix.hh:461
ConstIterator()
constructor
Definition: matrix.hh:405
bool operator==(const ConstIterator &it) const
equality
Definition: matrix.hh:455
Iterator class for sequential access.
Definition: matrix.hh:261
Iterator & operator--()
prefix decrement
Definition: matrix.hh:306
bool operator!=(const Iterator &it) const
inequality
Definition: matrix.hh:320
Iterator & operator=(Iterator &&other)
Move assignment.
Definition: matrix.hh:280
Iterator & operator++()
prefix increment
Definition: matrix.hh:298
Iterator()
constructor, no arguments
Definition: matrix.hh:264
window_type & operator*() const
dereferencing
Definition: matrix.hh:338
bool operator==(const Iterator &it) const
equality
Definition: matrix.hh:314
Iterator & operator=(Iterator &other)
Copy assignment.
Definition: matrix.hh:289
Iterator(B *data, size_type columns, size_type _i)
constructor
Definition: matrix.hh:274
window_type * operator->() const
arrow
Definition: matrix.hh:344
A Vector of blocks with different blocksizes.
Definition: matrix.hh:42
BlockVector< B, A > block_type
Same as value_type, here for historical reasons.
Definition: matrix.hh:65
DenseMatrixBase & operator=(const DenseMatrixBase &a)
assignment
Definition: matrix.hh:180
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: matrix.hh:48
Iterator beforeBegin() const
Definition: matrix.hh:383
void resize(size_type rows, size_type columns)
same effect as constructor with same argument
Definition: matrix.hh:151
DenseMatrixBase(size_type rows, size_type columns)
Definition: matrix.hh:93
Iterator end()
end Iterator
Definition: matrix.hh:369
reference operator[](size_type i)
random access to blocks
Definition: matrix.hh:239
BlockVector< B, A > value_type
Type of the elements of the outer vector, i.e., dynamic vectors of B.
Definition: matrix.hh:61
Iterator find(size_type i)
random access returning iterator (end if not contained)
Definition: matrix.hh:389
size_type N() const
number of blocks in the vector (are of variable size here)
Definition: matrix.hh:537
ConstIterator beforeEnd() const
Definition: matrix.hh:523
ConstIterator end() const
end ConstIterator
Definition: matrix.hh:516
ConstIterator rend() const
end ConstIterator
Definition: matrix.hh:529
Iterator beforeEnd()
Definition: matrix.hh:376
ConstIterator begin() const
begin ConstIterator
Definition: matrix.hh:510
A allocator_type
export the allocator type
Definition: matrix.hh:51
DenseMatrixBase()
Definition: matrix.hh:80
ConstIterator find(size_type i) const
random access returning iterator (end if not contained)
Definition: matrix.hh:395
A::size_type size_type
The size type for the index access.
Definition: matrix.hh:54
~DenseMatrixBase()
free dynamic memory
Definition: matrix.hh:140
Iterator begin()
begin Iterator
Definition: matrix.hh:363
DenseMatrixBase(const DenseMatrixBase &a)
copy constructor, has copy semantics
Definition: matrix.hh:114
A generic dynamic dense matrix.
Definition: matrix.hh:559
size_type cols_
Number of columns of the matrix.
Definition: matrix.hh:1086
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:990
A allocator_type
Export the allocator.
Definition: matrix.hh:569
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:1009
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: matrix.hh:954
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: matrix.hh:852
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:575
MatrixImp::DenseMatrixBase< T, A > data_
Abuse DenseMatrixBase as an engine for a 2d array ISTL-style.
Definition: matrix.hh:1079
Matrix transpose() const
Return the transpose of the matrix.
Definition: matrix.hh:743
void mtv(const X &x, Y &y) const
y = A^T x
Definition: matrix.hh:805
void umv(const X &x, Y &y) const
y += A x
Definition: matrix.hh:818
void mv(const X &x, Y &y) const
y = A x
Definition: matrix.hh:786
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:584
void setSize(size_type rows, size_type cols)
Change the matrix size.
Definition: matrix.hh:606
RowIterator beforeBegin()
Definition: matrix.hh:632
RowIterator beforeEnd()
Definition: matrix.hh:625
Matrix()
Create empty matrix.
Definition: matrix.hh:594
Matrix & operator-=(const Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:733
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:978
row_type::iterator ColIterator
Iterator for the entries of each row.
Definition: matrix.hh:581
ConstRowIterator beforeEnd() const
Definition: matrix.hh:651
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:664
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:618
const row_type operator[](size_type row) const
The const index operator.
Definition: matrix.hh:682
ConstRowIterator end() const
Get const iterator to one beyond last row.
Definition: matrix.hh:644
friend Y operator*(const Matrix< T > &m, const X &vec)
Generic matrix-vector multiplication.
Definition: matrix.hh:768
Matrix< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:703
row_type operator[](size_type row)
The index operator.
Definition: matrix.hh:671
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:835
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:719
ConstRowIterator begin() const
Get const iterator to first row.
Definition: matrix.hh:638
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: matrix.hh:937
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:612
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:563
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
static constexpr auto blocklevel
The number of nesting levels the matrix contains.
Definition: matrix.hh:591
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
T block_type
Export the type representing the components.
Definition: matrix.hh:566
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:1066
Matrix< T > & operator/=(const field_type &scalar)
Division by a scalar.
Definition: matrix.hh:709
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:753
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:869
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:886
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:572
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:972
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: matrix.hh:903
void umhv(const X &x, Y &y) const
y += A^H x
Definition: matrix.hh:920
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
ConstRowIterator beforeBegin() const
Definition: matrix.hh:658
Matrix(size_type rows, size_type cols)
Create uninitialized matrix of size rows x cols.
Definition: matrix.hh:599
MatrixImp::DenseMatrixBase< T, A >::Iterator RowIterator
Iterator over the matrix rows.
Definition: matrix.hh:578
Default exception class for range errors.
Definition: exceptions.hh:252
Traits for type conversions and type information.
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:11
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)