3 #ifndef DUNE_ISTL_MATRIX_HH 4 #define DUNE_ISTL_MATRIX_HH 13 #include <dune/common/ftraits.hh> 33 template<
class B,
class A=std::allocator<B> >
92 this->
n = rows*columns;
96 this->
p = allocator_.allocate(this->
n);
97 new (this->
p)B[this->
n];
114 columns_ = a.columns_;
118 this->
p = allocator_.allocate(this->
n);
119 new (this->
p)B[this->
n];
122 for (size_type i=0; i<this->
n; i++)
142 allocator_.deallocate(this->
p,this->
n);
147 void resize (size_type rows, size_type columns)
154 allocator_.deallocate(this->
p,this->
n);
158 this->
n = rows*columns;
161 this->
p = allocator_.allocate(this->
n);
162 new (this->
p)B[this->
n];
180 columns_ = a.columns_;
183 if (this->
n!=a.
n || rows_!=a.rows_)
190 allocator_.deallocate(this->
p,this->
n);
198 this->
p = allocator_.allocate(this->
n);
199 new (this->
p)B[this->
n];
212 for (size_type i=0; i<this->
n; i++)
237 #ifdef DUNE_ISTL_WITH_CHECKING 238 if (i>=rows_) DUNE_THROW(
ISTLError,
"index out of range");
246 #ifdef DUNE_ISTL_WITH_CHECKING 247 if (i<0 || i>=rows_) DUNE_THROW(
ISTLError,
"index out of range");
270 Iterator (B* data, size_type columns, size_type _i)
272 window_(data + _i*columns, columns)
355 mutable window_type window_;
367 return Iterator(this->
p, columns_, rows_);
374 return Iterator(this->
p, columns_, rows_-1);
387 return Iterator(this->
p, columns_, std::min(i,rows_));
410 window_(const_cast<B*>(data + _i * columns), columns)
415 : i(it.i), window_(it.window_.getptr(),it.window_.getsize())
422 window_.set(other.window_.getsize(),other.window_.getptr());
438 window_.setptr(window_.getptr()+window_.getsize());
446 window_.setptr(window_.getptr()-window_.getsize());
453 return window_.getptr() == it.window_.
getptr();
459 return window_.getptr() != it.window_.
getptr();
465 return window_.getptr() == it.window_.
getptr();
471 return window_.getptr() != it.window_.
getptr();
496 mutable window_type window_;
553 template<
class T,
class A=std::allocator<T> >
587 blocklevel = T::blocklevel+1
596 Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
603 void setSize(size_type rows, size_type cols) {
604 data_.resize(rows,cols);
611 return data_.begin();
624 return data_.beforeEnd();
631 return data_.beforeBegin();
637 return data_.begin();
641 ConstRowIterator
end()
const 650 return data_.beforeEnd();
657 return data_.beforeBegin();
669 #ifdef DUNE_ISTL_WITH_CHECKING 671 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
673 DUNE_THROW(
ISTLError,
"Row index out of range!");
680 #ifdef DUNE_ISTL_WITH_CHECKING 682 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
684 DUNE_THROW(
ISTLError,
"Row index out of range!");
690 size_type
N()
const {
695 size_type
M()
const {
717 #ifdef DUNE_ISTL_WITH_CHECKING 718 if(
N()!=b.
N() || M() != b.
M())
719 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
731 #ifdef DUNE_ISTL_WITH_CHECKING 732 if(
N()!=b.
N() || M() != b.
M())
733 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
742 for (size_type i=0; i<
N(); i++)
743 for (size_type j=0; j<M(); j++)
744 out[j][i] = (*
this)[i][j];
754 for (size_type i=0; i<out.
N(); i++ ) {
755 for ( size_type j=0; j<out.
M(); j++ )
756 for (size_type k=0; k<m1.
M(); k++)
757 out[i][j] += m1[i][k]*m2[k][j];
764 template <
class X,
class Y>
766 #ifdef DUNE_ISTL_WITH_CHECKING 767 if (m.
M()!=vec.size())
768 DUNE_THROW(
ISTLError,
"Vector size doesn't match the number of matrix columns!");
773 for (size_type i=0; i<out.size(); i++ ) {
774 for ( size_type j=0; j<vec.size(); j++ )
775 out[i] += m[i][j]*vec[j];
782 template <
class X,
class Y>
783 void mv(
const X& x, Y& y)
const 785 #ifdef DUNE_ISTL_WITH_CHECKING 786 if (x.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
787 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
790 for (size_type i=0; i<data_.N(); i++) {
792 for (size_type j=0; j<cols_; j++)
793 (*
this)[i][j].umv(x[j], y[i]);
800 template<
class X,
class Y>
801 void mtv (
const X& x, Y& y)
const 803 #ifdef DUNE_ISTL_WITH_CHECKING 804 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
805 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"index out of range");
808 for(size_type i=0; i<y.N(); ++i)
814 template <
class X,
class Y>
815 void umv(
const X& x, Y& y)
const 817 #ifdef DUNE_ISTL_WITH_CHECKING 818 if (x.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
819 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
822 for (size_type i=0; i<data_.N(); i++) {
824 for (size_type j=0; j<cols_; j++)
825 (*
this)[i][j].umv(x[j], y[i]);
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!");
839 ConstRowIterator endi=
end();
840 for (ConstRowIterator i=
begin(); i!=endi; ++i)
842 ConstColIterator endj = (*i).end();
843 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
844 (*j).mmv(x[j.index()],y[i.index()]);
849 template <
class X,
class Y>
850 void usmv(
const field_type& alpha,
const X& x, Y& y)
const 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!");
857 for (size_type i=0; i<data_.N(); i++) {
859 for (size_type j=0; j<cols_; j++)
860 (*
this)[i][j].usmv(alpha, x[j], y[i]);
867 template<
class X,
class Y>
868 void umtv (
const X& x, Y& y)
const 870 #ifdef DUNE_ISTL_WITH_CHECKING 871 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
872 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
874 ConstRowIterator endi=
end();
875 for (ConstRowIterator i=
begin(); i!=endi; ++i)
877 ConstColIterator endj = (*i).end();
878 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
879 (*j).umtv(x[i.index()],y[j.index()]);
884 template<
class X,
class Y>
885 void mmtv (
const X& x, Y& y)
const 887 #ifdef DUNE_ISTL_WITH_CHECKING 888 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
889 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
891 ConstRowIterator endi=
end();
892 for (ConstRowIterator i=
begin(); i!=endi; ++i)
894 ConstColIterator endj = (*i).end();
895 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
896 (*j).mmtv(x[i.index()],y[j.index()]);
901 template<
class X,
class Y>
902 void usmtv (
const field_type& alpha,
const X& x, Y& y)
const 904 #ifdef DUNE_ISTL_WITH_CHECKING 905 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
906 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
908 ConstRowIterator endi=
end();
909 for (ConstRowIterator i=
begin(); i!=endi; ++i)
911 ConstColIterator endj = (*i).end();
912 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
913 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
918 template<
class X,
class Y>
919 void umhv (
const X& x, Y& y)
const 921 #ifdef DUNE_ISTL_WITH_CHECKING 922 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
923 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
925 ConstRowIterator endi=
end();
926 for (ConstRowIterator i=
begin(); i!=endi; ++i)
928 ConstColIterator endj = (*i).end();
929 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
930 (*j).umhv(x[i.index()],y[j.index()]);
935 template<
class X,
class Y>
936 void mmhv (
const X& x, Y& y)
const 938 #ifdef DUNE_ISTL_WITH_CHECKING 939 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
940 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
942 ConstRowIterator endi=
end();
943 for (ConstRowIterator i=
begin(); i!=endi; ++i)
945 ConstColIterator endj = (*i).end();
946 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
947 (*j).mmhv(x[i.index()],y[j.index()]);
952 template<
class X,
class Y>
953 void usmhv (
const field_type& alpha,
const X& x, Y& y)
const 955 #ifdef DUNE_ISTL_WITH_CHECKING 956 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
957 if (y.N()!=M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
959 ConstRowIterator endi=
end();
960 for (ConstRowIterator i=
begin(); i!=endi; ++i)
962 ConstColIterator endj = (*i).end();
963 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
964 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
973 return std::sqrt(frobenius_norm2());
980 for (size_type i=0; i<
N(); ++i)
981 for (size_type j=0; j<M(); ++j)
982 sum += data_[i][j].frobenius_norm2();
988 typename std::enable_if<!has_nan<ft>::value,
int>::type = 0>
990 using real_type =
typename FieldTraits<ft>::real_type;
994 for (
auto const &x : *
this) {
996 for (
auto const &y : x)
997 sum += y.infinity_norm();
998 norm = max(sum, norm);
1005 typename std::enable_if<!has_nan<ft>::value,
int>::type = 0>
1007 using real_type =
typename FieldTraits<ft>::real_type;
1011 for (
auto const &x : *
this) {
1013 for (
auto const &y : x)
1014 sum += y.infinity_norm_real();
1015 norm = max(sum, norm);
1022 typename std::enable_if<has_nan<ft>::value,
int>::type = 0>
1024 using real_type =
typename FieldTraits<ft>::real_type;
1028 real_type isNaN = 1;
1029 for (
auto const &x : *
this) {
1031 for (
auto const &y : x)
1032 sum += y.infinity_norm();
1033 norm = max(sum, norm);
1037 return norm * isNaN;
1042 typename std::enable_if<has_nan<ft>::value,
int>::type = 0>
1044 using real_type =
typename FieldTraits<ft>::real_type;
1048 real_type isNaN = 1;
1049 for (
auto const &x : *
this) {
1051 for (
auto const &y : x)
1052 sum += y.infinity_norm_real();
1053 norm = max(sum, norm);
1057 return norm * isNaN;
1065 #ifdef DUNE_ISTL_WITH_CHECKING 1066 if (i<0 || i>=
N()) DUNE_THROW(
ISTLError,
"row index out of range");
1067 if (j<0 || i>=M()) DUNE_THROW(
ISTLError,
"column index out of range");
1069 DUNE_UNUSED_PARAMETER(i); DUNE_UNUSED_PARAMETER(j);
const window_type const_reference
Definition: matrix.hh:68
Iterator end()
end Iterator
Definition: matrix.hh:365
Matrix(size_type rows, size_type cols)
Create uninitialized matrix of size rows x cols.
Definition: matrix.hh:596
DenseMatrixBase & operator=(const DenseMatrixBase &a)
assignment
Definition: matrix.hh:176
Iterator & operator=(Iterator &other)
Copy assignment.
Definition: matrix.hh:285
block_vector_unmanaged & operator-=(const block_vector_unmanaged &y)
vector space subtraction
Definition: bvector.hh:100
void umhv(const X &x, Y &y) const
y += A^H x
Definition: matrix.hh:919
friend class ConstIterator
Definition: matrix.hh:351
window_type reference
Definition: matrix.hh:66
ConstIterator beforeEnd() const
Definition: matrix.hh:519
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:971
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:571
A generic dynamic dense matrix.
Definition: matrix.hh:554
A vector of blocks with memory management.
Definition: bvector.hh:312
ConstRowIterator beforeBegin() const
Definition: matrix.hh:655
bool operator==(const Iterator &it) const
equality
Definition: matrix.hh:310
ConstIterator & operator--()
prefix decrement
Definition: matrix.hh:443
A Vector of blocks with different blocksizes.
Definition: matrix.hh:34
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
Matrix()
Create empty matrix.
Definition: matrix.hh:591
size_type n
Definition: basearray.hh:256
ConstIterator()
constructor
Definition: matrix.hh:401
size_type getsize()
get size
Definition: bvector.hh:751
void set(size_type _n, B *_p)
set size and pointer
Definition: bvector.hh:726
size_type index() const
Definition: matrix.hh:346
Iterator implementation class.
Definition: basearray.hh:84
const row_type operator[](size_type row) const
The const index operator.
Definition: matrix.hh:679
ConstIterator(const B *data, size_type columns, size_type _i)
constructor from pointer
Definition: matrix.hh:408
A::size_type size_type
The size type for the index access.
Definition: matrix.hh:50
MatrixImp::DenseMatrixBase< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:568
B::field_type field_type
export the type representing the field
Definition: matrix.hh:44
T block_type
Export the type representing the components.
Definition: matrix.hh:562
void umv(const X &x, Y &y) const
y += A x
Definition: matrix.hh:815
RowIterator beforeBegin()
Definition: matrix.hh:629
B * getptr()
get pointer
Definition: bvector.hh:745
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:1063
ConstIterator end() const
end ConstIterator
Definition: matrix.hh:512
DenseMatrixBase()
Definition: matrix.hh:76
void setSize(size_type rows, size_type cols)
Change the matrix size.
Definition: matrix.hh:603
Matrix< T > & operator/=(const field_type &scalar)
Division by a scalar.
Definition: matrix.hh:706
size_type M() const
Return the number of columns.
Definition: matrix.hh:695
Iterator beforeBegin() const
Definition: matrix.hh:379
row_type operator[](size_type row)
The index operator.
Definition: matrix.hh:668
size_type N() const
Return the number of rows.
Definition: matrix.hh:690
RealIterator< const T > const_iterator
iterator class for sequential access
Definition: basearray.hh:204
base_array_unmanaged< T, A >::iterator Iterator
make iterators available as types
Definition: bvector.hh:64
ConstRowIterator begin() const
Get const iterator to first row.
Definition: matrix.hh:635
block_vector_unmanaged & operator+=(const block_vector_unmanaged &y)
vector space addition
Definition: bvector.hh:90
base_array_unmanaged< T, A >::const_iterator ConstIterator
make iterators available as types
Definition: bvector.hh:67
size_type N() const
number of blocks in the vector (are of variable size here)
Definition: matrix.hh:533
void mtv(const X &x, Y &y) const
y = A^T x
Definition: matrix.hh:801
RowIterator beforeEnd()
Definition: matrix.hh:622
size_type index() const
Definition: matrix.hh:487
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:977
MatrixImp::DenseMatrixBase< T, A >::Iterator RowIterator
Iterator over the matrix rows.
Definition: matrix.hh:574
Iterator beforeEnd()
Definition: matrix.hh:372
Iterator find(size_type i)
random access returning iterator (end if not contained)
Definition: matrix.hh:385
ConstIterator find(size_type i) const
random access returning iterator (end if not contained)
Definition: matrix.hh:391
Definition: basearray.hh:19
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: matrix.hh:902
ConstIterator & operator=(Iterator &other)
Definition: matrix.hh:426
BlockVector< B, A > block_type
Same as value_type, here for historical reasons.
Definition: matrix.hh:61
void setptr(B *_p)
set pointer only
Definition: bvector.hh:739
DenseMatrixBase(size_type rows, size_type columns)
Definition: matrix.hh:89
This file implements a vector space as a tensor product of a given vector space. The number of compon...
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:989
Matrix< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:700
ConstRowIterator beforeEnd() const
Definition: matrix.hh:648
void mv(const X &x, Y &y) const
y = A x
Definition: matrix.hh:783
Iterator & operator=(Iterator &&other)
Move assignment.
Definition: matrix.hh:276
B * p
Definition: basearray.hh:257
ConstIterator(const Iterator &it)
constructor from non_const iterator
Definition: matrix.hh:414
Iterator & operator--()
prefix decrement
Definition: matrix.hh:302
Matrix transpose() const
Return the transpose of the matrix.
Definition: matrix.hh:740
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:868
ConstIterator & operator++()
prefix increment
Definition: matrix.hh:435
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:580
DenseMatrixBase(const DenseMatrixBase &a)
copy constructor, has copy semantics
Definition: matrix.hh:110
ConstIterator & operator=(Iterator &&other)
Definition: matrix.hh:418
ConstIterator class for sequential access.
Definition: matrix.hh:397
An unmanaged vector of blocks.
Definition: bvector.hh:45
ConstIterator begin() const
begin ConstIterator
Definition: matrix.hh:506
void resize(size_type rows, size_type columns)
same effect as constructor with same argument
Definition: matrix.hh:147
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:1006
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
Iterator(B *data, size_type columns, size_type _i)
constructor
Definition: matrix.hh:270
friend Y operator*(const Matrix< T > &m, const X &vec)
Generic matrix-vector multiplication.
Definition: matrix.hh:765
Iterator begin()
begin Iterator
Definition: matrix.hh:359
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:885
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: matrix.hh:953
bool operator!=(const Iterator &it) const
inequality
Definition: matrix.hh:316
Iterator & operator++()
prefix increment
Definition: matrix.hh:294
size_type cols_
Number of columns of the matrix.
Definition: matrix.hh:1085
row_type::iterator ColIterator
Iterator for the entries of each row.
Definition: matrix.hh:577
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:750
~DenseMatrixBase()
free dynamic memory
Definition: matrix.hh:136
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: matrix.hh:936
Iterator()
constructor, no arguments
Definition: matrix.hh:260
ConstRowIterator end() const
Get const iterator to one beyond last row.
Definition: matrix.hh:641
RealIterator< T > iterator
iterator type for sequential access
Definition: basearray.hh:168
window_type & operator*() const
dereferencing
Definition: matrix.hh:334
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:833
MatrixImp::DenseMatrixBase< T, A > data_
Abuse DenseMatrixBase as an engine for a 2d array ISTL-style.
Definition: matrix.hh:1078
derive error class from the base class in common
Definition: istlexception.hh:16
reference operator[](size_type i)
random access to blocks
Definition: matrix.hh:235
window_type * operator->() const
arrow
Definition: matrix.hh:340
ConstIterator rend() const
end ConstIterator
Definition: matrix.hh:525
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
A allocator_type
export the allocator type
Definition: matrix.hh:47
BlockVector< B, A > value_type
Type of the elements of the outer vector, i.e., dynamic vectors of B.
Definition: matrix.hh:57
BlockVectorWindow< B, A > window_type
Definition: matrix.hh:64
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
Iterator class for sequential access.
Definition: matrix.hh:256
Definition: bvector.hh:29