Dune Core Modules (2.4.2)

matrix.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_MATRIX_HH
4#define DUNE_ISTL_MATRIX_HH
5
10#include <vector>
11#include <memory>
12
13#include <dune/istl/vbvector.hh>
15
16namespace Dune {
17
23 template<class T, class A=std::allocator<T> >
24 class Matrix
25 {
26 public:
27
29 typedef typename T::field_type field_type;
30
32 typedef T block_type;
33
35 typedef A allocator_type;
36
39
41 typedef typename A::size_type size_type;
42
45
48
51
54
55 enum {
57 blocklevel = T::blocklevel+1
58 };
59
61 Matrix() : data_(0), cols_(0)
62 {}
63
66 Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
67 {}
68
73 void setSize(size_type rows, size_type cols) {
74 data_.resize(rows,cols);
75 cols_ = cols;
76 }
77
80 {
81 return data_.begin();
82 }
83
86 {
87 return data_.end();
88 }
89
93 {
94 return data_.beforeEnd();
95 }
96
100 {
101 return data_.beforeBegin();
102 }
103
106 {
107 return data_.begin();
108 }
109
112 {
113 return data_.end();
114 }
115
119 {
120 return data_.beforeEnd();
121 }
122
126 {
127 return data_.beforeBegin();
128 }
129
132 {
133 data_ = t;
134 return *this;
135 }
136
139#ifdef DUNE_ISTL_WITH_CHECKING
140 if (row<0)
141 DUNE_THROW(ISTLError, "Can't access negative rows!");
142 if (row>=N())
143 DUNE_THROW(ISTLError, "Row index out of range!");
144#endif
145 return data_[row];
146 }
147
149 const row_type& operator[](size_type row) const {
150#ifdef DUNE_ISTL_WITH_CHECKING
151 if (row<0)
152 DUNE_THROW(ISTLError, "Can't access negative rows!");
153 if (row>=N())
154 DUNE_THROW(ISTLError, "Row index out of range!");
155#endif
156 return data_[row];
157 }
158
160 size_type N() const {
161 return data_.N();
162 }
163
165 size_type M() const {
166 return cols_;
167 }
168
171#ifdef DUNE_ISTL_WITH_CHECKING
172 if (M()==0)
173 DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
174#endif
175 size_type dim = 0;
176 for (size_type i=0; i<data_.N(); i++)
177 dim += data_[i][0].rowdim();
178
179 return dim;
180 }
181
184#ifdef DUNE_ISTL_WITH_CHECKING
185 if (N()==0)
186 DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
187#endif
188 size_type dim = 0;
189 for (size_type i=0; i<data_[0].size(); i++)
190 dim += data_[0][i].coldim();
191
192 return dim;
193 }
194
197#ifdef DUNE_ISTL_WITH_CHECKING
198 if (r<0 || r>=N())
199 DUNE_THROW(ISTLError, "Rowdim for nonexisting row " << r << " requested!");
200 if (M()==0)
201 DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
202#endif
203 return data_[r][0].rowdim();
204 }
205
208#ifdef DUNE_ISTL_WITH_CHECKING
209 if (c<0 || c>=M())
210 DUNE_THROW(ISTLError, "Coldim for nonexisting column " << c << " requested!");
211 if (N()==0)
212 DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
213#endif
214 return data_[0][c].coldim();
215 }
216
219 data_ *= scalar;
220 return (*this);
221 }
222
225 data_ /= scalar;
226 return (*this);
227 }
228
235#ifdef DUNE_ISTL_WITH_CHECKING
236 if(N()!=b.N() || M() != b.M())
237 DUNE_THROW(RangeError, "Matrix sizes do not match!");
238#endif
239 data_ += b.data_;
240 return (*this);
241 }
242
249#ifdef DUNE_ISTL_WITH_CHECKING
250 if(N()!=b.N() || M() != b.M())
251 DUNE_THROW(RangeError, "Matrix sizes do not match!");
252#endif
253 data_ -= b.data_;
254 return (*this);
255 }
256
259 Matrix out(N(), M());
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];
263
264 return out;
265 }
266
268 friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
269 Matrix<T> out(m1.N(), m2.M());
270 out = 0;
271
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];
276 }
277
278 return out;
279 }
280
282 template <class X, class Y>
283 friend Y operator*(const Matrix<T>& m, const X& vec) {
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!");
287#endif
288 Y out(m.N());
289 out = 0;
290
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];
294 }
295
296 return out;
297 }
298
300 template <class X, class Y>
301 void mv(const X& x, Y& y) const
302 {
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!");
306#endif
307
308 for (size_type i=0; i<data_.N(); i++) {
309 y[i]=0;
310 for (size_type j=0; j<cols_; j++)
311 (*this)[i][j].umv(x[j], y[i]);
312
313 }
314
315 }
316
318 template<class X, class Y>
319 void mtv (const X& x, Y& y) const
320 {
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");
324#endif
325
326 for(size_type i=0; i<y.N(); ++i)
327 y[i]=0;
328 umtv(x,y);
329 }
330
332 template <class X, class Y>
333 void umv(const X& x, Y& y) const
334 {
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!");
338#endif
339
340 for (size_type i=0; i<data_.N(); i++) {
341
342 for (size_type j=0; j<cols_; j++)
343 (*this)[i][j].umv(x[j], y[i]);
344
345 }
346
347 }
348
350 template<class X, class Y>
351 void mmv (const X& x, Y& y) const
352 {
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!");
356#endif
357 ConstRowIterator endi=end();
358 for (ConstRowIterator i=begin(); i!=endi; ++i)
359 {
360 ConstColIterator endj = (*i).end();
361 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
362 (*j).mmv(x[j.index()],y[i.index()]);
363 }
364 }
365
367 template <class X, class Y>
368 void usmv(const field_type& alpha, const X& x, Y& y) const
369 {
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!");
373#endif
374
375 for (size_type i=0; i<data_.N(); i++) {
376
377 for (size_type j=0; j<cols_; j++)
378 (*this)[i][j].usmv(alpha, x[j], y[i]);
379
380 }
381
382 }
383
385 template<class X, class Y>
386 void umtv (const X& x, Y& y) const
387 {
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!");
391#endif
392 ConstRowIterator endi=end();
393 for (ConstRowIterator i=begin(); i!=endi; ++i)
394 {
395 ConstColIterator endj = (*i).end();
396 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
397 (*j).umtv(x[i.index()],y[j.index()]);
398 }
399 }
400
402 template<class X, class Y>
403 void mmtv (const X& x, Y& y) const
404 {
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!");
408#endif
409 ConstRowIterator endi=end();
410 for (ConstRowIterator i=begin(); i!=endi; ++i)
411 {
412 ConstColIterator endj = (*i).end();
413 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
414 (*j).mmtv(x[i.index()],y[j.index()]);
415 }
416 }
417
419 template<class X, class Y>
420 void usmtv (const field_type& alpha, const X& x, Y& y) const
421 {
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!");
425#endif
426 ConstRowIterator endi=end();
427 for (ConstRowIterator i=begin(); i!=endi; ++i)
428 {
429 ConstColIterator endj = (*i).end();
430 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
431 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
432 }
433 }
434
436 template<class X, class Y>
437 void umhv (const X& x, Y& y) const
438 {
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!");
442#endif
443 ConstRowIterator endi=end();
444 for (ConstRowIterator i=begin(); i!=endi; ++i)
445 {
446 ConstColIterator endj = (*i).end();
447 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
448 (*j).umhv(x[i.index()],y[j.index()]);
449 }
450 }
451
453 template<class X, class Y>
454 void mmhv (const X& x, Y& y) const
455 {
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!");
459#endif
460 ConstRowIterator endi=end();
461 for (ConstRowIterator i=begin(); i!=endi; ++i)
462 {
463 ConstColIterator endj = (*i).end();
464 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
465 (*j).mmhv(x[i.index()],y[j.index()]);
466 }
467 }
468
470 template<class X, class Y>
471 void usmhv (const field_type& alpha, const X& x, Y& y) const
472 {
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!");
476#endif
477 ConstRowIterator endi=end();
478 for (ConstRowIterator i=begin(); i!=endi; ++i)
479 {
480 ConstColIterator endj = (*i).end();
481 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
482 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
483 }
484 }
485
486 //===== norms
487
489 typename FieldTraits<field_type>::real_type frobenius_norm () const
490 {
491 return std::sqrt(frobenius_norm2());
492 }
493
495 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
496 {
497 double sum=0;
498 for (size_type i=0; i<N(); ++i)
499 for (size_type j=0; j<M(); ++j)
500 sum += data_[i][j].frobenius_norm2();
501 return sum;
502 }
503
505 typename FieldTraits<field_type>::real_type infinity_norm () const
506 {
507 double max=0;
508 for (size_type i=0; i<N(); ++i) {
509 double sum=0;
510 for (size_type j=0; j<M(); j++)
511 sum += data_[i][j].infinity_norm();
512 max = std::max(max,sum);
513 }
514 return max;
515 }
516
518 typename FieldTraits<field_type>::real_type infinity_norm_real () const
519 {
520 double max=0;
521 for (size_type i=0; i<N(); ++i) {
522 double sum=0;
523 for (size_type j=0; j<M(); j++)
524 sum += data_[i][j].infinity_norm_real();
525 max = std::max(max,sum);
526 }
527 return max;
528 }
529
530 //===== query
531
533 bool exists (size_type i, size_type j) const
534 {
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");
538#else
540#endif
541 return true;
542 }
543
544 protected:
545
552
559 };
560
562} // end namespace Dune
563
564#endif
Definition: bvector.hh:584
derive error class from the base class in common
Definition: istlexception.hh:16
A generic dynamic dense matrix.
Definition: matrix.hh:25
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:505
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
@ blocklevel
The number of nesting levels the matrix contains.
Definition: matrix.hh:57
A allocator_type
Export the allocator.
Definition: matrix.hh:35
VariableBlockVector< T, A >::window_type row_type
The type implementing a matrix row.
Definition: matrix.hh:38
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: matrix.hh:471
void usmv(const field_type &alpha, const X &x, Y &y) const
Definition: matrix.hh:368
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:41
Matrix transpose() const
Return the transpose of the matrix.
Definition: matrix.hh:258
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
void mv(const X &x, Y &y) const
y = A x
Definition: matrix.hh:301
void setSize(size_type rows, size_type cols)
Change the matrix size.
Definition: matrix.hh:73
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:29
RowIterator beforeBegin()
Definition: matrix.hh:99
RowIterator beforeEnd()
Definition: matrix.hh:92
Matrix()
Create empty matrix.
Definition: matrix.hh:61
size_type coldim(size_type c) const
The number of scalar columns.
Definition: matrix.hh:207
Matrix & operator-=(const Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:248
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:495
VariableBlockVector< T, A >::Iterator RowIterator
Iterator over the matrix rows.
Definition: matrix.hh:44
row_type::iterator ColIterator
Iterator for the entries of each row.
Definition: matrix.hh:47
ConstRowIterator beforeEnd() const
Definition: matrix.hh:118
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:131
size_type rowdim() const
The number of scalar rows.
Definition: matrix.hh:170
size_type coldim() const
The number of scalar columns.
Definition: matrix.hh:183
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:85
ConstRowIterator end() const
Get const iterator to one beyond last row.
Definition: matrix.hh:111
friend Y operator*(const Matrix< T > &m, const X &vec)
Generic matrix-vector multiplication.
Definition: matrix.hh:283
Matrix< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:218
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:351
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:234
ConstRowIterator begin() const
Get const iterator to first row.
Definition: matrix.hh:105
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: matrix.hh:454
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
VariableBlockVector< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:50
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
VariableBlockVector< T, A > data_
Abuse VariableBlockVector as an engine for a 2d array ISTL-style.
Definition: matrix.hh:551
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:518
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
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
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: matrix.hh:533
Matrix< T > & operator/=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:224
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:268
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:386
const row_type & operator[](size_type row) const
The const index operator.
Definition: matrix.hh:149
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:403
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:489
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: matrix.hh:420
void umhv(const X &x, Y &y) const
y += A^H x
Definition: matrix.hh:437
size_type N() const
Return the number of rows.
Definition: matrix.hh:160
ConstRowIterator beforeBegin() const
Definition: matrix.hh:125
Matrix(size_type rows, size_type cols)
Create uninitialized matrix of size rows x cols.
Definition: matrix.hh:66
Default exception class for range errors.
Definition: exceptions.hh:279
ConstIterator class for sequential access.
Definition: vbvector.hh:647
Iterator class for sequential access.
Definition: vbvector.hh:525
A Vector of blocks with different blocksizes.
Definition: vbvector.hh:40
RealIterator< B > iterator
iterator type for sequential access
Definition: basearray.hh:156
RealIterator< const B > const_iterator
iterator class for sequential access
Definition: basearray.hh:195
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Dune namespace.
Definition: alignment.hh:10
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)