diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh index 041baf14203af30ff95c742490c86a72c233d429..3abe108c93682f487f9a4818ec27cbd86220af77 100644 --- a/dune/solvers/common/staticmatrixtools.hh +++ b/dune/solvers/common/staticmatrixtools.hh @@ -24,6 +24,18 @@ class StaticMatrix typedef Matrix Type; }; + template<typename FieldType, int n> + struct Promote<Dune::FieldMatrix<FieldType,n,n>, Dune::DiagonalMatrix<FieldType,n> > + { + typedef Dune::FieldMatrix<FieldType,n,n> Type; + }; + + template<typename FieldType, int n> + struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::FieldMatrix<FieldType,n,n> > + { + typedef Dune::FieldMatrix<FieldType,n,n> Type; + }; + template<typename FieldType, int n> struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n> > { @@ -47,7 +59,14 @@ class StaticMatrix } template <typename FieldType, int n> - static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a) +// static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a) + // the commented line should NOT be used as it may lead to ambiguities in which case the general method above will be used. + // an example is the line taken from massassembler.hh (line 83 at the time): + // StaticMatrix::addToDiagonal(localMatrix[i][i], values[i] * zi); + // where values[i] is of type FielVector<FieldType,1> and zi is a double. This call then does not exactly fit this specialization (without the implicit cast of FV<FT,1>) + // and hence some wild template voodoo decides which of the templates is to be taken - in this case with a gcc 4.4.5 that was the one above leading to a wrong factor n + // in the diagonal value + static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a) { x.scalar() += a; } @@ -58,7 +77,7 @@ class StaticMatrix template<class MatrixA, class MatrixB> static void axpy(MatrixA& A, const typename MatrixA::field_type a, const MatrixB& B) { - for(int i=0; i<B.N(); ++i) + for(size_t i=0; i<B.N(); ++i) { typename MatrixB::row_type::ConstIterator it=B[i].begin(); typename MatrixB::row_type::ConstIterator end=B[i].end(); @@ -76,13 +95,13 @@ class StaticMatrix // } template <typename FieldType, int n, class MatrixB> - static void axpy(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a, const MatrixB& y) + static void axpy(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a, const MatrixB& y) { x.scalar() += a * y[0][0]; } template <typename FieldType, int n, class MatrixB> - static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const FieldType a, const MatrixB& y) + static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const typename Dune::DiagonalMatrix<FieldType,n>::field_type a, const MatrixB& y) { for(int i=0; i<n; ++i) x.diagonal()[i] += a * y[i][i]; @@ -95,17 +114,20 @@ class StaticMatrix } + + + // add transformed matrix X = A^T*Y*B ****************************************************** template<class MatrixA, class MatrixB, class TransformationMatrix> static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2) { - for (int i=0; i<A.N(); ++i) + for (size_t i=0; i<A.N(); ++i) { typename MatrixA::row_type::ConstIterator jIt=A[i].begin(); typename MatrixA::row_type::ConstIterator jEnd=A[i].end(); for(; jIt!=jEnd; ++jIt) { - for (int k=0; k<B.N(); ++k) + for (size_t k=0; k<B.N(); ++k) { typename MatrixB::row_type::ConstIterator lIt=B[k].begin(); typename MatrixB::row_type::ConstIterator lEnd=B[k].end();