diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index 3ffe97f4c51f8470b6c976bdac8efa3a58dade61..8d647905e97d99835ec1506bf9aa31c616ccb3f7 100644 --- a/dune/matrix-vector/CMakeLists.txt +++ b/dune/matrix-vector/CMakeLists.txt @@ -7,4 +7,5 @@ install(FILES scalartraits.hh singlenonzerocolumnmatrix.hh singlenonzerorowmatrix.hh + tranpose.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/matrix-vector) \ No newline at end of file diff --git a/dune/matrix-vector/transpose.hh b/dune/matrix-vector/transpose.hh new file mode 100644 index 0000000000000000000000000000000000000000..627bceea1f4e04bf37130a3c2fae264cdde3f4ff --- /dev/null +++ b/dune/matrix-vector/transpose.hh @@ -0,0 +1,101 @@ +#ifndef DUNE_MATRIX_VECTOR_TRANSPOSE_HH +#define DUNE_MATRIX_VECTOR_TRANSPOSE_HH + +#include <dune/common/diagonalmatrix.hh> +#include <dune/common/fmatrix.hh> +#include <dune/istl/bcrsmatrix.hh> +#include <dune/istl/matrixindexset.hh> +#include <dune/istl/scaledidmatrix.hh> + +#include <dune/matrix-vector/axpy.hh> + +namespace Dune { + namespace MatrixVector { + + template <class A> + struct TransposeHelper; + + template <class MatrixType> + using Transposed = typename TransposeHelper<MatrixType>::TransposedType; + + template <class A> + struct TransposeHelper { + typedef A TransposedType; + + static void transpose(const A& a, TransposedType& aT) { + DUNE_UNUSED_PARAMETER(a); + DUNE_UNUSED_PARAMETER(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>> { + typedef Dune::FieldMatrix<T, n, m> MatrixType; + typedef Dune::FieldMatrix<T, m, n> TransposedType; + + static void transpose(const MatrixType& a, TransposedType& aT) { + for (int row = 0; row < m; ++row) + for (int col = 0; col < n; ++col) + aT[row][col] = a[col][row]; + } + }; + + template <class T, int n> + struct TransposeHelper<Dune::DiagonalMatrix<T, n>> { + typedef Dune::DiagonalMatrix<T, n> MatrixType; + typedef Dune::DiagonalMatrix<T, n> TransposedType; + + static void transpose(const MatrixType& a, TransposedType& aT) { aT = a; } + }; + + template <class T, int n> + struct TransposeHelper<Dune::ScaledIdentityMatrix<T, n>> { + typedef Dune::ScaledIdentityMatrix<T, n> MatrixType; + typedef Dune::ScaledIdentityMatrix<T, n> TransposedType; + + static void transpose(const MatrixType& a, TransposedType& aT) { aT = a; } + }; + + //! Specialization for Dune::BCRSMatrix Type + template <class A> + struct TransposeHelper<Dune::BCRSMatrix<A>> { + typedef Dune::BCRSMatrix<Transposed<A>> TransposedType; + + static void transpose(const Dune::BCRSMatrix<A>& a, TransposedType& aT) { + Dune::MatrixIndexSet idxSetaT(a.M(), a.N()); + + typedef typename Dune::BCRSMatrix<A>::ConstColIterator ColIterator; + + // add indices into transposed matrix + for (size_t 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 (size_t row = 0; row < a.N(); ++row) { + ColIterator col = a[row].begin(); + ColIterator end = a[row].end(); + + for (; col != end; ++col) + TransposeHelper<A>::transpose(*col, aT[col.index()][row]); + } + } + }; + + //! Compute the transposed of a matrix + template <class A> + void transpose(const A& a, Transposed<A>& aT) { + TransposeHelper<A>::transpose(a, aT); + } + } +} + +#endif