From 273fa630d87c5fe296efd24221c993eda7fb143f Mon Sep 17 00:00:00 2001 From: Elias Pipping <elias.pipping@fu-berlin.de> Date: Mon, 25 Jul 2016 12:53:34 +0200 Subject: [PATCH] Incorporate transpose from Arithmetic --- dune/matrix-vector/CMakeLists.txt | 1 + dune/matrix-vector/transpose.hh | 101 ++++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+) create mode 100644 dune/matrix-vector/transpose.hh diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index 3ffe97f..8d64790 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 0000000..627bcee --- /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 -- GitLab