dune-istl  2.2.0
matrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=8 sw=4 sts=4:
3 #ifndef DUNE_MATRIX_HH
4 #define DUNE_MATRIX_HH
5 
10 #include <vector>
11 #include <memory>
12 
13 #include <dune/istl/vbvector.hh>
14 
15 namespace Dune {
16 
22  template<class T, class A=std::allocator<T> >
23 class Matrix
24 {
25 public:
26 
28  typedef typename T::field_type field_type;
29 
31  typedef T block_type;
32 
34  typedef A allocator_type;
35 
38 
40  typedef typename A::size_type size_type;
41 
44 
46  typedef typename row_type::iterator ColIterator;
47 
50 
52  typedef typename row_type::const_iterator ConstColIterator;
53 
54  enum {
56  blocklevel = T::blocklevel+1
57  };
58 
60  Matrix() : data_(0), cols_(0)
61  {}
62 
65  Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols)
66  {}
67 
72  void setSize(size_type rows, size_type cols) {
73  data_.resize(rows,cols);
74  cols_ = cols;
75  }
76 
79  {
80  return data_.begin();
81  }
82 
85  {
86  return data_.end();
87  }
88 
92  {
93  return data_.beforeEnd();
94  }
95 
99  {
100  return data_.beforeBegin();
101  }
102 
105  {
106  return data_.begin();
107  }
108 
111  {
112  return data_.end();
113  }
114 
118  {
119  return data_.beforeEnd();
120  }
121 
125  {
126  return data_.beforeBegin();
127  }
128 
131  {
132  data_ = t;
133  return *this;
134  }
135 
138 #ifdef DUNE_ISTL_WITH_CHECKING
139  if (row<0)
140  DUNE_THROW(ISTLError, "Can't access negative rows!");
141  if (row>=N())
142  DUNE_THROW(ISTLError, "Row index out of range!");
143 #endif
144  return data_[row];
145  }
146 
149 #ifdef DUNE_ISTL_WITH_CHECKING
150  if (row<0)
151  DUNE_THROW(ISTLError, "Can't access negative rows!");
152  if (row>=N())
153  DUNE_THROW(ISTLError, "Row index out of range!");
154 #endif
155  return data_[row];
156  }
157 
159  size_type N() const {
160  return data_.N();
161  }
162 
164  size_type M() const {
165  return cols_;
166  }
167 
169  size_type rowdim() const {
170 #ifdef DUNE_ISTL_WITH_CHECKING
171  if (M()==0)
172  DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
173 #endif
174  size_type dim = 0;
175  for (size_type i=0; i<data_.N(); i++)
176  dim += data_[i][0].rowdim();
177 
178  return dim;
179  }
180 
182  size_type coldim() const {
183 #ifdef DUNE_ISTL_WITH_CHECKING
184  if (N()==0)
185  DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
186 #endif
187  size_type dim = 0;
188  for (size_type i=0; i<data_[0].size(); i++)
189  dim += data_[0][i].coldim();
190 
191  return dim;
192  }
193 
196 #ifdef DUNE_ISTL_WITH_CHECKING
197  if (r<0 || r>=N())
198  DUNE_THROW(ISTLError, "Rowdim for nonexisting row " << r << " requested!");
199  if (M()==0)
200  DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!");
201 #endif
202  return data_[r][0].rowdim();
203  }
204 
207 #ifdef DUNE_ISTL_WITH_CHECKING
208  if (c<0 || c>=M())
209  DUNE_THROW(ISTLError, "Coldim for nonexisting column " << c << " requested!");
210  if (N()==0)
211  DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!");
212 #endif
213  return data_[0][c].coldim();
214  }
215 
217  Matrix<T>& operator*=(const field_type& scalar) {
218  data_ *= scalar;
219  return (*this);
220  }
221 
223  Matrix<T>& operator/=(const field_type& scalar) {
224  data_ /= scalar;
225  return (*this);
226  }
227 
233  Matrix& operator+= (const Matrix& b) {
234 #ifdef DUNE_ISTL_WITH_CHECKING
235  if(N()!=b.N() || M() != b.M())
236  DUNE_THROW(RangeError, "Matrix sizes do not match!");
237 #endif
238  data_ += b.data_;
239  return (*this);
240  }
241 
247  Matrix& operator-= (const Matrix& b) {
248 #ifdef DUNE_ISTL_WITH_CHECKING
249  if(N()!=b.N() || M() != b.M())
250  DUNE_THROW(RangeError, "Matrix sizes do not match!");
251 #endif
252  data_ -= b.data_;
253  return (*this);
254  }
255 
257  Matrix transpose() const {
258  Matrix out(N(), M());
259  for (size_type i=0; i<M(); i++)
260  for (size_type j=0; j<N(); j++)
261  out[j][i] = (*this)[i][j];
262 
263  return out;
264  }
265 
269  template <class X, class Y>
270  Y transposedMult(const X& vec) DUNE_DEPRECATED;
271 
273  friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {
274  Matrix<T> out(m1.N(), m2.M());
275  out = 0;
276 
277  for (size_type i=0; i<out.N(); i++ ) {
278  for ( size_type j=0; j<out.M(); j++ )
279  for (size_type k=0; k<m1.M(); k++)
280  out[i][j] += m1[i][k]*m2[k][j];
281  }
282 
283  return out;
284  }
285 
287  template <class X, class Y>
288  friend Y operator*(const Matrix<T>& m, const X& vec) {
289 #ifdef DUNE_ISTL_WITH_CHECKING
290  if (m.M()!=vec.size())
291  DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix columns!");
292 #endif
293  Y out(m.N());
294  out = 0;
295 
296  for (size_type i=0; i<out.size(); i++ ) {
297  for ( size_type j=0; j<vec.size(); j++ )
298  out[i] += m[i][j]*vec[j];
299  }
300 
301  return out;
302  }
303 
305  template <class X, class Y>
306  void mv(const X& x, Y& y) const
307  {
308 #ifdef DUNE_ISTL_WITH_CHECKING
309  if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
310  if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
311 #endif
312 
313  for (size_type i=0; i<data_.N(); i++) {
314  y[i]=0;
315  for (size_type j=0; j<cols_; j++)
316  (*this)[i][j].umv(x[j], y[i]);
317 
318  }
319 
320  }
321 
323  template<class X, class Y>
324  void mtv (const X& x, Y& y) const
325  {
326 #ifdef DUNE_ISTL_WITH_CHECKING
327  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
328  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
329 #endif
330 
331  for(size_type i=0; i<y.N(); ++i)
332  y[i]=0;
333  umtv(x,y);
334  }
335 
337  template <class X, class Y>
338  void umv(const X& x, Y& y) const
339  {
340 #ifdef DUNE_ISTL_WITH_CHECKING
341  if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
342  if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
343 #endif
344 
345  for (size_type i=0; i<data_.N(); i++) {
346 
347  for (size_type j=0; j<cols_; j++)
348  (*this)[i][j].umv(x[j], y[i]);
349 
350  }
351 
352  }
353 
355  template<class X, class Y>
356  void mmv (const X& x, Y& y) const
357  {
358 #ifdef DUNE_ISTL_WITH_CHECKING
359  if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
360  if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
361 #endif
362  ConstRowIterator endi=end();
363  for (ConstRowIterator i=begin(); i!=endi; ++i)
364  {
365  ConstColIterator endj = (*i).end();
366  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
367  (*j).mmv(x[j.index()],y[i.index()]);
368  }
369  }
370 
372  template <class X, class Y>
373  void usmv(const field_type& alpha, const X& x, Y& y) const
374  {
375 #ifdef DUNE_ISTL_WITH_CHECKING
376  if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
377  if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
378 #endif
379 
380  for (size_type i=0; i<data_.N(); i++) {
381 
382  for (size_type j=0; j<cols_; j++)
383  (*this)[i][j].usmv(alpha, x[j], y[i]);
384 
385  }
386 
387  }
388 
390  template<class X, class Y>
391  void umtv (const X& x, Y& y) const
392  {
393 #ifdef DUNE_ISTL_WITH_CHECKING
394  if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
395  if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
396 #endif
397  ConstRowIterator endi=end();
398  for (ConstRowIterator i=begin(); i!=endi; ++i)
399  {
400  ConstColIterator endj = (*i).end();
401  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
402  (*j).umtv(x[i.index()],y[j.index()]);
403  }
404  }
405 
407  template<class X, class Y>
408  void mmtv (const X& x, Y& y) const
409  {
410 #ifdef DUNE_ISTL_WITH_CHECKING
411  if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
412  if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
413 #endif
414  ConstRowIterator endi=end();
415  for (ConstRowIterator i=begin(); i!=endi; ++i)
416  {
417  ConstColIterator endj = (*i).end();
418  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
419  (*j).mmtv(x[i.index()],y[j.index()]);
420  }
421  }
422 
424  template<class X, class Y>
425  void usmtv (const field_type& alpha, const X& x, Y& y) const
426  {
427 #ifdef DUNE_ISTL_WITH_CHECKING
428  if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
429  if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
430 #endif
431  ConstRowIterator endi=end();
432  for (ConstRowIterator i=begin(); i!=endi; ++i)
433  {
434  ConstColIterator endj = (*i).end();
435  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
436  (*j).usmtv(alpha,x[i.index()],y[j.index()]);
437  }
438  }
439 
441  template<class X, class Y>
442  void umhv (const X& x, Y& y) const
443  {
444 #ifdef DUNE_ISTL_WITH_CHECKING
445  if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
446  if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
447 #endif
448  ConstRowIterator endi=end();
449  for (ConstRowIterator i=begin(); i!=endi; ++i)
450  {
451  ConstColIterator endj = (*i).end();
452  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
453  (*j).umhv(x[i.index()],y[j.index()]);
454  }
455  }
456 
458  template<class X, class Y>
459  void mmhv (const X& x, Y& y) const
460  {
461 #ifdef DUNE_ISTL_WITH_CHECKING
462  if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
463  if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
464 #endif
465  ConstRowIterator endi=end();
466  for (ConstRowIterator i=begin(); i!=endi; ++i)
467  {
468  ConstColIterator endj = (*i).end();
469  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
470  (*j).mmhv(x[i.index()],y[j.index()]);
471  }
472  }
473 
475  template<class X, class Y>
476  void usmhv (const field_type& alpha, const X& x, Y& y) const
477  {
478 #ifdef DUNE_ISTL_WITH_CHECKING
479  if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
480  if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!");
481 #endif
482  ConstRowIterator endi=end();
483  for (ConstRowIterator i=begin(); i!=endi; ++i)
484  {
485  ConstColIterator endj = (*i).end();
486  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
487  (*j).usmhv(alpha,x[i.index()],y[j.index()]);
488  }
489  }
490 
491  //===== norms
492 
494  double frobenius_norm () const
495  {
496  return std::sqrt(frobenius_norm2());
497  }
498 
500  double frobenius_norm2 () const
501  {
502  double sum=0;
503  for (size_type i=0; i<N(); ++i)
504  for (size_type j=0; j<M(); ++j)
505  sum += data_[i][j].frobenius_norm2();
506  return sum;
507  }
508 
510  double infinity_norm () const
511  {
512  double max=0;
513  for (size_type i=0; i<N(); ++i) {
514  double sum=0;
515  for (size_type j=0; j<M(); j++)
516  sum += data_[i][j].infinity_norm();
517  max = std::max(max,sum);
518  }
519  return max;
520  }
521 
523  double infinity_norm_real () const
524  {
525  double max=0;
526  for (size_type i=0; i<N(); ++i) {
527  double sum=0;
528  for (size_type j=0; j<M(); j++)
529  sum += data_[i][j].infinity_norm_real();
530  max = std::max(max,sum);
531  }
532  return max;
533  }
534 
535  //===== query
536 
538  bool exists (size_type i, size_type j) const
539  {
540 #ifdef DUNE_ISTL_WITH_CHECKING
541  if (i<0 || i>=N()) DUNE_THROW(ISTLError,"row index out of range");
542  if (j<0 || i>=M()) DUNE_THROW(ISTLError,"column index out of range");
543 #endif
544  return true;
545  }
546 
547 protected:
548 
555 
562 };
563 
564 template<class T, class A>
565 template<class X, class Y>
566 inline Y Matrix<T,A>::transposedMult(const X& vec)
567 {
568 #ifdef DUNE_ISTL_WITH_CHECKING
569  if (N()!=vec.size())
570  DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix rows!");
571 #endif
572  Y out(M());
573  out = 0;
574 
575  for (size_type i=0; i<out.size(); i++ ) {
576  for ( size_type j=0; j<vec.size(); j++ )
577  out[i] += (*this)[j][i]*vec[j];
578  }
579 
580  return out;
581 }
582 
584 } // end namespace Dune
585 
586 #endif