3 #ifndef DUNE_ISTL_MATRIX_HH
4 #define DUNE_ISTL_MATRIX_HH
14 #include <dune/common/ftraits.hh>
23 template<
class T,
class A=std::allocator<T> >
73 void setSize(size_type rows, size_type cols) {
111 ConstRowIterator
end()
const
139 #ifdef DUNE_ISTL_WITH_CHECKING
141 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
143 DUNE_THROW(
ISTLError,
"Row index out of range!");
150 #ifdef DUNE_ISTL_WITH_CHECKING
152 DUNE_THROW(
ISTLError,
"Can't access negative rows!");
154 DUNE_THROW(
ISTLError,
"Row index out of range!");
160 size_type
N()
const {
165 size_type
M()
const {
171 #ifdef DUNE_ISTL_WITH_CHECKING
173 DUNE_THROW(
ISTLError,
"Can't compute rowdim() when there are no columns!");
176 for (size_type i=0; i<
data_.
N(); i++)
184 #ifdef DUNE_ISTL_WITH_CHECKING
186 DUNE_THROW(
ISTLError,
"Can't compute coldim() when there are no rows!");
189 for (size_type i=0; i<
data_[0].
size(); i++)
197 #ifdef DUNE_ISTL_WITH_CHECKING
199 DUNE_THROW(
ISTLError,
"Rowdim for nonexisting row " << r <<
" requested!");
201 DUNE_THROW(
ISTLError,
"Can't compute rowdim() when there are no columns!");
203 return data_[r][0].rowdim();
208 #ifdef DUNE_ISTL_WITH_CHECKING
210 DUNE_THROW(
ISTLError,
"Coldim for nonexisting column " << c <<
" requested!");
212 DUNE_THROW(
ISTLError,
"Can't compute coldim() when there are no rows!");
214 return data_[0][c].coldim();
235 #ifdef DUNE_ISTL_WITH_CHECKING
236 if(
N()!=b.
N() ||
M() != b.
M())
237 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
249 #ifdef DUNE_ISTL_WITH_CHECKING
250 if(
N()!=b.
N() ||
M() != b.
M())
251 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
260 for (size_type i=0; i<
M(); i++)
261 for (size_type j=0; j<
N(); j++)
262 out[j][i] = (*
this)[i][j];
272 for (size_type i=0; i<out.
N(); i++ ) {
273 for ( size_type j=0; j<out.
M(); j++ )
274 for (size_type k=0; k<m1.
M(); k++)
275 out[i][j] += m1[i][k]*m2[k][j];
282 template <
class X,
class Y>
284 #ifdef DUNE_ISTL_WITH_CHECKING
285 if (m.
M()!=vec.size())
286 DUNE_THROW(
ISTLError,
"Vector size doesn't match the number of matrix columns!");
291 for (size_type i=0; i<out.size(); i++ ) {
292 for ( size_type j=0; j<vec.size(); j++ )
293 out[i] += m[i][j]*vec[j];
300 template <
class X,
class Y>
301 void mv(
const X& x, Y& y)
const
303 #ifdef DUNE_ISTL_WITH_CHECKING
304 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
305 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
308 for (size_type i=0; i<
data_.
N(); i++) {
310 for (size_type j=0; j<
cols_; j++)
311 (*
this)[i][j].umv(x[j], y[i]);
318 template<
class X,
class Y>
319 void mtv (
const X& x, Y& y)
const
321 #ifdef DUNE_ISTL_WITH_CHECKING
322 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
323 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
326 for(size_type i=0; i<y.N(); ++i)
332 template <
class X,
class Y>
333 void umv(
const X& x, Y& y)
const
335 #ifdef DUNE_ISTL_WITH_CHECKING
336 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
337 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
340 for (size_type i=0; i<
data_.
N(); i++) {
342 for (size_type j=0; j<
cols_; j++)
343 (*
this)[i][j].umv(x[j], y[i]);
350 template<
class X,
class Y>
351 void mmv (
const X& x, Y& y)
const
353 #ifdef DUNE_ISTL_WITH_CHECKING
354 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
355 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
357 ConstRowIterator endi=
end();
358 for (ConstRowIterator i=
begin(); i!=endi; ++i)
360 ConstColIterator endj = (*i).end();
361 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
362 (*j).mmv(x[j.index()],y[i.index()]);
367 template <
class X,
class Y>
368 void usmv(
const field_type& alpha,
const X& x, Y& y)
const
370 #ifdef DUNE_ISTL_WITH_CHECKING
371 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
372 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
375 for (size_type i=0; i<
data_.
N(); i++) {
377 for (size_type j=0; j<
cols_; j++)
378 (*
this)[i][j].usmv(alpha, x[j], y[i]);
385 template<
class X,
class Y>
386 void umtv (
const X& x, Y& y)
const
388 #ifdef DUNE_ISTL_WITH_CHECKING
389 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
390 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
392 ConstRowIterator endi=
end();
393 for (ConstRowIterator i=
begin(); i!=endi; ++i)
395 ConstColIterator endj = (*i).end();
396 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
397 (*j).umtv(x[i.index()],y[j.index()]);
402 template<
class X,
class Y>
403 void mmtv (
const X& x, Y& y)
const
405 #ifdef DUNE_ISTL_WITH_CHECKING
406 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
407 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
409 ConstRowIterator endi=
end();
410 for (ConstRowIterator i=
begin(); i!=endi; ++i)
412 ConstColIterator endj = (*i).end();
413 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
414 (*j).mmtv(x[i.index()],y[j.index()]);
419 template<
class X,
class Y>
420 void usmtv (
const field_type& alpha,
const X& x, Y& y)
const
422 #ifdef DUNE_ISTL_WITH_CHECKING
423 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
424 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
426 ConstRowIterator endi=
end();
427 for (ConstRowIterator i=
begin(); i!=endi; ++i)
429 ConstColIterator endj = (*i).end();
430 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
431 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
436 template<
class X,
class Y>
437 void umhv (
const X& x, Y& y)
const
439 #ifdef DUNE_ISTL_WITH_CHECKING
440 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
441 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
443 ConstRowIterator endi=
end();
444 for (ConstRowIterator i=
begin(); i!=endi; ++i)
446 ConstColIterator endj = (*i).end();
447 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
448 (*j).umhv(x[i.index()],y[j.index()]);
453 template<
class X,
class Y>
454 void mmhv (
const X& x, Y& y)
const
456 #ifdef DUNE_ISTL_WITH_CHECKING
457 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
458 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
460 ConstRowIterator endi=
end();
461 for (ConstRowIterator i=
begin(); i!=endi; ++i)
463 ConstColIterator endj = (*i).end();
464 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
465 (*j).mmhv(x[i.index()],y[j.index()]);
470 template<
class X,
class Y>
471 void usmhv (
const field_type& alpha,
const X& x, Y& y)
const
473 #ifdef DUNE_ISTL_WITH_CHECKING
474 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
475 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"vector/matrix size mismatch!");
477 ConstRowIterator endi=
end();
478 for (ConstRowIterator i=
begin(); i!=endi; ++i)
480 ConstColIterator endj = (*i).end();
481 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
482 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
498 for (size_type i=0; i<
N(); ++i)
499 for (size_type j=0; j<
M(); ++j)
508 for (size_type i=0; i<
N(); ++i) {
510 for (size_type j=0; j<
M(); j++)
512 max = std::max(max,sum);
521 for (size_type i=0; i<
N(); ++i) {
523 for (size_type j=0; j<
M(); j++)
525 max = std::max(max,sum);
533 bool exists (size_type i, size_type j)
const
535 #ifdef DUNE_ISTL_WITH_CHECKING
536 if (i<0 || i>=
N()) DUNE_THROW(
ISTLError,
"row index out of range");
537 if (j<0 || i>=
M()) DUNE_THROW(
ISTLError,
"column index out of range");
539 DUNE_UNUSED_PARAMETER(i); DUNE_UNUSED_PARAMETER(j);
Row row
Definition: matrixmatrix.hh:345
VariableBlockVector< T, A > data_
Abuse VariableBlockVector as an engine for a 2d array ISTL-style.
Definition: matrix.hh:551
ConstRowIterator begin() const
Get const iterator to first row.
Definition: matrix.hh:105
size_type coldim(size_type c) const
The number of scalar columns.
Definition: matrix.hh:207
VariableBlockVector< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:38
derive error class from the base class in common
Definition: istlexception.hh:16
Matrix< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:218
void mv(const X &x, Y &y) const
y = A x
Definition: matrix.hh:301
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
Iterator begin()
begin Iterator
Definition: vbvector.hh:602
void umhv(const X &x, Y &y) const
y += A^H x
Definition: matrix.hh:437
Matrix & A
Definition: matrixmatrix.hh:216
Matrix< T > & operator/=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:224
Matrix transpose() const
Return the transpose of the matrix.
Definition: matrix.hh:258
ConstRowIterator end() const
Get const iterator to one beyond last row.
Definition: matrix.hh:111
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:351
A allocator_type
Export the allocator.
Definition: matrix.hh:35
VariableBlockVector< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:50
size_type coldim() const
The number of scalar columns.
Definition: matrix.hh:183
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:234
Iterator end()
end Iterator
Definition: vbvector.hh:608
Iterator beforeBegin() const
Definition: vbvector.hh:622
VariableBlockVector< T, A >::Iterator RowIterator
Iterator over the matrix rows.
Definition: matrix.hh:44
friend Y operator*(const Matrix< T > &m, const X &vec)
Generic matrix-vector multiplication.
Definition: matrix.hh:283
void mtv(const X &x, Y &y) const
y = A^T x
Definition: matrix.hh:319
void umv(const X &x, Y &y) const
y += A x
Definition: matrix.hh:333
The number of nesting levels the matrix contains.
Definition: matrix.hh:57
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
void resize(size_type _nblocks)
same effect as constructor with same argument
Definition: vbvector.hh:214
size_type N() const
Return the number of rows.
Definition: matrix.hh:160
A Vector of blocks with different blocksizes.
Definition: vbvector.hh:36
ConstRowIterator beforeBegin() const
Definition: matrix.hh:125
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: matrix.hh:420
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:85
Iterator beforeEnd()
Definition: vbvector.hh:615
void setSize(size_type rows, size_type cols)
Change the matrix size.
Definition: matrix.hh:73
RowIterator beforeBegin()
Definition: matrix.hh:99
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:495
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: matrix.hh:454
const row_type & operator[](size_type row) const
The const index operator.
Definition: matrix.hh:149
size_type cols_
Number of columns of the matrix.
Definition: matrix.hh:558
size_type rowdim(size_type r) const
The number of scalar rows.
Definition: matrix.hh:196
Matrix(size_type rows, size_type cols)
Create uninitialized matrix of size rows x cols.
Definition: matrix.hh:66
Matrix & operator-=(const Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:248
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:518
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:505
ConstRowIterator beforeEnd() const
Definition: matrix.hh:118
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:29
Matrix()
Create empty matrix.
Definition: matrix.hh:61
A generic dynamic dense matrix.
Definition: matrix.hh:24
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: matrix.hh:368
size_type rowdim() const
The number of scalar rows.
Definition: matrix.hh:170
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:41
row_type::iterator ColIterator
Iterator for the entries of each row.
Definition: matrix.hh:47
size_type N() const
number of blocks in the vector (are of variable size here)
Definition: vbvector.hh:756
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:131
row_type & operator[](size_type row)
The index operator.
Definition: matrix.hh:138
T block_type
Export the type representing the components.
Definition: matrix.hh:32
Definition: basearray.hh:19
size_type size() const
number of blocks in the array (are of size 1 here)
Definition: basearray.hh:236
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:386
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:268
RowIterator beforeEnd()
Definition: matrix.hh:92
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:533
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:403
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: matrix.hh:471
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:489