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>
37 template<
class B,
class A=std::allocator<B> >
48 using field_type =
typename Imp::BlockTraits<B>::field_type;
96 this->n = rows*columns;
100 this->p = allocator_.allocate(this->n);
101 new (this->p)B[this->n];
118 columns_ = a.columns_;
122 this->p = allocator_.allocate(this->n);
123 new (this->p)B[this->n];
146 allocator_.deallocate(this->p,this->n);
158 allocator_.deallocate(this->p,this->n);
162 this->n = rows*columns;
165 this->p = allocator_.allocate(this->n);
166 new (this->p)B[this->n];
184 columns_ = a.columns_;
187 if (this->n!=a.n || rows_!=a.rows_)
194 allocator_.deallocate(this->p,this->n);
202 this->p = allocator_.allocate(this->n);
203 new (this->p)B[this->n];
229 (
static_cast<Imp::block_vector_unmanaged<B,A>&
>(*this)) = k;
241 #ifdef DUNE_ISTL_WITH_CHECKING
242 if (i>=rows_) DUNE_THROW(
ISTLError,
"index out of range");
244 return window_type(this->p + i*columns_, columns_);
250 #ifdef DUNE_ISTL_WITH_CHECKING
251 if (i<0 || i>=rows_) DUNE_THROW(
ISTLError,
"index out of range");
253 return window_type(this->p + i*columns_, columns_);
276 window_(data + _i*columns, columns)
284 window_.set(other.window_.getsize(),other.window_.getptr());
293 window_.set(other.window_.getsize(),other.window_.getptr());
301 window_.setptr(window_.getptr()+window_.getsize());
309 window_.setptr(window_.getptr()-window_.getsize());
316 return window_.getptr() == it.window_.getptr();
322 return window_.getptr() != it.window_.getptr();
328 return window_.getptr() == it.window_.getptr();
334 return window_.getptr() != it.window_.getptr();
365 return Iterator(this->p, columns_, 0);
371 return Iterator(this->p, columns_, rows_);
378 return Iterator(this->p, columns_, rows_-1);
385 return Iterator(this->p, columns_, -1);
391 return Iterator(this->p, columns_, std::min(i,rows_));
414 window_(const_cast<B*>(data + _i * columns), columns)
419 : i(it.i), window_(it.window_.getptr(),it.window_.getsize())
426 window_.set(other.window_.getsize(),other.window_.getptr());
434 window_.set(other.window_.getsize(),other.window_.getptr());
442 window_.setptr(window_.getptr()+window_.getsize());
450 window_.setptr(window_.getptr()-window_.getsize());
457 return window_.getptr() == it.window_.getptr();
463 return window_.getptr() != it.window_.getptr();
469 return window_.getptr() == it.window_.getptr();
475 return window_.getptr() != it.window_.getptr();
557 template<
class T,
class A=std::allocator<T> >
590 [[deprecated(
"Use free function blockLevel(). Will be removed after 2.8.")]]
607 data_.resize(rows,cols);
614 return data_.begin();
627 return data_.beforeEnd();
634 return data_.beforeBegin();
640 return data_.begin();
653 return data_.beforeEnd();
660 return data_.beforeBegin();
672 #ifdef DUNE_ISTL_WITH_CHECKING
674 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
676 DUNE_THROW(
ISTLError,
"Row index out of range!");
683 #ifdef DUNE_ISTL_WITH_CHECKING
685 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
687 DUNE_THROW(
ISTLError,
"Row index out of range!");
720 #ifdef DUNE_ISTL_WITH_CHECKING
721 if(
N()!=b.
N() ||
M() != b.
M())
722 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
734 #ifdef DUNE_ISTL_WITH_CHECKING
735 if(
N()!=b.
N() ||
M() != b.
M())
736 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
747 out[j][i] = (*
this)[i][j];
760 out[i][j] += m1[i][k]*m2[k][j];
767 template <
class X,
class Y>
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!");
776 for (
size_type i=0; i<out.size(); i++ ) {
778 out[i] += m[i][j]*vec[j];
785 template <
class X,
class Y>
786 void mv(
const X& x, Y& y)
const
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!");
796 auto&& xj = Impl::asVector(x[j]);
797 auto&& yi = Impl::asVector(y[i]);
798 Impl::asMatrix((*
this)[i][j]).umv(xj, yi);
804 template<
class X,
class Y>
805 void mtv (
const X& x, Y& y)
const
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");
817 template <
class X,
class Y>
818 void umv(
const X& x, Y& y)
const
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!");
827 auto&& xj = Impl::asVector(x[j]);
828 auto&& yi = Impl::asVector(y[i]);
829 Impl::asMatrix((*
this)[i][j]).umv(xj, yi);
834 template<
class X,
class Y>
835 void mmv (
const X& x, Y& y)
const
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!");
844 auto&& xj = Impl::asVector(x[j]);
845 auto&& yi = Impl::asVector(y[i]);
846 Impl::asMatrix((*
this)[i][j]).mmv(xj, yi);
851 template <
class X,
class Y>
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!");
861 auto&& xj = Impl::asVector(x[j]);
862 auto&& yi = Impl::asVector(y[i]);
863 Impl::asMatrix((*
this)[i][j]).usmv(alpha, xj, yi);
868 template<
class X,
class Y>
869 void umtv (
const X& x, Y& y)
const
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!");
878 auto&& xi = Impl::asVector(x[i]);
879 auto&& yj = Impl::asVector(y[j]);
880 Impl::asMatrix((*
this)[i][j]).umtv(xi, yj);
885 template<
class X,
class Y>
886 void mmtv (
const X& x, Y& y)
const
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!");
895 auto&& xi = Impl::asVector(x[i]);
896 auto&& yj = Impl::asVector(y[j]);
897 Impl::asMatrix((*
this)[i][j]).mmtv(xi, yj);
902 template<
class X,
class Y>
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!");
912 auto&& xi = Impl::asVector(x[i]);
913 auto&& yj = Impl::asVector(y[j]);
914 Impl::asMatrix((*
this)[i][j]).usmtv(alpha, xi, yj);
919 template<
class X,
class Y>
920 void umhv (
const X& x, Y& y)
const
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!");
929 auto&& xi = Impl::asVector(x[i]);
930 auto&& yj = Impl::asVector(y[j]);
931 Impl::asMatrix((*
this)[i][j]).umhv(xi,yj);
936 template<
class X,
class Y>
937 void mmhv (
const X& x, Y& y)
const
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!");
946 auto&& xi = Impl::asVector(x[i]);
947 auto&& yj = Impl::asVector(y[j]);
948 Impl::asMatrix((*
this)[i][j]).mmhv(xi,yj);
953 template<
class X,
class Y>
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!");
963 auto&& xi = Impl::asVector(x[i]);
964 auto&& yj = Impl::asVector(y[j]);
965 Impl::asMatrix((*
this)[i][j]).usmhv(alpha,xi,yj);
980 typename FieldTraits<field_type>::real_type sum=0;
983 sum += Impl::asMatrix(
data_[i][j]).frobenius_norm2();
989 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
991 using real_type =
typename FieldTraits<ft>::real_type;
995 for (
auto const &x : *
this) {
997 for (
auto const &y : x)
998 sum += Impl::asMatrix(y).infinity_norm();
999 norm = max(sum, norm);
1008 typename std::enable_if<!HasNaN<ft>::value,
int>::type = 0>
1010 using real_type =
typename FieldTraits<ft>::real_type;
1014 for (
auto const &x : *
this) {
1016 for (
auto const &y : x)
1017 sum += Impl::asMatrix(y).infinity_norm_real();
1018 norm = max(sum, norm);
1025 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
1027 using real_type =
typename FieldTraits<ft>::real_type;
1031 real_type isNaN = 1;
1032 for (
auto const &x : *
this) {
1034 for (
auto const &y : x)
1035 sum += Impl::asMatrix(y).infinity_norm();
1036 norm = max(sum, norm);
1040 return norm * (isNaN / isNaN);
1045 typename std::enable_if<HasNaN<ft>::value,
int>::type = 0>
1047 using real_type =
typename FieldTraits<ft>::real_type;
1051 real_type isNaN = 1;
1052 for (
auto const &x : *
this) {
1054 for (
auto const &y : x)
1055 sum += Impl::asMatrix(y).infinity_norm_real();
1056 norm = max(sum, norm);
1060 return norm * (isNaN / isNaN);
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");
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...
Definition: allocator.hh:9
A vector of blocks with memory management.
Definition: bvector.hh:393
derive error class from the base class in common
Definition: istlexception.hh:17
A Vector of blocks with different blocksizes.
Definition: matrix.hh:42
Imp::BlockVectorWindow< B, A > window_type
Definition: matrix.hh:68
BlockVector< B, A > block_type
Same as value_type, here for historical reasons.
Definition: matrix.hh:65
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
const window_type const_reference
Definition: matrix.hh:72
DenseMatrixBase(size_type rows, size_type columns)
Definition: matrix.hh:93
DenseMatrixBase & operator=(const DenseMatrixBase &a)
assignment
Definition: matrix.hh:180
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
window_type reference
Definition: matrix.hh:70
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
Iterator class for sequential access.
Definition: matrix.hh:261
size_type index() const
Definition: matrix.hh:350
Iterator(Iterator &other)=default
Iterator(Iterator &&other)=default
Iterator & operator++()
prefix increment
Definition: matrix.hh:298
bool operator!=(const Iterator &it) const
inequality
Definition: matrix.hh:320
Iterator & operator--()
prefix decrement
Definition: matrix.hh:306
window_type & operator*() const
dereferencing
Definition: matrix.hh:338
Iterator & operator=(Iterator &&other)
Move assignment.
Definition: matrix.hh:280
Iterator()
constructor, no arguments
Definition: matrix.hh:264
window_type * operator->() const
arrow
Definition: matrix.hh:344
bool operator==(const Iterator &it) const
equality
Definition: matrix.hh:314
Iterator(B *data, size_type columns, size_type _i)
constructor
Definition: matrix.hh:274
Iterator & operator=(Iterator &other)
Copy assignment.
Definition: matrix.hh:289
ConstIterator class for sequential access.
Definition: matrix.hh:402
const window_type & operator*() const
dereferencing
Definition: matrix.hh:479
ConstIterator & operator++()
prefix increment
Definition: matrix.hh:439
ConstIterator & operator--()
prefix decrement
Definition: matrix.hh:447
ConstIterator(const B *data, size_type columns, size_type _i)
constructor from pointer
Definition: matrix.hh:412
ConstIterator(const Iterator &it)
constructor from non_const iterator
Definition: matrix.hh:418
ConstIterator & operator=(Iterator &other)
Definition: matrix.hh:430
bool operator!=(const ConstIterator &it) const
inequality
Definition: matrix.hh:461
ConstIterator & operator=(Iterator &&other)
Definition: matrix.hh:422
ConstIterator()
constructor
Definition: matrix.hh:405
bool operator==(const ConstIterator &it) const
equality
Definition: matrix.hh:455
size_type index() const
Definition: matrix.hh:491
const window_type * operator->() const
arrow
Definition: matrix.hh:485
A generic dynamic dense matrix.
Definition: matrix.hh:559
size_type cols_
Number of columns of the matrix.
Definition: matrix.hh:1086
A allocator_type
Export the allocator.
Definition: matrix.hh:569
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:972
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
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:664
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
Matrix< T > & operator/=(const field_type &scalar)
Division by a scalar.
Definition: matrix.hh:709
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
bool exists([[maybe_unused]] size_type i, [[maybe_unused]] size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:1066
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
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:990
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 Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:733
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
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
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:1009
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
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:753
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
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:869
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:978
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
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:719
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< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:703
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