diff --git a/dune/solvers/common/arithmetic.hh b/dune/solvers/common/arithmetic.hh
index 6461e881b246f4a7465e8ec3d73600672cf5ab4a..618c65ebb97da6d7c080555bb34fb9413f0ab8a0 100644
--- a/dune/solvers/common/arithmetic.hh
+++ b/dune/solvers/common/arithmetic.hh
@@ -6,6 +6,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/diagonalmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/common/typetraits.hh>
@@ -20,7 +21,7 @@ namespace Arithmetic
/** \brief Class to identify scalar types
*
- * Specialize this class for all type that can be used
+ * Specialize this class for all types that can be used
* like scalar quantities.
*/
template<class T>
@@ -68,7 +69,7 @@ namespace Arithmetic
/** \brief Class to identify matrix types and extract information
*
- * Specialize this class for all type that can be used like a matrix.
+ * Specialize this class for all types that can be used like a matrix.
*/
template<class T>
struct MatrixTraits
@@ -102,6 +103,12 @@ namespace Arithmetic
enum { cols=n};
};
+ template<class T>
+ struct MatrixTraits<Dune::BCRSMatrix<T> >
+ {
+ enum { isMatrix=true };
+ };
+
/** \brief Internal helper class for product operations
*
@@ -317,6 +324,84 @@ namespace Arithmetic
}
};
+ //! Helper class for computing the cross product
+ template <class T, int n>
+ struct CrossProductHelper
+ {
+ static Dune::FieldVector<T,n> crossProduct(const Dune::FieldVector<T,n>& a, const Dune::FieldVector<T,n>& b)
+ {
+ DUNE_THROW(Dune::Exception, "You can only call crossProduct with dim==3");
+ }
+ };
+
+ //! Specialisation for n=3
+ template <class T>
+ struct CrossProductHelper<T,3>
+ {
+ static Dune::FieldVector<T,3> crossProduct(const Dune::FieldVector<T,3>& a, const Dune::FieldVector<T,3>& b)
+ {
+ Dune::FieldVector<T,3> r;
+ r[0] = a[1]*b[2] - a[2]*b[1];
+ r[1] = a[2]*b[0] - a[0]*b[2];
+ r[2] = a[0]*b[1] - a[1]*b[0];
+ return r;
+ }
+ };
+
+ template<class A, class TransposedA>
+ struct TransposeHelper
+ {
+ static void transpose(const A& a, TransposedA& aT)
+ {
+ DUNE_THROW(Dune::Exception, "Not implemented for general matrix types!");
+ }
+ };
+
+ //! Specialization for Dune::FieldMatrix
+ template<class T, int n, int m>
+ struct TransposeHelper<Dune::FieldMatrix<T,n,m>, Dune::FieldMatrix<T,m,n> >
+ {
+ static void transpose(const Dune::FieldMatrix<T,n,m>& a, Dune::FieldMatrix<T,m,n>& aT)
+ {
+ for (int row = 0; row < m; ++row)
+ for (int col = 0 ; col < n; ++col)
+ aT[row][col] = a[col][row];
+ }
+ };
+
+ //! Specialization for Dune::BCRSMatrix Type
+ template<class A, class TransposedA>
+ struct TransposeHelper<Dune::BCRSMatrix<A>, Dune::BCRSMatrix<TransposedA> >
+ {
+ static void transpose(const Dune::BCRSMatrix<A>& a, Dune::BCRSMatrix<TransposedA>& aT)
+ {
+ Dune::MatrixIndexSet idxSetaT(a.M(),a.N());
+
+ typedef typename Dune::BCRSMatrix<A>::ConstColIterator ColIterator;
+
+ // add indices into transposed matrix
+ for (int row = 0; row < a.N(); ++row)
+ {
+ ColIterator col = a[row].begin();
+ ColIterator end = a[row].end();
+
+ for ( ; col != end; ++col)
+ idxSetaT.add(col.index(),row);
+ }
+
+ idxSetaT.exportIdx(aT);
+
+ for (int row = 0; row < a.N(); ++row)
+ {
+ ColIterator col = a[row].begin();
+ ColIterator end = a[row].end();
+
+ for ( ; col != end; ++col)
+ TransposeHelper<A,TransposedA>::transpose(*col, aT[col.index()][row]);
+ }
+ }
+ };
+
/** \brief Add a product to some matrix or vector
*
* This function computes a+=b*c.
@@ -350,6 +435,19 @@ namespace Arithmetic
{
ProductHelper<A,B,C,ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar, ScalarTraits<C>::isScalar>::subtractProduct(a,b,c);
}
+ //! Compute the cross product of two vectors. Only works for n==3
+ template<class T, int n>
+ Dune::FieldVector<T,n> crossProduct(const Dune::FieldVector<T,n>& a, const Dune::FieldVector<T,n>& b)
+ {
+ return CrossProductHelper<T,n>::crossProduct(a,b);
+ }
+
+ //! Compute the transposed of a matrix
+ template <class A, class TransposedA>
+ void transpose(const A& a, TransposedA& aT)
+ {
+ TransposeHelper<A,TransposedA>::transpose(a,aT);
+ }
/** \brief Add a scaled product to some matrix or vector
*
@@ -387,12 +485,17 @@ namespace Arithmetic
ScaledProductHelper<A,B,C,D,ScalarTraits<A>::isScalar, ScalarTraits<C>::isScalar, ScalarTraits<D>::isScalar>::subtractProduct(a,b,c,d);
}
- namespace {
- template <class MatrixType, class VectorType>
- typename VectorType::field_type
- AxyWithTemporary(const MatrixType &A,
- const VectorType &x,
- const VectorType &y)
+ /** \brief Internal helper class for Matrix operations
+ *
+ */
+ template<class OperatorType, bool isMatrix>
+ struct OperatorHelper
+ {
+ template <class VectorType>
+ static typename VectorType::field_type
+ Axy(const OperatorType &A,
+ const VectorType &x,
+ const VectorType &y)
{
VectorType tmp(y.size());
tmp = 0.0;
@@ -400,11 +503,25 @@ namespace Arithmetic
return tmp * y;
}
- template <class MatrixType, class VectorType>
- typename VectorType::field_type
- AxyWithoutTemporary(const MatrixType &A,
- const VectorType &x,
- const VectorType &y)
+ template <class VectorType>
+ static typename VectorType::field_type
+ bmAxy(const OperatorType &A, const VectorType &b,
+ const VectorType &x, const VectorType &y)
+ {
+ VectorType tmp = b;
+ subtractProduct(tmp, A, x);
+ return tmp * y;
+ }
+ };
+
+ template<class MatrixType>
+ struct OperatorHelper<MatrixType, true>
+ {
+ template <class VectorType>
+ static typename VectorType::field_type
+ Axy(const MatrixType &A,
+ const VectorType &x,
+ const VectorType &y)
{
assert(x.N() == A.M());
assert(y.N() == A.N());
@@ -423,53 +540,11 @@ namespace Arithmetic
}
return outer;
}
- }
-
- //! Compute \f$(Ax,y)\f$
- template <class MatrixType, class VectorType>
- typename VectorType::field_type
- Axy(const MatrixType &A,
- const VectorType &x,
- const VectorType &y)
- {
- return AxyWithTemporary(A, x, y);
- }
-
- //! Compute \f$(Ax,y)\f$
- template <int n, typename field_type, class VectorType>
- typename VectorType::field_type
- Axy(const Dune::FieldMatrix<field_type, n, n> &A,
- const VectorType &x,
- const VectorType &y)
- {
- return AxyWithoutTemporary(A, x, y);
- }
-
- //! Compute \f$(Ax,y)\f$
- template <typename block_type, class VectorType>
- typename VectorType::field_type
- Axy(const Dune::BCRSMatrix<block_type> &A,
- const VectorType &x,
- const VectorType &y)
- {
- return AxyWithoutTemporary(A, x, y);
- }
-
- namespace {
- template <class MatrixType, class VectorType>
- typename VectorType::field_type
- bmAxyWithTemporary(const MatrixType &A, const VectorType &b,
- const VectorType &x, const VectorType &y)
- {
- VectorType tmp = b;
- subtractProduct(tmp, A, x);
- return tmp * y;
- }
- template <class MatrixType, class VectorType>
- typename VectorType::field_type
- bmAxyWithoutTemporary(const MatrixType &A, const VectorType &b,
- const VectorType &x, const VectorType &y)
+ template <class VectorType>
+ static typename VectorType::field_type
+ bmAxy(const MatrixType &A, const VectorType &b,
+ const VectorType &x, const VectorType &y)
{
assert(x.N() == A.M());
assert(y.N() == A.N());
@@ -489,33 +564,29 @@ namespace Arithmetic
}
return outer;
}
- }
-
- //! Compute \f$(b-Ax,y)\f$
- template <class MatrixType, class VectorType>
- typename VectorType::field_type
- bmAxy(const MatrixType &A, const VectorType &b,
- const VectorType &x, const VectorType &y)
- {
- return bmAxyWithTemporary(A, b, x, y);
- }
+ };
- //! Compute \f$(b-Ax,y)\f$
- template <int n, typename field_type, class VectorType>
+ //! Compute \f$(Ax,y)\f$
+ template <class OperatorType, class VectorType>
typename VectorType::field_type
- bmAxy(const Dune::FieldMatrix<field_type, n, n> &A, const VectorType &b,
- const VectorType &x, const VectorType &y)
+ Axy(const OperatorType &A,
+ const VectorType &x,
+ const VectorType &y)
{
- return bmAxyWithoutTemporary(A, b, x, y);
+ return OperatorHelper<OperatorType,
+ MatrixTraits<OperatorType>::isMatrix>::Axy(A, x, y);
}
//! Compute \f$(b-Ax,y)\f$
- template <typename block_type, class VectorType>
+ template <class OperatorType, class VectorType>
typename VectorType::field_type
- bmAxy(const Dune::BCRSMatrix<block_type> &A, const VectorType &b,
- const VectorType &x, const VectorType &y)
+ bmAxy(const OperatorType &A,
+ const VectorType &b,
+ const VectorType &x,
+ const VectorType &y)
{
- return bmAxyWithoutTemporary(A, b, x, y);
+ return OperatorHelper<OperatorType,
+ MatrixTraits<OperatorType>::isMatrix>::bmAxy(A, b, x, y);
}
};