diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh
index 1b69bf6f2cdc6367fc6af1ad1df20b918ff281aa..7fcc02e40e8222ca583e9f2b294f09c651f68a09 100644
--- a/dune/solvers/common/staticmatrixtools.hh
+++ b/dune/solvers/common/staticmatrixtools.hh
@@ -14,6 +14,7 @@
#include <dune/matrix-vector/ldlt.hh>
#include <dune/matrix-vector/promote.hh>
#include <dune/matrix-vector/scalartraits.hh>
+#include <dune/matrix-vector/transformmatrix.hh>
#include "arithmetic.hh"
@@ -34,291 +35,22 @@ namespace StaticMatrix
}
// 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)
- {
- Dune::MatrixVector::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,
- Dune::MatrixVector::ScalarTraits<MatrixA>::isScalar, Dune::MatrixVector::ScalarTraits<MatrixB>::isScalar, Dune::MatrixVector::ScalarTraits<TransformationMatrix>::isScalar>
- ::addTransformedMatrix(A,T1,B,T2);
+ Dune::MatrixVector::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,
- Dune::MatrixVector::ScalarTraits<MatrixA>::isScalar, Dune::MatrixVector::ScalarTraits<MatrixB>::isScalar, Dune::MatrixVector::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;
- }
- }
- }
- }
+ Dune::MatrixVector::transformMatrix(A,T1,B,T2);
}
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);
+ Dune::MatrixVector::transformMatrixPattern(A,T1,B,T2);
}
diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index b4a6ca4284c98a5bf5fb8af7ea7c3c5dbd34a00c..60313b54ef929c02b3a1e0b2dcddce0e3030439f 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh
@@ -13,9 +13,7 @@
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dune/matrix-vector/axpy.hh>
-
-#include "dune/solvers/common/staticmatrixtools.hh"
-
+#include <dune/matrix-vector/transformmatrix.hh>
/** \brief Restriction and prolongation operator for standard multigrid
*
@@ -700,7 +698,7 @@ public:
if(TransferOperatorType::block_type::rows==1)
Dune::MatrixVector::addProduct(cm, (*im)[0][0] * (*jm)[0][0], *m);
else
- StaticMatrix::addTransformedMatrix(cm, *im, *m, *jm);
+ Dune::MatrixVector::addTransformedMatrix(cm, *im, *m, *jm);
}
}
}
@@ -767,7 +765,7 @@ public:
if(TransferOperatorType::block_type::rows==1)
Dune::MatrixVector::addProduct(cm, (*im)[0][0] * (*jm)[0][0], *m);
else
- StaticMatrix::addTransformedMatrix(cm, *im, *m, *jm);
+ Dune::MatrixVector::addTransformedMatrix(cm, *im, *m, *jm);
}
}
}