From b0feb4264935424780520a2f8536269fd1ae6886 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 25 Jul 2016 15:18:55 +0200
Subject: [PATCH] Incorporate matrix transforms from StaticMatrix
---
dune/matrix-vector/CMakeLists.txt | 1 +
dune/matrix-vector/transformmatrix.hh | 358 ++++++++++++++++++++++++++
2 files changed, 359 insertions(+)
create mode 100644 dune/matrix-vector/transformmatrix.hh
diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt
index c3d9227..4d090e0 100644
--- a/dune/matrix-vector/CMakeLists.txt
+++ b/dune/matrix-vector/CMakeLists.txt
@@ -13,4 +13,5 @@ install(FILES
singlenonzerocolumnmatrix.hh
singlenonzerorowmatrix.hh
tranpose.hh
+ transformmatrix.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/matrix-vector)
\ No newline at end of file
diff --git a/dune/matrix-vector/transformmatrix.hh b/dune/matrix-vector/transformmatrix.hh
new file mode 100644
index 0000000..3ce403e
--- /dev/null
+++ b/dune/matrix-vector/transformmatrix.hh
@@ -0,0 +1,358 @@
+#ifndef DUNE_MATRIX_VECTOR_ADDTRANSFORMED_HH
+#define DUNE_MATRIX_VECTOR_ADDTRANSFORMED_HH
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+#include "axpy.hh"
+#include "scalartraits.hh"
+
+namespace Dune {
+namespace MatrixVector {
+
+ // add transformed matrix A += T1^t*B*T2
+ // ******************************************************
+ template <class MatrixA, class MatrixB, class TransformationMatrix,
+ bool AisScalar, bool BisScalar, bool TisScalar>
+ struct TransformMatrixHelper {
+ static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1,
+ const MatrixB& B,
+ const TransformationMatrix& T2) {
+ for (size_t i = 0; i < A.N(); ++i) {
+ typename MatrixA::row_type::Iterator jIt = A[i].begin();
+ typename MatrixA::row_type::Iterator jEnd = A[i].end();
+ for (; jIt != jEnd; ++jIt) {
+ for (size_t k = 0; k < B.N(); ++k) {
+ if (T1.exists(k, i)) {
+ typename MatrixB::row_type::ConstIterator lIt = B[k].begin();
+ typename MatrixB::row_type::ConstIterator lEnd = B[k].end();
+ for (; lIt != lEnd; ++lIt)
+ if (T2.exists(lIt.index(), jIt.index()))
+ *jIt += T1[k][i] * (*lIt) * T2[lIt.index()][jIt.index()];
+ }
+ }
+ }
+ }
+ }
+
+ // static void transformMatrix(MatrixA& A, const
+ // TransformationMatrix& T1, const MatrixB& B, const
+ // TransformationMatrix& T2)
+ // {
+ // A = 0.0;
+ // for (int k=0; k<B.N(); ++k)
+ // {
+ // typename MatrixB::row_type::ConstIterator
+ // lIt=B[k].begin();
+ // typename MatrixB::row_type::ConstIterator
+ // lEnd=B[k].end();
+ // for (; lIt!=lEnd; ++lIt)
+ // {
+ // typename
+ // TransformationMatrix::row_type::ConstIterator
+ // iIt=T1[k].begin();
+ // typename
+ // TransformationMatrix::row_type::ConstIterator
+ // iEnd=T1[k].end();
+ // for (; iIt!=iEnd; ++iIt)
+ // {
+ // typename
+ // TransformationMatrix::row_type::ConstIterator
+ // jIt=T2[lIt.index()].begin();
+ // typename
+ // TransformationMatrix::row_type::ConstIterator
+ // jEnd=T2[lIt.index()].end();
+ // for (; jIt!=jEnd; ++jIt)
+ // A[iIt.index()][jIt.index()] = (*iIt) *
+ // (*lIt) * (*jIt);
+ // }
+ // }
+ // }
+ // }
+ };
+
+ template <class K1, class K2, class K3, int n, int m>
+ struct TransformMatrixHelper<
+ Dune::FieldMatrix<K1, m, m>, Dune::FieldMatrix<K3, n, n>,
+ Dune::FieldMatrix<K2, n, m>, false, false, false> {
+ typedef Dune::FieldMatrix<K1, m, m> MatrixA;
+ typedef Dune::FieldMatrix<K3, n, n> MatrixB;
+ typedef Dune::FieldMatrix<K2, n, m> TransformationMatrix;
+
+ static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1,
+ const MatrixB& B,
+ const TransformationMatrix& T2) {
+ Dune::FieldMatrix<K1, m, n> T1transposedB;
+ T1transposedB = 0;
+ for (size_t i = 0; i < MatrixA::rows; ++i) {
+ for (size_t k = 0; k < MatrixB::rows; ++k) {
+ if (T1[k][i] != 0)
+ for (size_t l = 0; l < MatrixB::cols; ++l)
+ T1transposedB[i][l] += T1[k][i] * B[k][l];
+ }
+ }
+ for (size_t k = 0; k < TransformationMatrix::rows; ++k) {
+ for (size_t l = 0; l < TransformationMatrix::cols; ++l) {
+ if (T2[k][l] != 0)
+ for (size_t i = 0; i < MatrixA::rows; ++i)
+ A[i][l] += T1transposedB[i][k] * T2[k][l];
+ }
+ }
+ }
+
+ // static void transformMatrix(MatrixA& A, const
+ // TransformationMatrix& T1, const MatrixB& B, const
+ // TransformationMatrix& T2)
+ // {
+ // A = 0.0;
+ // for (int k=0; k<B.N(); ++k)
+ // for (int l=0; l<B.M(); ++l)
+ // for (int i=0; i<T1.M(); ++i)
+ // for (int j=0; j<T2.M(); ++j)
+ // A[i][j] = T1[k][i] * B[k][l] * T2[l][j];
+ // }
+ };
+
+ template <class MatrixA, class MatrixB, class ScalarTransform, bool AisScalar,
+ bool BisScalar>
+ struct TransformMatrixHelper<MatrixA, MatrixB, ScalarTransform, AisScalar,
+ BisScalar, true> {
+ static void addTransformedMatrix(MatrixA& A, const ScalarTransform& T1,
+ const MatrixB& B,
+ const ScalarTransform& T2) {
+ addProduct(A, T1 * T2, B);
+ }
+ };
+
+ template <class MatrixA, class ScalarB, class TransformationMatrix,
+ bool AisScalar, bool TisScalar>
+ struct TransformMatrixHelper<MatrixA, ScalarB, TransformationMatrix,
+ AisScalar, true, TisScalar> {
+ static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1,
+ const ScalarB& B,
+ const TransformationMatrix& T2) {
+ for (size_t k = 0; k < TransformationMatrix::rows; ++k) {
+ typename TransformationMatrix::ConstColIterator Skj = T2[k].begin();
+ typename TransformationMatrix::ConstColIterator SkEnd = T2[k].end();
+ for (; Skj != SkEnd; ++Skj) {
+ typename TransformationMatrix::ConstColIterator Tki = T1[k].begin();
+ typename TransformationMatrix::ConstColIterator TkEnd = T1[k].end();
+ for (; Tki != TkEnd; Tki++)
+ if (A.exists(Tki.index(), Skj.index()))
+ A[Tki.index()][Skj.index()] += (*Tki) * B * (*Skj);
+ }
+ }
+ }
+ };
+
+ template <class FieldType, int n, class ScalarB, class TransformationMatrix,
+ bool AisScalar, bool TisScalar>
+ struct TransformMatrixHelper<Dune::ScaledIdentityMatrix<FieldType, n>,
+ ScalarB, TransformationMatrix, AisScalar, true,
+ TisScalar> {
+ typedef Dune::ScaledIdentityMatrix<FieldType, n> MatrixA;
+ static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1,
+ const ScalarB& B,
+ const TransformationMatrix& T2) {
+ TransformMatrixHelper<FieldType, ScalarB,
+ typename TransformationMatrix::field_type, true,
+ true, true>::addTransformedMatrix(A.scalar(),
+ T1.scalar(), B,
+ T2.scalar());
+ }
+ };
+
+ template <class ScalarA, class ScalarB, class ScalarTransform>
+ struct TransformMatrixHelper<ScalarA, ScalarB, ScalarTransform, true, true,
+ true> {
+ static void addTransformedMatrix(ScalarA& A, const ScalarTransform& T1,
+ const ScalarB& B,
+ const ScalarTransform& T2) {
+ A += T1 * B * T2;
+ }
+ };
+
+ template <class MatrixA, typename FieldType, int n,
+ class TransformationMatrix>
+ struct TransformMatrixHelper<MatrixA, Dune::DiagonalMatrix<FieldType, n>,
+ TransformationMatrix, false, false, false> {
+ static void addTransformedMatrix(
+ MatrixA& A, const TransformationMatrix& T1,
+ const Dune::DiagonalMatrix<FieldType, n>& B,
+ const TransformationMatrix& T2) {
+ for (size_t k = 0; k < n; ++k) {
+ typename TransformationMatrix::ConstColIterator Skj = T2[k].begin();
+ typename TransformationMatrix::ConstColIterator SkEnd = T2[k].end();
+ for (; Skj != SkEnd; ++Skj) {
+ typename TransformationMatrix::ConstColIterator Tki = T1[k].begin();
+ typename TransformationMatrix::ConstColIterator TkEnd = T1[k].end();
+ for (; Tki != TkEnd; Tki++)
+ if (A.exists(Tki.index(), Skj.index())) {
+ A[Tki.index()][Skj.index()] += (*Tki) * B.diagonal(k) * (*Skj);
+ }
+ }
+ }
+ }
+ };
+
+ template <class MatrixA, typename FieldType, int n,
+ class TransformationMatrix>
+ struct TransformMatrixHelper<MatrixA,
+ Dune::ScaledIdentityMatrix<FieldType, n>,
+ TransformationMatrix, false, false, false> {
+ static void addTransformedMatrix(
+ MatrixA& A, const TransformationMatrix& T1,
+ const Dune::ScaledIdentityMatrix<FieldType, n>& B,
+ const TransformationMatrix& T2) {
+ TransformMatrixHelper<MatrixA, FieldType, TransformationMatrix, false,
+ true, false>::addTransformedMatrix(A, T1,
+ B.scalar(), T2);
+ }
+ };
+
+ template <typename FieldType, int n, class TransformationMatrix>
+ struct TransformMatrixHelper<Dune::DiagonalMatrix<FieldType, n>,
+ Dune::DiagonalMatrix<FieldType, n>,
+ TransformationMatrix, false, false, false> {
+ static void addTransformedMatrix(
+ Dune::DiagonalMatrix<FieldType, n>& A, const TransformationMatrix& T1,
+ const Dune::DiagonalMatrix<FieldType, n>& B,
+ const TransformationMatrix& T2) {
+ for (int k = 0; k < n; k++) {
+ typename TransformationMatrix::ConstColIterator Tki = T1[k].begin();
+ typename TransformationMatrix::ConstColIterator TkEnd = T1[k].end();
+ for (; Tki != TkEnd; ++Tki)
+ A.diagonal(Tki.index()) +=
+ (*Tki) * B.diagonal(k) * T2[k][Tki.index()];
+ }
+ }
+ };
+
+ template <typename FieldType, int n, class TransformationMatrix>
+ struct TransformMatrixHelper<Dune::ScaledIdentityMatrix<FieldType, n>,
+ Dune::ScaledIdentityMatrix<FieldType, n>,
+ TransformationMatrix, false, false, false> {
+ static void addTransformedMatrix(
+ Dune::ScaledIdentityMatrix<FieldType, n>& A,
+ const TransformationMatrix& T1,
+ const Dune::ScaledIdentityMatrix<FieldType, n>& B,
+ const TransformationMatrix& T2) {
+ for (int k = 0; k < n; k++)
+ A.scalar() += T1[k][0] * B.scalar() * T2[k][0];
+ }
+ };
+
+ template <typename FieldType, int n>
+ struct TransformMatrixHelper<Dune::ScaledIdentityMatrix<FieldType, n>,
+ Dune::ScaledIdentityMatrix<FieldType, n>,
+ Dune::ScaledIdentityMatrix<FieldType, n>, false,
+ false, false> {
+ static void addTransformedMatrix(
+ Dune::ScaledIdentityMatrix<FieldType, n>& A,
+ const Dune::ScaledIdentityMatrix<FieldType, n>& T1,
+ const Dune::ScaledIdentityMatrix<FieldType, n>& B,
+ const Dune::ScaledIdentityMatrix<FieldType, n>& T2) {
+ TransformMatrixHelper<FieldType, FieldType, FieldType, true, true,
+ true>::addTransformedMatrix(A.scalar(), T1.scalar(),
+ B.scalar(),
+ T2.scalar());
+ }
+ };
+
+ template <class MatrixA, class MatrixB, class TransformationMatrix>
+ void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1,
+ const MatrixB& B, const TransformationMatrix& T2) {
+ TransformMatrixHelper<
+ MatrixA, MatrixB, TransformationMatrix, ScalarTraits<MatrixA>::isScalar,
+ ScalarTraits<MatrixB>::isScalar,
+ ScalarTraits<TransformationMatrix>::isScalar>::addTransformedMatrix(A,
+ T1,
+ B,
+ T2);
+ }
+
+ template <class MatrixA, class MatrixB, class TransformationMatrix>
+ void transformMatrix(MatrixA& A, const TransformationMatrix& T1,
+ const MatrixB& B, const TransformationMatrix& T2) {
+ A = 0;
+ TransformMatrixHelper<
+ MatrixA, MatrixB, TransformationMatrix, ScalarTraits<MatrixA>::isScalar,
+ ScalarTraits<MatrixB>::isScalar,
+ ScalarTraits<TransformationMatrix>::isScalar>::addTransformedMatrix(A,
+ T1,
+ B,
+ T2);
+ }
+
+ template <class MatrixBlockA, class MatrixB, class TransformationMatrix>
+ static void transformMatrix(Dune::BCRSMatrix<MatrixBlockA>& A,
+ const TransformationMatrix& T1, const MatrixB& B,
+ const TransformationMatrix& T2) {
+ transformMatrixPattern(A, T1, B, T2);
+ A = 0.0;
+ for (int k = 0; k < B.N(); ++k) {
+ typename MatrixB::row_type::ConstIterator BklIt = B[k].begin();
+ typename MatrixB::row_type::ConstIterator BklEnd = B[k].end();
+ for (; BklIt != BklEnd; ++BklIt) {
+ int l = BklIt.index();
+
+ typename TransformationMatrix::row_type::ConstIterator T1kiIt =
+ T1[k].begin();
+ typename TransformationMatrix::row_type::ConstIterator T1kiEnd =
+ T1[k].end();
+ for (; T1kiIt != T1kiEnd; ++T1kiIt) {
+ int i = T1kiIt.index();
+
+ typename TransformationMatrix::row_type::ConstIterator T2ljIt =
+ T2[l].begin();
+ typename TransformationMatrix::row_type::ConstIterator T2ljEnd =
+ T2[l].end();
+ for (; T2ljIt != T2ljEnd; ++T2ljIt) {
+ int j = T2ljIt.index();
+ MatrixBlockA Aij;
+ transformMatrix(Aij, *T1kiIt, *BklIt, *T2ljIt);
+ A[i][j] += Aij;
+ }
+ }
+ }
+ }
+ }
+
+ template <class MatrixBlockA, class MatrixB, class TransformationMatrix>
+ static void transformMatrixPattern(Dune::BCRSMatrix<MatrixBlockA>& A,
+ const TransformationMatrix& T1,
+ const MatrixB& B,
+ const TransformationMatrix& T2) {
+ Dune::MatrixIndexSet indices(T1.M(), T2.M());
+ for (int k = 0; k < B.N(); ++k) {
+ typename MatrixB::row_type::ConstIterator BklIt = B[k].begin();
+ typename MatrixB::row_type::ConstIterator BklEnd = B[k].end();
+ for (; BklIt != BklEnd; ++BklIt) {
+ int l = BklIt.index();
+
+ typename TransformationMatrix::row_type::ConstIterator T1kiIt =
+ T1[k].begin();
+ typename TransformationMatrix::row_type::ConstIterator T1kiEnd =
+ T1[k].end();
+ for (; T1kiIt != T1kiEnd; ++T1kiIt) {
+ int i = T1kiIt.index();
+
+ typename TransformationMatrix::row_type::ConstIterator T2ljIt =
+ T2[l].begin();
+ typename TransformationMatrix::row_type::ConstIterator T2ljEnd =
+ T2[l].end();
+ for (; T2ljIt != T2ljEnd; ++T2ljIt) {
+ int j = T2ljIt.index();
+ indices.add(i, j);
+ }
+ }
+ }
+ }
+ indices.exportIdx(A);
+ }
+}
+}
+
+#endif
--
GitLab