DUNE PDELab (git)

matrix.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_MATRIX_HH
6#define DUNE_ISTL_MATRIX_HH
7
12#include <cmath>
13#include <memory>
14
19
20#include <dune/istl/bvector.hh>
21#include <dune/istl/istlexception.hh>
23
24namespace Dune {
25
26namespace MatrixImp
27{
39 template<class B, class A=std::allocator<B> >
40 class DenseMatrixBase : public Imp::block_vector_unmanaged<B,typename A::size_type>
41 // this derivation gives us all the blas level 1 and norms
42 // on the large array. However, access operators have to be
43 // overwritten.
44 {
45 public:
46
47 //===== type definitions and constants
48
50 using field_type = typename Imp::BlockTraits<B>::field_type;
51
53 typedef A allocator_type;
54
56 typedef typename A::size_type size_type;
57
64
68
69 // just a shorthand
70 typedef Imp::BlockVectorWindow<B,A> window_type;
71
72 typedef window_type reference;
73
74 typedef const window_type const_reference;
75
76
77 //===== constructors and such
78
82 DenseMatrixBase () : Imp::block_vector_unmanaged<B,size_type>()
83 {
84 // nothing is known ...
85 rows_ = 0;
86 columns_ = 0;
87 }
88
95 DenseMatrixBase (size_type rows, size_type columns) : Imp::block_vector_unmanaged<B,size_type>()
96 {
97 // and we can allocate the big array in the base class
98 this->n = rows*columns;
99 columns_ = columns;
100 if (this->n>0)
101 {
102 this->p = allocator_.allocate(this->n);
103 new (this->p)B[this->n];
104 }
105 else
106 {
107 this->n = 0;
108 this->p = 0;
109 }
110
111 // we can allocate the windows now
112 rows_ = rows;
113 }
114
117 {
118 // allocate the big array in the base class
119 this->n = a.n;
120 columns_ = a.columns_;
121 if (this->n>0)
122 {
123 // allocate and construct objects
124 this->p = allocator_.allocate(this->n);
125 new (this->p)B[this->n];
126
127 // copy data
128 for (size_type i=0; i<this->n; i++)
129 this->p[i]=a.p[i];
130 }
131 else
132 {
133 this->n = 0;
134 this->p = nullptr;
135 }
136
137 // we can allocate the windows now
138 rows_ = a.rows_;
139 }
140
143 {
144 if (this->n>0) {
145 size_type i=this->n;
146 while (i)
147 this->p[--i].~B();
148 allocator_.deallocate(this->p,this->n);
149 }
150 }
151
153 void resize (size_type rows, size_type columns)
154 {
155 // deconstruct objects and deallocate memory if necessary
156 if (this->n>0) {
157 size_type i=this->n;
158 while (i)
159 this->p[--i].~B();
160 allocator_.deallocate(this->p,this->n);
161 }
162
163 // and we can allocate the big array in the base class
164 this->n = rows*columns;
165 if (this->n>0)
166 {
167 this->p = allocator_.allocate(this->n);
168 new (this->p)B[this->n];
169 }
170 else
171 {
172 this->n = 0;
173 this->p = nullptr;
174 }
175
176 // we can allocate the windows now
177 rows_ = rows;
178 columns_ = columns;
179 }
180
183 {
184 if (&a!=this) // check if this and a are different objects
185 {
186 columns_ = a.columns_;
187 // reallocate arrays if necessary
188 // Note: still the block sizes may vary !
189 if (this->n!=a.n || rows_!=a.rows_)
190 {
191 // deconstruct objects and deallocate memory if necessary
192 if (this->n>0) {
193 size_type i=this->n;
194 while (i)
195 this->p[--i].~B();
196 allocator_.deallocate(this->p,this->n);
197 }
198
199 // allocate the big array in the base class
200 this->n = a.n;
201 if (this->n>0)
202 {
203 // allocate and construct objects
204 this->p = allocator_.allocate(this->n);
205 new (this->p)B[this->n];
206 }
207 else
208 {
209 this->n = 0;
210 this->p = nullptr;
211 }
212
213 // Copy number of rows
214 rows_ = a.rows_;
215 }
216
217 // and copy the data
218 for (size_type i=0; i<this->n; i++)
219 this->p[i]=a.p[i];
220 }
221
222 return *this;
223 }
224
225
226 //===== assignment from scalar
227
230 {
231 (static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
232 return *this;
233 }
234
235
236 //===== access to components
237 // has to be overwritten from base class because it must
238 // return access to the windows
239
242 {
243#ifdef DUNE_ISTL_WITH_CHECKING
244 if (i>=rows_) DUNE_THROW(ISTLError,"index out of range");
245#endif
246 return window_type(this->p + i*columns_, columns_);
247 }
248
250 const_reference operator[] (size_type i) const
251 {
252#ifdef DUNE_ISTL_WITH_CHECKING
253 if (i<0 || i>=rows_) DUNE_THROW(ISTLError,"index out of range");
254#endif
255 return window_type(this->p + i*columns_, columns_);
256 }
257
258 // forward declaration
259 class ConstIterator;
260
263 {
264 public:
267 : window_(nullptr,0)
268 {
269 i = 0;
270 }
271
272 Iterator (Iterator& other) = default;
273 Iterator (Iterator&& other) = default;
274
276 Iterator (B* data, size_type columns, size_type _i)
277 : i(_i),
278 window_(data + _i*columns, columns)
279 {}
280
283 {
284 i = other.i;
285 // Do NOT use window_.operator=, because that copies the window content, not just the window!
286 window_.set(other.window_.getsize(),other.window_.getptr());
287 return *this;
288 }
289
292 {
293 i = other.i;
294 // Do NOT use window_.operator=, because that copies the window content, not just the window!
295 window_.set(other.window_.getsize(),other.window_.getptr());
296 return *this;
297 }
298
301 {
302 ++i;
303 window_.setptr(window_.getptr()+window_.getsize());
304 return *this;
305 }
306
309 {
310 --i;
311 window_.setptr(window_.getptr()-window_.getsize());
312 return *this;
313 }
314
316 bool operator== (const Iterator& it) const
317 {
318 return window_.getptr() == it.window_.getptr();
319 }
320
322 bool operator!= (const Iterator& it) const
323 {
324 return window_.getptr() != it.window_.getptr();
325 }
326
328 bool operator== (const ConstIterator& it) const
329 {
330 return window_.getptr() == it.window_.getptr();
331 }
332
334 bool operator!= (const ConstIterator& it) const
335 {
336 return window_.getptr() != it.window_.getptr();
337 }
338
340 window_type& operator* () const
341 {
342 return window_;
343 }
344
346 window_type* operator-> () const
347 {
348 return &window_;
349 }
350
351 // return index corresponding to pointer
352 size_type index () const
353 {
354 return i;
355 }
356
357 friend class ConstIterator;
358
359 private:
360 size_type i;
361 mutable window_type window_;
362 };
363
366 {
367 return Iterator(this->p, columns_, 0);
368 }
369
372 {
373 return Iterator(this->p, columns_, rows_);
374 }
375
379 {
380 return Iterator(this->p, columns_, rows_-1);
381 }
382
386 {
387 return Iterator(this->p, columns_, -1);
388 }
389
392 {
393 return Iterator(this->p, columns_, std::min(i,rows_));
394 }
395
398 {
399 return ConstIterator(this->p, columns_, std::min(i,rows_));
400 }
401
404 {
405 public:
408 : window_(nullptr,0)
409 {
410 i = 0;
411 }
412
414 ConstIterator (const B* data, size_type columns, size_type _i)
415 : i(_i),
416 window_(const_cast<B*>(data + _i * columns), columns)
417 {}
418
421 : i(it.i), window_(it.window_.getptr(),it.window_.getsize())
422 {}
423
424 ConstIterator& operator=(Iterator&& other)
425 {
426 i = other.i;
427 // Do NOT use window_.operator=, because that copies the window content, not just the window!
428 window_.set(other.window_.getsize(),other.window_.getptr());
429 return *this;
430 }
431
432 ConstIterator& operator=(Iterator& other)
433 {
434 i = other.i;
435 // Do NOT use window_.operator=, because that copies the window content, not just the window!
436 window_.set(other.window_.getsize(),other.window_.getptr());
437 return *this;
438 }
439
442 {
443 ++i;
444 window_.setptr(window_.getptr()+window_.getsize());
445 return *this;
446 }
447
450 {
451 --i;
452 window_.setptr(window_.getptr()-window_.getsize());
453 return *this;
454 }
455
457 bool operator== (const ConstIterator& it) const
458 {
459 return window_.getptr() == it.window_.getptr();
460 }
461
463 bool operator!= (const ConstIterator& it) const
464 {
465 return window_.getptr() != it.window_.getptr();
466 }
467
469 bool operator== (const Iterator& it) const
470 {
471 return window_.getptr() == it.window_.getptr();
472 }
473
475 bool operator!= (const Iterator& it) const
476 {
477 return window_.getptr() != it.window_.getptr();
478 }
479
481 const window_type& operator* () const
482 {
483 return window_;
484 }
485
487 const window_type* operator-> () const
488 {
489 return &window_;
490 }
491
492 // return index corresponding to pointer
493 size_type index () const
494 {
495 return i;
496 }
497
498 friend class Iterator;
499
500 private:
501 size_type i;
502 mutable window_type window_;
503 };
504
507
510
513 {
514 return ConstIterator(this->p, columns_, 0);
515 }
516
519 {
520 return ConstIterator(this->p, columns_, rows_);
521 }
522
526 {
527 return ConstIterator(this->p, columns_, rows_-1);
528 }
529
532 {
533 return ConstIterator(this->p, columns_, -1);
534 }
535
536 //===== sizes
537
539 size_type N () const
540 {
541 return rows_;
542 }
543
544
545 private:
546 size_type rows_; // number of matrix rows
547 size_type columns_; // number of matrix columns
548
549 A allocator_;
550 };
551
552} // namespace MatrixImp
553
559 template<class T, class A=std::allocator<T> >
560 class Matrix
561 {
562 public:
563
565 using field_type = typename Imp::BlockTraits<T>::field_type;
566
568 typedef T block_type;
569
571 typedef A allocator_type;
572
574 typedef typename MatrixImp::DenseMatrixBase<T,A>::window_type row_type;
575
577 typedef typename A::size_type size_type;
578
581
583 typedef typename row_type::iterator ColIterator;
584
587
589 typedef typename row_type::const_iterator ConstColIterator;
590
592 Matrix() : data_(0,0), cols_(0)
593 {}
594
597 Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
598 {}
599
604 void setSize(size_type rows, size_type cols) {
605 data_.resize(rows,cols);
606 cols_ = cols;
607 }
608
611 {
612 return data_.begin();
613 }
614
617 {
618 return data_.end();
619 }
620
624 {
625 return data_.beforeEnd();
626 }
627
631 {
632 return data_.beforeBegin();
633 }
634
637 {
638 return data_.begin();
639 }
640
643 {
644 return data_.end();
645 }
646
650 {
651 return data_.beforeEnd();
652 }
653
657 {
658 return data_.beforeBegin();
659 }
660
663 {
664 data_ = t;
665 return *this;
666 }
667
670#ifdef DUNE_ISTL_WITH_CHECKING
671 if (row<0)
672 DUNE_THROW(ISTLError, "Can't access negative rows!");
673 if (row>=N())
674 DUNE_THROW(ISTLError, "Row index out of range!");
675#endif
676 return data_[row];
677 }
678
680 const row_type operator[](size_type row) const {
681#ifdef DUNE_ISTL_WITH_CHECKING
682 if (row<0)
683 DUNE_THROW(ISTLError, "Can't access negative rows!");
684 if (row>=N())
685 DUNE_THROW(ISTLError, "Row index out of range!");
686#endif
687 return data_[row];
688 }
689
691 size_type N() const {
692 return data_.N();
693 }
694
696 size_type M() const {
697 return cols_;
698 }
699
702 data_ *= scalar;
703 return (*this);
704 }
705
708 data_ /= scalar;
709 return (*this);
710 }
711
718#ifdef DUNE_ISTL_WITH_CHECKING
719 if(N()!=b.N() || M() != b.M())
720 DUNE_THROW(RangeError, "Matrix sizes do not match!");
721#endif
722 data_ += b.data_;
723 return (*this);
724 }
725
732#ifdef DUNE_ISTL_WITH_CHECKING
733 if(N()!=b.N() || M() != b.M())
734 DUNE_THROW(RangeError, "Matrix sizes do not match!");
735#endif
736 data_ -= b.data_;
737 return (*this);
738 }
739
742 Matrix out(M(), N());
743 for (size_type i=0; i<N(); i++)
744 for (size_type j=0; j<M(); j++)
745 out[j][i] = (*this)[i][j];
746
747 return out;
748 }
749
751 friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
752 Matrix<T> out(m1.N(), m2.M());
753 out = 0;
754
755 for (size_type i=0; i<out.N(); i++ ) {
756 for ( size_type j=0; j<out.M(); j++ )
757 for (size_type k=0; k<m1.M(); k++)
758 out[i][j] += m1[i][k]*m2[k][j];
759 }
760
761 return out;
762 }
763
765 template <class X, class Y>
766 friend Y operator*(const Matrix<T>& m, const X& vec) {
767#ifdef DUNE_ISTL_WITH_CHECKING
768 if (m.M()!=vec.size())
769 DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix columns!");
770#endif
771 Y out(m.N());
772 out = 0;
773
774 for (size_type i=0; i<out.size(); i++ ) {
775 for ( size_type j=0; j<vec.size(); j++ )
776 out[i] += m[i][j]*vec[j];
777 }
778
779 return out;
780 }
781
783 template <class X, class Y>
784 void mv(const X& x, Y& y) const
785 {
786#ifdef DUNE_ISTL_WITH_CHECKING
787 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
788 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
789#endif
790 for (size_type i=0; i<data_.N(); i++) {
791 y[i]=0;
792 for (size_type j=0; j<cols_; j++)
793 {
794 auto&& xj = Impl::asVector(x[j]);
795 auto&& yi = Impl::asVector(y[i]);
796 Impl::asMatrix((*this)[i][j]).umv(xj, yi);
797 }
798 }
799 }
800
802 template<class X, class Y>
803 void mtv (const X& x, Y& y) const
804 {
805#ifdef DUNE_ISTL_WITH_CHECKING
806 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
807 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
808#endif
809 for(size_type i=0; i<y.N(); ++i)
810 y[i]=0;
811 umtv(x,y);
812 }
813
815 template <class X, class Y>
816 void umv(const X& x, Y& y) const
817 {
818#ifdef DUNE_ISTL_WITH_CHECKING
819 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
820 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
821#endif
822 for (size_type i=0; i<data_.N(); i++)
823 for (size_type j=0; j<cols_; j++)
824 {
825 auto&& xj = Impl::asVector(x[j]);
826 auto&& yi = Impl::asVector(y[i]);
827 Impl::asMatrix((*this)[i][j]).umv(xj, yi);
828 }
829 }
830
832 template<class X, class Y>
833 void mmv (const X& x, Y& y) const
834 {
835#ifdef DUNE_ISTL_WITH_CHECKING
836 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
837 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
838#endif
839 for (size_type i=0; i<data_.N(); i++)
840 for (size_type j=0; j<cols_; j++)
841 {
842 auto&& xj = Impl::asVector(x[j]);
843 auto&& yi = Impl::asVector(y[i]);
844 Impl::asMatrix((*this)[i][j]).mmv(xj, yi);
845 }
846 }
847
849 template <class X, class Y>
850 void usmv(const field_type& alpha, const X& x, Y& y) const
851 {
852#ifdef DUNE_ISTL_WITH_CHECKING
853 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
854 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
855#endif
856 for (size_type i=0; i<data_.N(); i++)
857 for (size_type j=0; j<cols_; j++)
858 {
859 auto&& xj = Impl::asVector(x[j]);
860 auto&& yi = Impl::asVector(y[i]);
861 Impl::asMatrix((*this)[i][j]).usmv(alpha, xj, yi);
862 }
863 }
864
866 template<class X, class Y>
867 void umtv (const X& x, Y& y) const
868 {
869#ifdef DUNE_ISTL_WITH_CHECKING
870 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
871 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
872#endif
873 for (size_type i=0; i<data_.N(); i++)
874 for (size_type j=0; j<cols_; j++)
875 {
876 auto&& xi = Impl::asVector(x[i]);
877 auto&& yj = Impl::asVector(y[j]);
878 Impl::asMatrix((*this)[i][j]).umtv(xi, yj);
879 }
880 }
881
883 template<class X, class Y>
884 void mmtv (const X& x, Y& y) const
885 {
886#ifdef DUNE_ISTL_WITH_CHECKING
887 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
888 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
889#endif
890 for (size_type i=0; i<data_.N(); i++)
891 for (size_type j=0; j<cols_; j++)
892 {
893 auto&& xi = Impl::asVector(x[i]);
894 auto&& yj = Impl::asVector(y[j]);
895 Impl::asMatrix((*this)[i][j]).mmtv(xi, yj);
896 }
897 }
898
900 template<class X, class Y>
901 void usmtv (const field_type& alpha, const X& x, Y& y) const
902 {
903#ifdef DUNE_ISTL_WITH_CHECKING
904 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
905 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
906#endif
907 for (size_type i=0; i<data_.N(); i++)
908 for (size_type j=0; j<cols_; j++)
909 {
910 auto&& xi = Impl::asVector(x[i]);
911 auto&& yj = Impl::asVector(y[j]);
912 Impl::asMatrix((*this)[i][j]).usmtv(alpha, xi, yj);
913 }
914 }
915
917 template<class X, class Y>
918 void umhv (const X& x, Y& y) const
919 {
920#ifdef DUNE_ISTL_WITH_CHECKING
921 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
922 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
923#endif
924 for (size_type i=0; i<data_.N(); i++)
925 for (size_type j=0; j<cols_; j++)
926 {
927 auto&& xi = Impl::asVector(x[i]);
928 auto&& yj = Impl::asVector(y[j]);
929 Impl::asMatrix((*this)[i][j]).umhv(xi,yj);
930 }
931 }
932
934 template<class X, class Y>
935 void mmhv (const X& x, Y& y) const
936 {
937#ifdef DUNE_ISTL_WITH_CHECKING
938 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
939 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
940#endif
941 for (size_type i=0; i<data_.N(); i++)
942 for (size_type j=0; j<cols_; j++)
943 {
944 auto&& xi = Impl::asVector(x[i]);
945 auto&& yj = Impl::asVector(y[j]);
946 Impl::asMatrix((*this)[i][j]).mmhv(xi,yj);
947 }
948 }
949
951 template<class X, class Y>
952 void usmhv (const field_type& alpha, const X& x, Y& y) const
953 {
954#ifdef DUNE_ISTL_WITH_CHECKING
955 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
956 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
957#endif
958 for (size_type i=0; i<data_.N(); i++)
959 for (size_type j=0; j<cols_; j++)
960 {
961 auto&& xi = Impl::asVector(x[i]);
962 auto&& yj = Impl::asVector(y[j]);
963 Impl::asMatrix((*this)[i][j]).usmhv(alpha,xi,yj);
964 }
965 }
966
967 //===== norms
968
970 typename FieldTraits<field_type>::real_type frobenius_norm () const
971 {
972 return std::sqrt(frobenius_norm2());
973 }
974
976 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
977 {
978 typename FieldTraits<field_type>::real_type sum=0;
979 for (size_type i=0; i<this->N(); i++)
980 for (size_type j=0; j<this->M(); j++)
981 sum += Impl::asMatrix(data_[i][j]).frobenius_norm2();
982 return sum;
983 }
984
986 template <typename ft = field_type,
987 typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
988 typename FieldTraits<ft>::real_type infinity_norm() const {
989 using real_type = typename FieldTraits<ft>::real_type;
990 using std::max;
991
992 real_type norm = 0;
993 for (auto const &x : *this) {
994 real_type sum = 0;
995 for (auto const &y : x)
996 sum += Impl::asMatrix(y).infinity_norm();
997 norm = max(sum, norm);
998 isNaN += sum;
999 }
1000
1001 return norm;
1002 }
1003
1005 template <typename ft = field_type,
1006 typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
1007 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1008 using real_type = typename FieldTraits<ft>::real_type;
1009 using std::max;
1010
1011 real_type norm = 0;
1012 for (auto const &x : *this) {
1013 real_type sum = 0;
1014 for (auto const &y : x)
1015 sum += Impl::asMatrix(y).infinity_norm_real();
1016 norm = max(sum, norm);
1017 }
1018 return norm;
1019 }
1020
1022 template <typename ft = field_type,
1023 typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1024 typename FieldTraits<ft>::real_type infinity_norm() const {
1025 using real_type = typename FieldTraits<ft>::real_type;
1026 using std::max;
1027
1028 real_type norm = 0;
1029 real_type isNaN = 1;
1030 for (auto const &x : *this) {
1031 real_type sum = 0;
1032 for (auto const &y : x)
1033 sum += Impl::asMatrix(y).infinity_norm();
1034 norm = max(sum, norm);
1035 isNaN += sum;
1036 }
1037
1038 return norm * (isNaN / isNaN);
1039 }
1040
1042 template <typename ft = field_type,
1043 typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1044 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1045 using real_type = typename FieldTraits<ft>::real_type;
1046 using std::max;
1047
1048 real_type norm = 0;
1049 real_type isNaN = 1;
1050 for (auto const &x : *this) {
1051 real_type sum = 0;
1052 for (auto const &y : x)
1053 sum += Impl::asMatrix(y).infinity_norm_real();
1054 norm = max(sum, norm);
1055 isNaN += sum;
1056 }
1057
1058 return norm * (isNaN / isNaN);
1059 }
1060
1061 //===== query
1062
1064 bool exists ([[maybe_unused]] size_type i, [[maybe_unused]] size_type j) const
1065 {
1066#ifdef DUNE_ISTL_WITH_CHECKING
1067 if (i<0 || i>=N()) DUNE_THROW(ISTLError,"row index out of range");
1068 if (j<0 || i>=M()) DUNE_THROW(ISTLError,"column index out of range");
1069#endif
1070 return true;
1071 }
1072
1073 protected:
1074
1078
1085 };
1086
1087 template<class T, class A>
1088 struct FieldTraits< Matrix<T, A> >
1089 {
1090 using field_type = typename Matrix<T, A>::field_type;
1091 using real_type = typename FieldTraits<field_type>::real_type;
1092 };
1093
1095} // end namespace Dune
1096
1097#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:392
derive error class from the base class in common
Definition: istlexception.hh:19
ConstIterator class for sequential access.
Definition: matrix.hh:404
const window_type * operator->() const
arrow
Definition: matrix.hh:487
const window_type & operator*() const
dereferencing
Definition: matrix.hh:481
ConstIterator & operator++()
prefix increment
Definition: matrix.hh:441
ConstIterator(const B *data, size_type columns, size_type _i)
constructor from pointer
Definition: matrix.hh:414
ConstIterator & operator--()
prefix decrement
Definition: matrix.hh:449
ConstIterator(const Iterator &it)
constructor from non_const iterator
Definition: matrix.hh:420
bool operator!=(const ConstIterator &it) const
inequality
Definition: matrix.hh:463
ConstIterator()
constructor
Definition: matrix.hh:407
bool operator==(const ConstIterator &it) const
equality
Definition: matrix.hh:457
Iterator class for sequential access.
Definition: matrix.hh:263
Iterator & operator--()
prefix decrement
Definition: matrix.hh:308
bool operator!=(const Iterator &it) const
inequality
Definition: matrix.hh:322
Iterator & operator=(Iterator &&other)
Move assignment.
Definition: matrix.hh:282
Iterator & operator++()
prefix increment
Definition: matrix.hh:300
Iterator()
constructor, no arguments
Definition: matrix.hh:266
window_type & operator*() const
dereferencing
Definition: matrix.hh:340
bool operator==(const Iterator &it) const
equality
Definition: matrix.hh:316
Iterator & operator=(Iterator &other)
Copy assignment.
Definition: matrix.hh:291
Iterator(B *data, size_type columns, size_type _i)
constructor
Definition: matrix.hh:276
window_type * operator->() const
arrow
Definition: matrix.hh:346
A Vector of blocks with different blocksizes.
Definition: matrix.hh:44
BlockVector< B, A > block_type
Same as value_type, here for historical reasons.
Definition: matrix.hh:67
DenseMatrixBase & operator=(const DenseMatrixBase &a)
assignment
Definition: matrix.hh:182
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: matrix.hh:50
Iterator beforeBegin() const
Definition: matrix.hh:385
void resize(size_type rows, size_type columns)
same effect as constructor with same argument
Definition: matrix.hh:153
DenseMatrixBase(size_type rows, size_type columns)
Definition: matrix.hh:95
Iterator end()
end Iterator
Definition: matrix.hh:371
reference operator[](size_type i)
random access to blocks
Definition: matrix.hh:241
BlockVector< B, A > value_type
Type of the elements of the outer vector, i.e., dynamic vectors of B.
Definition: matrix.hh:63
Iterator find(size_type i)
random access returning iterator (end if not contained)
Definition: matrix.hh:391
size_type N() const
number of blocks in the vector (are of variable size here)
Definition: matrix.hh:539
ConstIterator beforeEnd() const
Definition: matrix.hh:525
ConstIterator end() const
end ConstIterator
Definition: matrix.hh:518
ConstIterator rend() const
end ConstIterator
Definition: matrix.hh:531
Iterator beforeEnd()
Definition: matrix.hh:378
ConstIterator begin() const
begin ConstIterator
Definition: matrix.hh:512
A allocator_type
export the allocator type
Definition: matrix.hh:53
DenseMatrixBase()
Definition: matrix.hh:82
ConstIterator find(size_type i) const
random access returning iterator (end if not contained)
Definition: matrix.hh:397
A::size_type size_type
The size type for the index access.
Definition: matrix.hh:56
~DenseMatrixBase()
free dynamic memory
Definition: matrix.hh:142
Iterator begin()
begin Iterator
Definition: matrix.hh:365
DenseMatrixBase(const DenseMatrixBase &a)
copy constructor, has copy semantics
Definition: matrix.hh:116
A generic dynamic dense matrix.
Definition: matrix.hh:561
size_type cols_
Number of columns of the matrix.
Definition: matrix.hh:1084
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:988
A allocator_type
Export the allocator.
Definition: matrix.hh:571
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:1007
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: matrix.hh:952
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: matrix.hh:850
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:577
MatrixImp::DenseMatrixBase< T, A > data_
Abuse DenseMatrixBase as an engine for a 2d array ISTL-style.
Definition: matrix.hh:1077
Matrix transpose() const
Return the transpose of the matrix.
Definition: matrix.hh:741
void mtv(const X &x, Y &y) const
y = A^T x
Definition: matrix.hh:803
void umv(const X &x, Y &y) const
y += A x
Definition: matrix.hh:816
void mv(const X &x, Y &y) const
y = A x
Definition: matrix.hh:784
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:586
void setSize(size_type rows, size_type cols)
Change the matrix size.
Definition: matrix.hh:604
RowIterator beforeBegin()
Definition: matrix.hh:630
RowIterator beforeEnd()
Definition: matrix.hh:623
Matrix()
Create empty matrix.
Definition: matrix.hh:592
Matrix & operator-=(const Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:731
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:976
row_type::iterator ColIterator
Iterator for the entries of each row.
Definition: matrix.hh:583
ConstRowIterator beforeEnd() const
Definition: matrix.hh:649
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:662
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:616
const row_type operator[](size_type row) const
The const index operator.
Definition: matrix.hh:680
ConstRowIterator end() const
Get const iterator to one beyond last row.
Definition: matrix.hh:642
friend Y operator*(const Matrix< T > &m, const X &vec)
Generic matrix-vector multiplication.
Definition: matrix.hh:766
Matrix< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:701
row_type operator[](size_type row)
The index operator.
Definition: matrix.hh:669
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:833
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:717
ConstRowIterator begin() const
Get const iterator to first row.
Definition: matrix.hh:636
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: matrix.hh:935
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:610
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
size_type M() const
Return the number of columns.
Definition: matrix.hh:696
T block_type
Export the type representing the components.
Definition: matrix.hh:568
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:1064
Matrix< T > & operator/=(const field_type &scalar)
Division by a scalar.
Definition: matrix.hh:707
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:751
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:867
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:884
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:574
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:970
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: matrix.hh:901
void umhv(const X &x, Y &y) const
y += A^H x
Definition: matrix.hh:918
size_type N() const
Return the number of rows.
Definition: matrix.hh:691
ConstRowIterator beforeBegin() const
Definition: matrix.hh:656
Matrix(size_type rows, size_type cols)
Create uninitialized matrix of size rows x cols.
Definition: matrix.hh:597
MatrixImp::DenseMatrixBase< T, A >::Iterator RowIterator
Iterator over the matrix rows.
Definition: matrix.hh:580
Default exception class for range errors.
Definition: exceptions.hh:254
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:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:627
Dune namespace.
Definition: alignedallocator.hh:13
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 23, 23:29, 2024)