Code owners
Assign users and groups as approvers for specific file changes. Learn more.
transformmatrix.hh 13.76 KiB
#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 TransformationMatrix1, class MatrixB,
class TransformationMatrix2, bool AisScalar, bool T1isScalar,
bool BisScalar, bool T2isScalar>
struct TransformMatrixHelper {
static void addTransformedMatrix(MatrixA& A,
const TransformationMatrix1& T1,
const MatrixB& B,
const TransformationMatrix2& T2) {
for (size_t i = 0; i < A.N(); ++i)
for (auto jIt = A[i].begin(); jIt != A[i].end(); ++jIt)
for (size_t k = 0; k < B.N(); ++k)
if (T1.exists(k, i))
for (auto lIt = B[k].begin(); lIt != B[k].end(); ++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<K2, n, m>,
Dune::FieldMatrix<K3, n, n>, Dune::FieldMatrix<K2, n, m>, false, 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 ScalarTransform1, class MatrixB,
class ScalarTransform2, bool AisScalar, bool BisScalar>
struct TransformMatrixHelper<MatrixA, ScalarTransform1, MatrixB,
ScalarTransform2, AisScalar, true, BisScalar,
true> {
static void addTransformedMatrix(MatrixA& A, const ScalarTransform1& T1,
const MatrixB& B,
const ScalarTransform2& T2) {
addProduct(A, T1 * T2, B);
}
};
template <class MatrixA, class TransformationMatrix1, class ScalarB,
class TransformationMatrix2, bool AisScalar, bool T1isScalar,
bool T2isScalar>
struct TransformMatrixHelper<MatrixA, TransformationMatrix1, ScalarB,
TransformationMatrix2, AisScalar, T1isScalar,
true, T2isScalar> {
static void addTransformedMatrix(MatrixA& A,
const TransformationMatrix1& T1,
const ScalarB& B,
const TransformationMatrix2& T2) {
assert(TransformationMatrix1::rows == TransformationMatrix2::rows);
for (size_t k = 0; k < TransformationMatrix1::rows; ++k)
for (auto Skj = T2[k].begin(); Skj != T2[k].end(); ++Skj)
for (auto Tki = T1[k].begin(); Tki != T1[k].end(); Tki++)
if (A.exists(Tki.index(), Skj.index()))
A[Tki.index()][Skj.index()] += (*Tki) * B * (*Skj);
}
};
template <class FieldType, int n, class TransformationMatrix1, class ScalarB, class TransformationMatrix2, bool AisScalar, bool T1isScalar,
bool T2isScalar>
struct TransformMatrixHelper<
Dune::ScaledIdentityMatrix<FieldType, n>, TransformationMatrix1, ScalarB,
TransformationMatrix2, AisScalar, T1isScalar, true, T2isScalar> {
typedef Dune::ScaledIdentityMatrix<FieldType, n> MatrixA;
static void addTransformedMatrix(MatrixA& A,
const TransformationMatrix1& T1,
const ScalarB& B,
const TransformationMatrix2& T2) {
TransformMatrixHelper<
FieldType, typename TransformationMatrix1::field_type, ScalarB,
typename TransformationMatrix2::field_type, true, true, true,
true>::addTransformedMatrix(A.scalar(), T1.scalar(), B, T2.scalar());
}
};
template <class ScalarA, class ScalarTransform1, class ScalarB,
class ScalarTransform2>
struct TransformMatrixHelper<ScalarA, ScalarTransform1, ScalarB,
ScalarTransform2, true, true, true, true> {
static void addTransformedMatrix(ScalarA& A, const ScalarTransform1& T1,
const ScalarB& B,
const ScalarTransform2& T2) {
A += T1 * B * T2;
}
};
template <class MatrixA, class TransformationMatrix1, typename FieldType,
int n, class TransformationMatrix2>
struct TransformMatrixHelper<
MatrixA, TransformationMatrix1, Dune::DiagonalMatrix<FieldType, n>,
TransformationMatrix2, false, false, false, false> {
static void addTransformedMatrix(
MatrixA& A, const TransformationMatrix1& T1,
const Dune::DiagonalMatrix<FieldType, n>& B,
const TransformationMatrix2& T2) {
for (size_t k = 0; k < n; ++k)
for (auto Skj = T2[k].begin(); Skj != T2[k].end(); ++Skj)
for (auto Tki = T1[k].begin(); Tki != T1[k].end(); Tki++)
if (A.exists(Tki.index(), Skj.index()))
A[Tki.index()][Skj.index()] += (*Tki) * B.diagonal(k) * (*Skj);
}
};
template <class MatrixA, class TransformationMatrix1, typename FieldType,
int n, class TransformationMatrix2>
struct TransformMatrixHelper<
MatrixA, TransformationMatrix1, Dune::ScaledIdentityMatrix<FieldType, n>,
TransformationMatrix2, false, false, false, false> {
static void addTransformedMatrix(
MatrixA& A, const TransformationMatrix1& T1,
const Dune::ScaledIdentityMatrix<FieldType, n>& B,
const TransformationMatrix2& T2) {
TransformMatrixHelper<MatrixA, TransformationMatrix1, FieldType,
TransformationMatrix2, false, false, true,
false>::addTransformedMatrix(A, T1, B.scalar(), T2);
}
};
template <typename FieldType, int n, class TransformationMatrix1,
class TransformationMatrix2>
struct TransformMatrixHelper<
Dune::DiagonalMatrix<FieldType, n>, TransformationMatrix1,
Dune::DiagonalMatrix<FieldType, n>, TransformationMatrix2, false, false,
false, false> {
static void addTransformedMatrix(
Dune::DiagonalMatrix<FieldType, n>& A, const TransformationMatrix1& T1,
const Dune::DiagonalMatrix<FieldType, n>& B,
const TransformationMatrix2& T2) { for (size_t k = 0; k < n; k++)
for (auto Tki = T1[k].begin(); Tki != T1[k].end(); ++Tki)
A.diagonal(Tki.index()) +=
(*Tki) * B.diagonal(k) * T2[k][Tki.index()];
}
};
template <typename FieldType, int n, class TransformationMatrix1,
class TransformationMatrix2>
struct TransformMatrixHelper<
Dune::ScaledIdentityMatrix<FieldType, n>, TransformationMatrix1,
Dune::ScaledIdentityMatrix<FieldType, n>, TransformationMatrix2, false,
false, false, false> {
static void addTransformedMatrix(
Dune::ScaledIdentityMatrix<FieldType, n>& A,
const TransformationMatrix1& T1,
const Dune::ScaledIdentityMatrix<FieldType, n>& B,
const TransformationMatrix2& T2) {
for (size_t 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>,
Dune::ScaledIdentityMatrix<FieldType, n>, false,
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, FieldType, true, true, true,
true>::addTransformedMatrix(A.scalar(), T1.scalar(), B.scalar(),
T2.scalar());
}
};
template <class MatrixA, class TransformationMatrix1, class MatrixB,
class TransformationMatrix2>
void addTransformedMatrix(MatrixA& A, const TransformationMatrix1& T1,
const MatrixB& B, const TransformationMatrix2& T2) {
TransformMatrixHelper<
MatrixA, TransformationMatrix1, MatrixB, TransformationMatrix2,
ScalarTraits<MatrixA>::isScalar,
ScalarTraits<TransformationMatrix1>::isScalar,
ScalarTraits<MatrixB>::isScalar,
ScalarTraits<
TransformationMatrix2>::isScalar>::addTransformedMatrix(A, T1, B,
T2);
}
template <class MatrixA, class TransformationMatrix1, class MatrixB,
class TransformationMatrix2>
void transformMatrix(MatrixA& A, const TransformationMatrix1& T1,
const MatrixB& B, const TransformationMatrix2& T2) {
A = 0;
TransformMatrixHelper<
MatrixA, TransformationMatrix1, MatrixB, TransformationMatrix2,
ScalarTraits<MatrixA>::isScalar,
ScalarTraits<TransformationMatrix1>::isScalar,
ScalarTraits<MatrixB>::isScalar,
ScalarTraits<
TransformationMatrix2>::isScalar>::addTransformedMatrix(A, T1, B,
T2);
}
template <class MatrixBlockA, class TransformationMatrix1, class MatrixB,
class TransformationMatrix2>
static void transformMatrix(Dune::BCRSMatrix<MatrixBlockA>& A,
const TransformationMatrix1& T1, const MatrixB& B,
const TransformationMatrix2& T2) {
transformMatrixPattern(A, T1, B, T2);
A = 0.0;
for (size_t k = 0; k < B.N(); ++k)
for (auto BkIt = B[k].begin(); BkIt != B[k].end(); ++BkIt) {
size_t l = BkIt.index();
for (auto T1kIt = T1[k].begin(); T1kIt != T1[k].end(); ++T1kIt)
for (auto T2lIt = T2[l].begin(); T2lIt != T2[l].end(); ++T2lIt) {
MatrixBlockA Aij;
transformMatrix(Aij, *T1kIt, *BkIt, *T2lIt);
A[T1kIt.index()][T2lIt.index()] += Aij;
}
}
}
template <class MatrixBlockA, class TransformationMatrix1, class MatrixB,
class TransformationMatrix2>
static void transformMatrixPattern(Dune::BCRSMatrix<MatrixBlockA>& A,
const TransformationMatrix1& T1,
const MatrixB& B,
const TransformationMatrix2& T2) {
Dune::MatrixIndexSet indices(T1.M(), T2.M());
for (size_t k = 0; k < B.N(); ++k)
for (auto BIt = B[k].begin(); BIt != B[k].end(); ++BIt) {
size_t l = BIt.index();
for (auto T1kIt = T1[k].begin(); T1kIt != T1[k].end(); ++T1kIt)
for (auto T2lIt = T2[l].begin(); T2lIt != T2[l].end(); ++T2lIt)
indices.add(T1kIt.index(), T2lIt.index());
}
indices.exportIdx(A);
}
}
}
#endif