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>
14 #include <dune/common/ftraits.hh>
15 
16 namespace 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 
47  typedef typename row_type::iterator ColIterator;
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 
170  size_type rowdim() const {
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 
183  size_type coldim() const {
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 
218  Matrix<T>& operator*=(const field_type& scalar) {
219  data_ *= scalar;
220  return (*this);
221  }
222 
224  Matrix<T>& operator/=(const field_type& scalar) {
225  data_ /= scalar;
226  return (*this);
227  }
228 
234  Matrix& operator+= (const Matrix& b) {
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 
248  Matrix& operator-= (const Matrix& b) {
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 
258  Matrix transpose() const {
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
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
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: matrix.hh:489
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
Matrix & operator=(const field_type &t)
Assignment from scalar.
Definition: matrix.hh:131
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
Matrix< T > & operator/=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:224
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
row_type & operator[](size_type row)
The index operator.
Definition: matrix.hh:138
const row_type & operator[](size_type row) const
The const index operator.
Definition: matrix.hh:149
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 Matrix &b)
Subtract the entries of another matrix from this one.
Definition: matrix.hh:248
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
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: matrix.hh:505
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
void mmv(const X &x, Y &y) const
y -= A x
Definition: matrix.hh:351
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
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:518
VariableBlockVector< T, A > data_
Abuse VariableBlockVector as an engine for a 2d array ISTL-style.
Definition: matrix.hh:551
friend Matrix< T > operator*(const Matrix< T > &m1, const Matrix< T > &m2)
Generic matrix multiplication.
Definition: matrix.hh:268
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
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
void umtv(const X &x, Y &y) const
y += A^T x
Definition: matrix.hh:386
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: matrix.hh:495
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: matrix.hh:403
Matrix & operator+=(const Matrix &b)
Add the entries of another matrix to this one.
Definition: matrix.hh:234
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< T > & operator*=(const field_type &scalar)
Multiplication with a scalar.
Definition: matrix.hh:218
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.80.0 (Apr 18, 22:30, 2024)