dune-istl  2.2.1
diagonalmatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 #ifndef DUNE_DIAGONAL_MATRIX_HH
4 #define DUNE_DIAGONAL_MATRIX_HH
5 
10 #include<cmath>
11 #include<cstddef>
12 #include<complex>
13 #include<iostream>
14 #include<memory>
15 #include <dune/common/exceptions.hh>
16 #include <dune/common/fmatrix.hh>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/typetraits.hh>
19 #include <dune/common/genericiterator.hh>
20 
21 
22 
23 namespace Dune {
24 
25 template< class K, int n > class DiagonalRowVectorConst;
26 template< class K, int n > class DiagonalRowVector;
27 template< class DiagonalMatrixType > class DiagonalMatrixWrapper;
28 template< class C, class T, class R> class ContainerWrapperIterator;
29 
44 template<class K, int n>
46 {
48 
49  public:
50  //===== type definitions and constants
51 
53  typedef K field_type;
54 
56  typedef K block_type;
57 
59  typedef std::size_t size_type;
60 
62  enum {
65  };
66 
72 
74  enum {
76  rows = n,
78  cols = n
79  };
80 
81 
82 
83  //===== constructors
84 
87 
89  DiagonalMatrix (const K& k)
90  : diag_(k)
91  {}
92 
94  DiagonalMatrix (const FieldVector<K,n>& diag)
95  : diag_(diag)
96  {}
97 
98 
101  {
102  diag_ = k;
103  return *this;
104  }
105 
107  bool identical(const DiagonalMatrix<K,n>& other) const
108  {
109  return (this==&other);
110  }
111 
112  //===== iterator interface to rows of the matrix
120  typedef typename row_type::Iterator ColIterator;
121 
124  {
125  return Iterator(WrapperType(this),0);
126  }
127 
130  {
131  return Iterator(WrapperType(this),n);
132  }
133 
137  {
138  return Iterator(WrapperType(this),n-1);
139  }
140 
144  {
145  return Iterator(WrapperType(this),-1);
146  }
147 
148 
157 
160  {
161  return ConstIterator(WrapperType(this),0);
162  }
163 
166  {
167  return ConstIterator(WrapperType(this),n);
168  }
169 
173  {
174  return ConstIterator(WrapperType(this),n-1);
175  }
176 
180  {
181  return ConstIterator(WrapperType(this),-1);
182  }
183 
184 
185 
186  //===== vector space arithmetic
187 
190  {
191  diag_ += y.diag_;
192  return *this;
193  }
194 
197  {
198  diag_ -= y.diag_;
199  return *this;
200  }
201 
204  {
205  diag_ += k;
206  return *this;
207  }
208 
211  {
212  diag_ -= k;
213  return *this;
214  }
215 
218  {
219  diag_ *= k;
220  return *this;
221  }
222 
225  {
226  diag_ /= k;
227  return *this;
228  }
229 
230  //===== comparison ops
231 
233  bool operator==(const DiagonalMatrix& other) const
234  {
235  return diag_==other.diagonal();
236  }
237 
239  bool operator!=(const DiagonalMatrix& other) const
240  {
241  return diag_!=other.diagonal();
242  }
243 
244 
245  //===== linear maps
246 
248  template<class X, class Y>
249  void mv (const X& x, Y& y) const
250  {
251 #ifdef DUNE_FMatrix_WITH_CHECKING
252  if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
253  if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
254 #endif
255  for (size_type i=0; i<n; ++i)
256  y[i] = diag_[i] * x[i];
257  }
258 
260  template<class X, class Y>
261  void mtv (const X& x, Y& y) const
262  {
263  mv(x, y);
264  }
265 
267  template<class X, class Y>
268  void umv (const X& x, Y& y) const
269  {
270 #ifdef DUNE_FMatrix_WITH_CHECKING
271  if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
272  if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
273 #endif
274  for (size_type i=0; i<n; ++i)
275  y[i] += diag_[i] * x[i];
276  }
277 
279  template<class X, class Y>
280  void umtv (const X& x, Y& y) const
281  {
282 #ifdef DUNE_FMatrix_WITH_CHECKING
283  if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
284  if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
285 #endif
286  for (size_type i=0; i<n; ++i)
287  y[i] += diag_[i] * x[i];
288  }
289 
291  template<class X, class Y>
292  void umhv (const X& x, Y& y) const
293  {
294 #ifdef DUNE_FMatrix_WITH_CHECKING
295  if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
296  if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
297 #endif
298  for (size_type i=0; i<n; i++)
299  y[i] += conjugateComplex(diag_[i])*x[i];
300  }
301 
303  template<class X, class Y>
304  void mmv (const X& x, Y& y) const
305  {
306 #ifdef DUNE_FMatrix_WITH_CHECKING
307  if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
308  if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
309 #endif
310  for (size_type i=0; i<n; ++i)
311  y[i] -= diag_[i] * x[i];
312  }
313 
315  template<class X, class Y>
316  void mmtv (const X& x, Y& y) const
317  {
318 #ifdef DUNE_FMatrix_WITH_CHECKING
319  if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
320  if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
321 #endif
322  for (size_type i=0; i<n; ++i)
323  y[i] -= diag_[i] * x[i];
324  }
325 
327  template<class X, class Y>
328  void mmhv (const X& x, Y& y) const
329  {
330 #ifdef DUNE_FMatrix_WITH_CHECKING
331  if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
332  if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
333 #endif
334  for (size_type i=0; i<n; i++)
335  y[i] -= conjugateComplex(diag_[i])*x[i];
336  }
337 
339  template<class X, class Y>
340  void usmv (const K& alpha, const X& x, Y& y) const
341  {
342 #ifdef DUNE_FMatrix_WITH_CHECKING
343  if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
344  if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
345 #endif
346  for (size_type i=0; i<n; i++)
347  y[i] += alpha * diag_[i] * x[i];
348  }
349 
351  template<class X, class Y>
352  void usmtv (const K& alpha, const X& x, Y& y) const
353  {
354 #ifdef DUNE_FMatrix_WITH_CHECKING
355  if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
356  if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
357 #endif
358  for (size_type i=0; i<n; i++)
359  y[i] += alpha * diag_[i] * x[i];
360  }
361 
363  template<class X, class Y>
364  void usmhv (const K& alpha, const X& x, Y& y) const
365  {
366 #ifdef DUNE_FMatrix_WITH_CHECKING
367  if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
368  if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
369 #endif
370  for (size_type i=0; i<n; i++)
371  y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
372  }
373 
374  //===== norms
375 
377  double frobenius_norm () const
378  {
379  return diag_.two_norm();
380  }
381 
383  double frobenius_norm2 () const
384  {
385  return diag_.two_norm2();
386  }
387 
389  double infinity_norm () const
390  {
391  return diag_.infinity_norm();
392  }
393 
395  double infinity_norm_real () const
396  {
397  return diag_.infinity_norm_real();
398  }
399 
400 
401 
402  //===== solve
403 
405  template<class V>
406  void solve (V& x, const V& b) const
407  {
408  for (int i=0; i<n; i++)
409  x[i] = b[i]/diag_[i];
410  }
411 
413  void invert()
414  {
415  for (int i=0; i<n; i++)
416  diag_[i] = 1/diag_[i];
417  }
418 
420  K determinant () const
421  {
422  K det = diag_[0];
423  for (int i=1; i<n; i++)
424  det *= diag_[i];
425  return det;
426  }
427 
428 
429 
430  //===== sizes
431 
433  size_type N () const
434  {
435  return n;
436  }
437 
439  size_type M () const
440  {
441  return n;
442  }
443 
444 
445 
446  //===== query
447 
449  bool exists (size_type i, size_type j) const
450  {
451 #ifdef DUNE_FMatrix_WITH_CHECKING
452  if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range");
453  if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range");
454 #endif
455  return i==j;
456  }
457 
458 
459 
461  friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
462  {
463  for (size_type i=0; i<n; i++) {
464  for (size_type j=0; j<n; j++)
465  s << ((i==j) ? a.diag_[i] : 0) << " ";
466  s << std::endl;
467  }
468  return s;
469  }
470 
473  {
474  return reference(const_cast<K*>(&diag_[i]), i);
475  }
476 
479  {
480  return const_reference(const_cast<K*>(&diag_[i]), i);
481  }
482 
484  const K& diagonal(size_type i) const
485  {
486  return diag_[i];
487  }
488 
491  {
492  return diag_[i];
493  }
494 
496  const FieldVector<K,n>& diagonal() const
497  {
498  return diag_;
499  }
500 
502  FieldVector<K,n>& diagonal()
503  {
504  return diag_;
505  }
506 
507  private:
508 
509  // the data, a FieldVector storing the diagonal
510  FieldVector<K,n> diag_;
511 };
512 
513 #ifndef DOXYGEN // hide specialization
514 
516 template< class K >
517 class DiagonalMatrix<K, 1> : public FieldMatrix<K, 1, 1>
518 {
519  typedef FieldMatrix<K,1,1> Base;
520  public:
522  typedef typename Base::size_type size_type;
523 
525  enum {
528  blocklevel = 1
529  };
530 
531  typedef typename Base::row_type row_type;
532 
533  typedef typename Base::row_reference row_reference;
534  typedef typename Base::const_row_reference const_row_reference;
535 
537  enum {
540  rows = 1,
543  cols = 1
544  };
545 
546 
549  {}
550 
552  DiagonalMatrix(const K& scalar)
553  {
554  (*this)[0][0] = scalar;
555  }
556 
558  template<typename T>
559  DiagonalMatrix(const T& t)
560  {
561  DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
562  }
563 
565  const K& diagonal(size_type) const
566  {
567  return (*this)[0][0];
568  }
569 
571  K& diagonal(size_type)
572  {
573  return (*this)[0][0];
574  }
575 
577  const FieldVector<K,1>& diagonal() const
578  {
579  return (*this)[0];
580  }
581 
583  FieldVector<K,1>& diagonal()
584  {
585  return (*this)[0];
586  }
587 
588 };
589 #endif
590 
591 
592 template<class DiagonalMatrixType>
593 class DiagonalMatrixWrapper
594 {
595  typedef typename DiagonalMatrixType::reference reference;
596  typedef typename DiagonalMatrixType::const_reference const_reference;
597  typedef typename DiagonalMatrixType::field_type K;
598  typedef DiagonalRowVector<K, DiagonalMatrixType::rows> row_type;
599  typedef std::size_t size_type;
600  typedef DiagonalMatrixWrapper< DiagonalMatrixType> MyType;
601 
602  friend class ContainerWrapperIterator<const MyType, reference, reference>;
603  friend class ContainerWrapperIterator<const MyType, const_reference, const_reference>;
604 
605  public:
606 
608  mat_(0)
609  {}
610 
611  DiagonalMatrixWrapper(const DiagonalMatrixType* mat) :
612  mat_(const_cast<DiagonalMatrixType*>(mat))
613  {}
614 
615  size_type realIndex(int i) const
616  {
617  return i;
618  }
619 
620  row_type* pointer(int i) const
621  {
622  row_ = row_type(&(mat_->diagonal(i)), i);
623  return &row_;
624  }
625 
626  bool identical(const DiagonalMatrixWrapper& other) const
627  {
628  return mat_==other.mat_;
629  }
630 
631  private:
632 
633  mutable DiagonalMatrixType* mat_;
634  mutable row_type row_;
635 };
636 
640 template< class K, int n >
641 class DiagonalRowVectorConst
642 {
643  template<class DiagonalMatrixType>
644  friend class DiagonalMatrixWrapper;
645  friend class ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&>;
646 
647 public:
648  // remember size of vector
649  enum { dimension = n };
650 
651  // standard constructor and everything is sufficient ...
652 
653  //===== type definitions and constants
654 
656  typedef K field_type;
657 
659  typedef K block_type;
660 
662  typedef std::size_t size_type;
663 
665  enum {
668  };
669 
671  enum {
673  size = n
674  };
675 
678  p_(0),
679  row_(0)
680  {}
681 
683  explicit DiagonalRowVectorConst (K* p, int col) :
684  p_(p),
685  row_(col)
686  {}
687 
688  //===== access to components
689 
691  const K& operator[] (size_type i) const
692  {
693 #ifdef DUNE_FMatrix_WITH_CHECKING
694  if (i!=row_)
695  DUNE_THROW(FMatrixError,"index is not contained in pattern");
696 #endif
697  return *p_;
698  }
699 
700  // check if row is identical to other row (not only identical values)
701  // since this is a proxy class we need to check equality of the stored pointer
702  bool identical(const DiagonalRowVectorConst<K,n>& other) const
703  {
704  return ((p_ == other.p_) and (row_ == other.row_));
705  }
706 
711 
714  {
715  return ConstIterator(*this,0);
716  }
717 
720  {
721  return ConstIterator(*this,1);
722  }
723 
727  {
728  return ConstIterator(*this,0);
729  }
730 
734  {
735  return ConstIterator(*this,-1);
736  }
737 
739  bool operator== (const DiagonalRowVectorConst& y) const
740  {
741  return ((p_==y.p_) and (row_==y.row_));
742  }
743 
744  //===== sizes
745 
747  size_type N () const
748  {
749  return n;
750  }
751 
753  size_type dim () const
754  {
755  return n;
756  }
757 
760  {
761  return row_;
762  }
763 
765  const K& diagonal() const
766  {
767  return *p_;
768  }
769 
770 protected:
771 
772  size_type realIndex(int i) const
773  {
774  return rowIndex();
775  }
776 
777  K* pointer(size_type i) const
778  {
779  return const_cast<K*>(p_);
780  }
781 
783  {
784  return this;
785  }
786 
787  // the data, very simply a pointer to the diagonal value and the row number
788  K* p_;
790 };
791 
792 template< class K, int n >
793 class DiagonalRowVector : public DiagonalRowVectorConst<K,n>
794 {
795  template<class DiagonalMatrixType>
796  friend class DiagonalMatrixWrapper;
797  friend class ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&>;
798 
799 public:
800  // standard constructor and everything is sufficient ...
801 
802  //===== type definitions and constants
803 
805  typedef K field_type;
806 
808  typedef K block_type;
809 
811  typedef std::size_t size_type;
812 
815  {}
816 
818  explicit DiagonalRowVector (K* p, int col) : DiagonalRowVectorConst<K,n>(p, col)
819  {}
820 
821  //===== assignment from scalar
824  {
825  *p_ = k;
826  return *this;
827  }
828 
829  //===== access to components
830 
833  {
834 #ifdef DUNE_FMatrix_WITH_CHECKING
835  if (i!=row_)
836  DUNE_THROW(FMatrixError,"index is contained in pattern");
837 #endif
838  return *p_;
839  }
840 
845 
848  {
849  return Iterator(*this, 0);
850  }
851 
854  {
855  return Iterator(*this, 1);
856  }
857 
861  {
862  return Iterator(*this, 0);
863  }
864 
868  {
869  return Iterator(*this, -1);
870  }
871 
876 
888 
889 protected:
890 
892  {
893  return this;
894  }
895 
896 private:
897 
900 };
901 
902 
903 // implement type traits
904 template<class K, int n>
905 struct const_reference< DiagonalRowVector<K,n> >
906 {
908 };
909 
910 template<class K, int n>
911 struct const_reference< DiagonalRowVectorConst<K,n> >
912 {
914 };
915 
916 template<class K, int n>
917 struct mutable_reference< DiagonalRowVector<K,n> >
918 {
920 };
921 
922 template<class K, int n>
923 struct mutable_reference< DiagonalRowVectorConst<K,n> >
924 {
926 };
927 
928 
929 
952 template<class CW, class T, class R>
953 class ContainerWrapperIterator : public BidirectionalIteratorFacade<ContainerWrapperIterator<CW,T,R>,T, R, int>
954 {
955  typedef typename remove_const<CW>::type NonConstCW;
956 
957  friend class ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type>;
958  friend class ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type>;
959 
960  typedef ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type> MyType;
961  typedef ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type> MyConstType;
962 
963  public:
964 
965  // Constructors needed by the facade iterators.
967  containerWrapper_(),
968  position_(0)
969  {}
970 
971  ContainerWrapperIterator(CW containerWrapper, int position) :
972  containerWrapper_(containerWrapper),
973  position_(position)
974  {}
975 
976  template<class OtherContainerWrapperIteratorType>
977  ContainerWrapperIterator(OtherContainerWrapperIteratorType& other):
978  containerWrapper_(other.containerWrapper_),
979  position_(other.position_)
980  {}
981 
983  containerWrapper_(other.containerWrapper_),
984  position_(other.position_)
985  {}
986 
988  containerWrapper_(other.containerWrapper_),
989  position_(other.position_)
990  {}
991 
992  template<class OtherContainerWrapperIteratorType>
993  ContainerWrapperIterator& operator=(OtherContainerWrapperIteratorType& other)
994  {
995  containerWrapper_ = other.containerWrapper_;
996  position_ = other.position_;
997  }
998 
999  // This operator is needed since we can not get the address of the
1000  // temporary object returned by dereference
1001  T* operator->() const
1002  {
1003  return containerWrapper_.pointer(position_);
1004  }
1005 
1006  // Methods needed by the forward iterator
1007  bool equals(const MyType& other) const
1008  {
1009  return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1010  }
1011 
1012  bool equals(const MyConstType& other) const
1013  {
1014  return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1015  }
1016 
1017  R dereference() const
1018  {
1019  return *containerWrapper_.pointer(position_);
1020  }
1021 
1022  void increment()
1023  {
1024  ++position_;
1025  }
1026 
1027  // Additional function needed by BidirectionalIterator
1028  void decrement()
1029  {
1030  --position_;
1031  }
1032 
1033  // Additional function needed by RandomAccessIterator
1034  R elementAt(int i) const
1035  {
1036  return *containerWrapper_.pointer(position_+i);
1037  }
1038 
1039  void advance(int n)
1040  {
1041  position_=position_+n;
1042  }
1043 
1044  template<class OtherContainerWrapperIteratorType>
1045  std::ptrdiff_t distanceTo(OtherContainerWrapperIteratorType& other) const
1046  {
1047  assert(containerWrapper_.identical(other));
1048  return other.position_ - position_;
1049  }
1050 
1051  std::ptrdiff_t index() const
1052  {
1053  return containerWrapper_.realIndex(position_);
1054  }
1055 
1056  private:
1057  NonConstCW containerWrapper_;
1058  size_t position_;
1059 };
1060 
1061 
1062 
1063 template<class M, class K, int n>
1064 void istl_assign_to_fmatrix(DenseMatrix<M>& fm, const DiagonalMatrix<K,n>& s)
1065 {
1066  fm = K();
1067  for(int i=0; i<n; ++i)
1068  fm[i][i] = s.diagonal()[i];
1069 }
1070 /* @} */
1071 } // end namespace
1072 #endif