3 #ifndef DUNE_DIAGONAL_MATRIX_HH
4 #define DUNE_DIAGONAL_MATRIX_HH
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>
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;
44 template<
class K,
int n>
109 return (
this==&other);
248 template<
class X,
class Y>
249 void mv (
const X& x, Y& y)
const
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");
256 y[i] = diag_[i] * x[i];
260 template<
class X,
class Y>
261 void mtv (
const X& x, Y& y)
const
267 template<
class X,
class Y>
268 void umv (
const X& x, Y& y)
const
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");
275 y[i] += diag_[i] * x[i];
279 template<
class X,
class Y>
280 void umtv (
const X& x, Y& y)
const
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");
287 y[i] += diag_[i] * x[i];
291 template<
class X,
class Y>
292 void umhv (
const X& x, Y& y)
const
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");
299 y[i] += conjugateComplex(diag_[i])*x[i];
303 template<
class X,
class Y>
304 void mmv (
const X& x, Y& y)
const
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");
311 y[i] -= diag_[i] * x[i];
315 template<
class X,
class Y>
316 void mmtv (
const X& x, Y& y)
const
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");
323 y[i] -= diag_[i] * x[i];
327 template<
class X,
class Y>
328 void mmhv (
const X& x, Y& y)
const
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");
335 y[i] -= conjugateComplex(diag_[i])*x[i];
339 template<
class X,
class Y>
340 void usmv (
const K& alpha,
const X& x, Y& y)
const
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");
347 y[i] += alpha * diag_[i] * x[i];
351 template<
class X,
class Y>
352 void usmtv (
const K& alpha,
const X& x, Y& y)
const
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");
359 y[i] += alpha * diag_[i] * x[i];
363 template<
class X,
class Y>
364 void usmhv (
const K& alpha,
const X& x, Y& y)
const
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");
371 y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
379 return diag_.two_norm();
385 return diag_.two_norm2();
391 return diag_.infinity_norm();
397 return diag_.infinity_norm_real();
408 for (
int i=0; i<n; i++)
409 x[i] = b[i]/diag_[i];
415 for (
int i=0; i<n; i++)
416 diag_[i] = 1/diag_[i];
423 for (
int i=1; i<n; i++)
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");
461 friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
465 s << ((i==j) ? a.diag_[i] : 0) <<
" ";
474 return reference(const_cast<K*>(&diag_[i]), i);
510 FieldVector<K,n> diag_;
513 #ifndef DOXYGEN // hide specialization
517 class DiagonalMatrix<K, 1> :
public FieldMatrix<K, 1, 1>
519 typedef FieldMatrix<K,1,1> Base;
522 typedef typename Base::size_type
size_type;
531 typedef typename Base::row_type
row_type;
533 typedef typename Base::row_reference row_reference;
534 typedef typename Base::const_row_reference const_row_reference;
554 (*this)[0][0] = scalar;
561 DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*
this, t);
567 return (*
this)[0][0];
573 return (*
this)[0][0];
577 const FieldVector<K,1>&
diagonal()
const
592 template<
class DiagonalMatrixType>
595 typedef typename DiagonalMatrixType::reference reference;
596 typedef typename DiagonalMatrixType::const_reference const_reference;
597 typedef typename DiagonalMatrixType::field_type K;
599 typedef std::size_t size_type;
612 mat_(const_cast<DiagonalMatrixType*>(mat))
622 row_ =
row_type(&(mat_->diagonal(i)), i);
628 return mat_==other.mat_;
633 mutable DiagonalMatrixType* mat_;
634 mutable row_type row_;
640 template<
class K,
int n >
643 template<
class DiagonalMatrixType>
693 #ifdef DUNE_FMatrix_WITH_CHECKING
695 DUNE_THROW(FMatrixError,
"index is not contained in pattern");
779 return const_cast<K*
>(
p_);
792 template<
class K,
int n >
795 template<
class DiagonalMatrixType>
834 #ifdef DUNE_FMatrix_WITH_CHECKING
836 DUNE_THROW(FMatrixError,
"index is contained in pattern");
904 template<
class K,
int n>
910 template<
class K,
int n>
916 template<
class K,
int n>
922 template<
class K,
int n>
952 template<
class CW,
class T,
class R>
972 containerWrapper_(containerWrapper),
976 template<
class OtherContainerWrapperIteratorType>
978 containerWrapper_(other.containerWrapper_),
979 position_(other.position_)
983 containerWrapper_(other.containerWrapper_),
984 position_(other.position_)
988 containerWrapper_(other.containerWrapper_),
989 position_(other.position_)
992 template<
class OtherContainerWrapperIteratorType>
995 containerWrapper_ = other.containerWrapper_;
996 position_ = other.position_;
1003 return containerWrapper_.pointer(position_);
1009 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1014 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
1019 return *containerWrapper_.pointer(position_);
1036 return *containerWrapper_.pointer(position_+i);
1041 position_=position_+n;
1044 template<
class OtherContainerWrapperIteratorType>
1045 std::ptrdiff_t
distanceTo(OtherContainerWrapperIteratorType& other)
const
1047 assert(containerWrapper_.identical(other));
1048 return other.position_ - position_;
1053 return containerWrapper_.realIndex(position_);
1057 NonConstCW containerWrapper_;
1063 template<
class M,
class K,
int n>
1067 for(
int i=0; i<n; ++i)
1068 fm[i][i] = s.diagonal()[i];