Go to the documentation of this file.
3 #ifndef DUNE_ISTL_MATRIX_HH
4 #define DUNE_ISTL_MATRIX_HH
13 #include <dune/common/ftraits.hh>
14 #include <dune/common/typetraits.hh>
15 #include <dune/common/scalarvectorview.hh>
16 #include <dune/common/scalarmatrixview.hh>
36 template<
class B,
class A=std::allocator<B> >
47 using field_type =
typename Imp::BlockTraits<B>::field_type;
95 this->n = rows*columns;
99 this->p = allocator_.allocate(this->n);
100 new (this->p)B[this->n];
117 columns_ = a.columns_;
121 this->p = allocator_.allocate(this->n);
122 new (this->p)B[this->n];
145 allocator_.deallocate(this->p,this->n);
157 allocator_.deallocate(this->p,this->n);
161 this->n = rows*columns;
164 this->p = allocator_.allocate(this->n);
165 new (this->p)B[this->n];
183 columns_ = a.columns_;
186 if (this->n!=a.n || rows_!=a.rows_)
193 allocator_.deallocate(this->p,this->n);
201 this->p = allocator_.allocate(this->n);
202 new (this->p)B[this->n];
228 (
static_cast<Imp::block_vector_unmanaged<B,A>&
>(*this)) = k;
240 #ifdef DUNE_ISTL_WITH_CHECKING
241 if (i>=rows_) DUNE_THROW(
ISTLError,
"index out of range");
243 return window_type(this->p + i*columns_, columns_);
249 #ifdef DUNE_ISTL_WITH_CHECKING
250 if (i<0 || i>=rows_) DUNE_THROW(
ISTLError,
"index out of range");
252 return window_type(this->p + i*columns_, columns_);
275 window_(data + _i*columns, columns)
283 window_.set(other.window_.getsize(),other.window_.getptr());
292 window_.set(other.window_.getsize(),other.window_.getptr());
300 window_.setptr(window_.getptr()+window_.getsize());
308 window_.setptr(window_.getptr()-window_.getsize());
315 return window_.getptr() == it.window_.getptr();
321 return window_.getptr() != it.window_.getptr();
327 return window_.getptr() == it.window_.getptr();
333 return window_.getptr() != it.window_.getptr();
364 return Iterator(this->p, columns_, 0);
370 return Iterator(this->p, columns_, rows_);
377 return Iterator(this->p, columns_, rows_-1);
384 return Iterator(this->p, columns_, -1);
390 return Iterator(this->p, columns_, std::min(i,rows_));
396 return ConstIterator(this->p, columns_, std::min(i,rows_));
413 window_(const_cast<B*>(data + _i * columns), columns)
418 : i(it.i), window_(it.window_.getptr(),it.window_.getsize())
425 window_.set(other.window_.getsize(),other.window_.getptr());
433 window_.set(other.window_.getsize(),other.window_.getptr());
441 window_.setptr(window_.getptr()+window_.getsize());
449 window_.setptr(window_.getptr()-window_.getsize());
456 return window_.getptr() == it.window_.getptr();
462 return window_.getptr() != it.window_.getptr();
468 return window_.getptr() == it.window_.getptr();
474 return window_.getptr() != it.window_.getptr();
511 return ConstIterator(this->p, columns_, 0);
515 ConstIterator
end ()
const
517 return ConstIterator(this->p, columns_, rows_);
524 return ConstIterator(this->p, columns_, rows_-1);
530 return ConstIterator(this->p, columns_, -1);
556 template<
class T,
class A=std::allocator<T> >
589 static constexpr
unsigned int blocklevel = Imp::BlockTraits<T>::blockLevel()+1;
605 data_.resize(rows,cols);
612 return data_.begin();
625 return data_.beforeEnd();
632 return data_.beforeBegin();
638 return data_.begin();
651 return data_.beforeEnd();
658 return data_.beforeBegin();
670 #ifdef DUNE_ISTL_WITH_CHECKING
672 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
674 DUNE_THROW(
ISTLError,
"Row index out of range!");
681 #ifdef DUNE_ISTL_WITH_CHECKING
683 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
685 DUNE_THROW(
ISTLError,
"Row index out of range!");
718 #ifdef DUNE_ISTL_WITH_CHECKING
719 if(
N()!=b.
N() ||
M() != b.
M())
720 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
732 #ifdef DUNE_ISTL_WITH_CHECKING
733 if(
N()!=b.
N() ||
M() != b.
M())
734 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
745 out[j][i] = (*
this)[i][j];
758 out[i][j] += m1[i][k]*m2[k][j];
765 template <
class X,
class Y>
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!");
774 for (
size_type i=0; i<out.size(); i++ ) {
776 out[i] += m[i][j]*vec[j];
783 template <
class X,
class Y>
784 void mv(
const X& x, Y& y)
const
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!");
794 auto&& xj = Impl::asVector(x[j]);
795 auto&& yi = Impl::asVector(y[i]);
796 Impl::asMatrix((*
this)[i][j]).umv(xj, yi);
802 template<
class X,
class Y>
803 void mtv (
const X& x, Y& y)
const
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");
815 template <
class X,
class Y>
816 void umv(
const X& x, Y& y)
const
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!");
825 auto&& xj = Impl::asVector(x[j]);
826 auto&& yi = Impl::asVector(y[i]);
827 Impl::asMatrix((*
this)[i][j]).umv(xj, yi);
832 template<
class X,
class Y>
833 void mmv (
const X& x, Y& y)
const
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!");
842 auto&& xj = Impl::asVector(x[j]);
843 auto&& yi = Impl::asVector(y[i]);
844 Impl::asMatrix((*
this)[i][j]).mmv(xj, yi);
849 template <
class X,
class Y>
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!");
859 auto&& xj = Impl::asVector(x[j]);
860 auto&& yi = Impl::asVector(y[i]);
861 Impl::asMatrix((*
this)[i][j]).usmv(alpha, xj, yi);
866 template<
class X,
class Y>
867 void umtv (
const X& x, Y& y)
const
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!");
876 auto&& xi = Impl::asVector(x[i]);
877 auto&& yj = Impl::asVector(y[j]);
878 Impl::asMatrix((*
this)[i][j]).umtv(xi, yj);
883 template<
class X,
class Y>
884 void mmtv (
const X& x, Y& y)
const
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!");
893 auto&& xi = Impl::asVector(x[i]);
894 auto&& yj = Impl::asVector(y[j]);
895 Impl::asMatrix((*
this)[i][j]).mmtv(xi, yj);
900 template<
class X,
class Y>
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!");
910 auto&& xi = Impl::asVector(x[i]);
911 auto&& yj = Impl::asVector(y[j]);
912 Impl::asMatrix((*
this)[i][j]).usmtv(alpha, xi, yj);
917 template<
class X,
class Y>
918 void umhv (
const X& x, Y& y)
const
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!");
927 auto&& xi = Impl::asVector(x[i]);
928 auto&& yj = Impl::asVector(y[j]);
929 Impl::asMatrix((*
this)[i][j]).umhv(xi,yj);
934 template<
class X,
class Y>
935 void mmhv (
const X& x, Y& y)
const
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!");
944 auto&& xi = Impl::asVector(x[i]);
945 auto&& yj = Impl::asVector(y[j]);
946 Impl::asMatrix((*
this)[i][j]).mmhv(xi,yj);
951 template<
class X,
class Y>
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!");
961 auto&& xi = Impl::asVector(x[i]);
962 auto&& yj = Impl::asVector(y[j]);
963 Impl::asMatrix((*
this)[i][j]).usmhv(alpha,xi,yj);
978 typename FieldTraits<field_type>::real_type sum=0;
981 sum += Impl::asMatrix(
data_[i][j]).frobenius_norm2();
987 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
989 using real_type =
typename FieldTraits<ft>::real_type;
993 for (
auto const &x : *
this) {
995 for (
auto const &y : x)
996 sum += Impl::asMatrix(y).infinity_norm();
997 norm = max(sum, norm);
1006 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
1008 using real_type =
typename FieldTraits<ft>::real_type;
1012 for (
auto const &x : *
this) {
1014 for (
auto const &y : x)
1015 sum += Impl::asMatrix(y).infinity_norm_real();
1016 norm = max(sum, norm);
1023 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
1025 using real_type =
typename FieldTraits<ft>::real_type;
1029 real_type isNaN = 1;
1030 for (
auto const &x : *
this) {
1032 for (
auto const &y : x)
1033 sum += Impl::asMatrix(y).infinity_norm();
1034 norm = max(sum, norm);
1038 return norm * (isNaN / isNaN);
1043 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
1045 using real_type =
typename FieldTraits<ft>::real_type;
1049 real_type isNaN = 1;
1050 for (
auto const &x : *
this) {
1052 for (
auto const &y : x)
1053 sum += Impl::asMatrix(y).infinity_norm_real();
1054 norm = max(sum, norm);
1058 return norm * (isNaN / isNaN);
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");
1070 DUNE_UNUSED_PARAMETER(i); DUNE_UNUSED_PARAMETER(j);
const row_type operator[](size_type row) const
The const index operator.
Definition: matrix.hh:680
A allocator_type
Export the allocator.
Definition: matrix.hh:568
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:867
ConstIterator(const B *data, size_type columns, size_type _i)
constructor from pointer
Definition: matrix.hh:411
Iterator & operator=(Iterator &&other)
Move assignment.
Definition: matrix.hh:279
bool operator!=(const Iterator &it) const
inequality
Definition: matrix.hh:319
ConstIterator end() const
end ConstIterator
Definition: matrix.hh:515
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:583
void umhv(const X &x, Y &y) const
y += A^H x
Definition: matrix.hh:918
Matrix(size_type rows, size_type cols)
Create uninitialized matrix of size rows x cols.
Definition: matrix.hh:597
Iterator class for sequential access.
Definition: matrix.hh:259
size_type M() const
Return the number of columns.
Definition: matrix.hh:696
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:586
static constexpr unsigned int blocklevel
The number of nesting levels the matrix contains.
Definition: matrix.hh:589
RowIterator beforeBegin()
Definition: matrix.hh:630
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:751
const window_type & operator*() const
dereferencing
Definition: matrix.hh:478
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:616
ConstIterator()
constructor
Definition: matrix.hh:404
ConstIterator & operator=(Iterator &&other)
Definition: matrix.hh:421
Iterator beforeEnd()
Definition: matrix.hh:375
Imp::BlockVectorWindow< B, A > window_type
Definition: matrix.hh:67
Iterator & operator--()
prefix decrement
Definition: matrix.hh:305
row_type operator[](size_type row)
The index operator.
Definition: matrix.hh:669
BlockVector< B, A > block_type
Same as value_type, here for historical reasons.
Definition: matrix.hh:64
ConstIterator(const Iterator &it)
constructor from non_const iterator
Definition: matrix.hh:417
reference operator[](size_type i)
random access to blocks
Definition: matrix.hh:238
Matrix()
Create empty matrix.
Definition: matrix.hh:592
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:976
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:884
MatrixImp::DenseMatrixBase< T, A >::Iterator RowIterator
Iterator over the matrix rows.
Definition: matrix.hh:577
DenseMatrixBase(const DenseMatrixBase &a)
copy constructor, has copy semantics
Definition: matrix.hh:113
size_type index() const
Definition: matrix.hh:349
MatrixImp::DenseMatrixBase< T, A > data_
Abuse DenseMatrixBase as an engine for a 2d array ISTL-style.
Definition: matrix.hh:1079
ConstRowIterator begin() const
Get const iterator to first row.
Definition: matrix.hh:636
void umv(const X &x, Y &y) const
y += A x
Definition: matrix.hh:816
Matrix< T > & operator/=(const field_type &scalar)
Division by a scalar.
Definition: matrix.hh:707
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:562
A generic dynamic dense matrix.
Definition: matrix.hh:557
size_type cols_
Number of columns of the matrix.
Definition: matrix.hh:1086
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:970
void mtv(const X &x, Y &y) const
y = A^T x
Definition: matrix.hh:803
window_type reference
Definition: matrix.hh:69
Iterator iterator
Export the iterator type using std naming rules.
Definition: matrix.hh:503
Iterator beforeBegin() const
Definition: matrix.hh:382
ConstIterator beforeEnd() const
Definition: matrix.hh:522
Iterator find(size_type i)
random access returning iterator (end if not contained)
Definition: matrix.hh:388
Matrix< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:701
Iterator()
constructor, no arguments
Definition: matrix.hh:263
DenseMatrixBase & operator=(const DenseMatrixBase &a)
assignment
Definition: matrix.hh:179
A::size_type size_type
The size type for the index access.
Definition: matrix.hh:53
row_type::iterator ColIterator
Iterator for the entries of each row.
Definition: matrix.hh:580
const window_type * operator->() const
arrow
Definition: matrix.hh:484
void setSize(size_type rows, size_type cols)
Change the matrix size.
Definition: matrix.hh:604
size_type N() const
Return the number of rows.
Definition: matrix.hh:691
bool operator==(const Iterator &it) const
equality
Definition: matrix.hh:313
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:1007
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:717
Iterator & operator++()
prefix increment
Definition: matrix.hh:297
Iterator(B *data, size_type columns, size_type _i)
constructor
Definition: matrix.hh:273
Iterator end()
end Iterator
Definition: matrix.hh:368
derive error class from the base class in common
Definition: istlexception.hh:16
BlockVector< B, A > value_type
Type of the elements of the outer vector, i.e., dynamic vectors of B.
Definition: matrix.hh:60
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:833
ConstRowIterator beforeBegin() const
Definition: matrix.hh:656
window_type * operator->() const
arrow
Definition: matrix.hh:343
Definition: allocator.hh:7
friend Y operator*(const Matrix< T > &m, const X &vec)
Generic matrix-vector multiplication.
Definition: matrix.hh:766
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 type
Definition: matrix.hh:50
const typedef window_type const_reference
Definition: matrix.hh:71
void resize(size_type rows, size_type columns)
same effect as constructor with same argument
Definition: matrix.hh:150
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:571
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: matrix.hh:850
~DenseMatrixBase()
free dynamic memory
Definition: matrix.hh:139
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
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:574
RowIterator beforeEnd()
Definition: matrix.hh:623
window_type & operator*() const
dereferencing
Definition: matrix.hh:337
A Vector of blocks with different blocksizes.
Definition: matrix.hh:37
Iterator begin()
begin Iterator
Definition: matrix.hh:362
ConstRowIterator end() const
Get const iterator to one beyond last row.
Definition: matrix.hh:642
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:662
bool operator!=(const ConstIterator &it) const
inequality
Definition: matrix.hh:460
ConstIterator find(size_type i) const
random access returning iterator (end if not contained)
Definition: matrix.hh:394
ConstIterator begin() const
begin ConstIterator
Definition: matrix.hh:509
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: matrix.hh:952
ConstIterator rend() const
end ConstIterator
Definition: matrix.hh:528
ConstIterator const_iterator
Export the const iterator type using std naming rules.
Definition: matrix.hh:506
Matrix transpose() const
Return the transpose of the matrix.
Definition: matrix.hh:741
DenseMatrixBase()
Definition: matrix.hh:79
typename Imp::BlockTraits< T >::field_type field_type
export the type representing the field
Definition: matrix.hh:47
ConstIterator & operator--()
prefix decrement
Definition: matrix.hh:446
Iterator & operator=(Iterator &other)
Copy assignment.
Definition: matrix.hh:288
A vector of blocks with memory management.
Definition: bvector.hh:402
T block_type
Export the type representing the components.
Definition: matrix.hh:565
ConstIterator & operator=(Iterator &other)
Definition: matrix.hh:429
ConstRowIterator beforeEnd() const
Definition: matrix.hh:649
size_type index() const
Definition: matrix.hh:490
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Matrix & operator-=(const Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:731
DenseMatrixBase(size_type rows, size_type columns)
Definition: matrix.hh:92
size_type N() const
number of blocks in the vector (are of variable size here)
Definition: matrix.hh:536
bool operator==(const ConstIterator &it) const
equality
Definition: matrix.hh:454
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:1064
ConstIterator class for sequential access.
Definition: matrix.hh:400
void mv(const X &x, Y &y) const
y = A x
Definition: matrix.hh:784
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: matrix.hh:901
ConstIterator & operator++()
prefix increment
Definition: matrix.hh:438