3 #ifndef DUNE_SCALED_IDENTITY_MATRIX_HH
4 #define DUNE_SCALED_IDENTITY_MATRIX_HH
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/fmatrix.hh>
18 #include <dune/common/diagonalmatrix.hh>
19 #include <dune/common/ftraits.hh>
26 template<
class K,
int n>
29 typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType;
84 return (
this==&other);
89 typedef ContainerWrapperIterator<const WrapperType, reference, reference>
Iterator;
100 return Iterator(WrapperType(
this),0);
106 return Iterator(WrapperType(
this),n);
113 return Iterator(WrapperType(
this),n-1);
120 return Iterator(WrapperType(
this),-1);
125 typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference>
ConstIterator;
207 return p_==other.
scalar();
213 return p_!=other.
scalar();
219 template<
class X,
class Y>
220 void mv (
const X& x, Y& y)
const
222 #ifdef DUNE_FMatrix_WITH_CHECKING
223 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
224 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
231 template<
class X,
class Y>
232 void mtv (
const X& x, Y& y)
const
238 template<
class X,
class Y>
239 void umv (
const X& x, Y& y)
const
241 #ifdef DUNE_FMatrix_WITH_CHECKING
242 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
243 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
250 template<
class X,
class Y>
251 void umtv (
const X& x, Y& y)
const
253 #ifdef DUNE_FMatrix_WITH_CHECKING
254 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
255 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
262 template<
class X,
class Y>
263 void umhv (
const X& x, Y& y)
const
265 #ifdef DUNE_FMatrix_WITH_CHECKING
266 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
267 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
270 y[i] += conjugateComplex(p_)*x[i];
274 template<
class X,
class Y>
275 void mmv (
const X& x, Y& y)
const
277 #ifdef DUNE_FMatrix_WITH_CHECKING
278 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
279 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
286 template<
class X,
class Y>
287 void mmtv (
const X& x, Y& y)
const
289 #ifdef DUNE_FMatrix_WITH_CHECKING
290 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
291 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
298 template<
class X,
class Y>
299 void mmhv (
const X& x, Y& y)
const
301 #ifdef DUNE_FMatrix_WITH_CHECKING
302 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
303 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
306 y[i] -= conjugateComplex(p_)*x[i];
310 template<
class X,
class Y>
311 void usmv (
const K& alpha,
const X& x, Y& y)
const
313 #ifdef DUNE_FMatrix_WITH_CHECKING
314 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
315 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
318 y[i] += alpha * p_ * x[i];
322 template<
class X,
class Y>
323 void usmtv (
const K& alpha,
const X& x, Y& y)
const
325 #ifdef DUNE_FMatrix_WITH_CHECKING
326 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
327 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
330 y[i] += alpha * p_ * x[i];
334 template<
class X,
class Y>
335 void usmhv (
const K& alpha,
const X& x, Y& y)
const
337 #ifdef DUNE_FMatrix_WITH_CHECKING
338 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
339 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
342 y[i] += alpha * conjugateComplex(p_) * x[i];
350 return fvmeta::sqrt(n*p_*p_);
368 return fvmeta::absreal(p_);
378 for (
int i=0; i<n; i++)
391 return std::pow(p_,n);
413 #ifdef DUNE_FMatrix_WITH_CHECKING
414 if (i<0 || i>=n) DUNE_THROW(FMatrixError,
"row index out of range");
415 if (j<0 || j>=n) DUNE_THROW(FMatrixError,
"column index out of range");
423 friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
427 s << ((i==j) ? a.p_ : 0) <<
" ";
436 return reference(const_cast<K*>(&p_), i);
477 template<
class M,
class K,
int n>
481 for(
int i=0; i<n; ++i)