diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh
index 70bd60caa9eee6b690ebbe1baa25545cc364672f..c9dfabe0679de3f370020b670bb89dc310e1b06d 100644
--- a/dune/solvers/common/staticmatrixtools.hh
+++ b/dune/solvers/common/staticmatrixtools.hh
@@ -3,16 +3,16 @@
#ifndef STATIC_MATRIX_TOOL_HH
#define STATIC_MATRIX_TOOL_HH
-#include "dune/common/diagonalmatrix.hh"
#include "dune/common/fmatrix.hh"
+#include "dune/common/diagonalmatrix.hh"
#include "dune/istl/scaledidmatrix.hh"
#include "dune/istl/bcrsmatrix.hh"
#include "dune/istl/matrixindexset.hh"
-class StaticMatrix
-{
- public:
+#include "arithmetic.hh"
+namespace StaticMatrix
+{
// type promotion (smallest matrix that can hold the sum of two matrices *******************
template<class MatrixA, class MatrixB>
struct Promote
@@ -51,7 +51,6 @@ class StaticMatrix
};
-
// add scalar to matrix diagonal ***********************************************************
template<class Matrix>
static void addToDiagonal(Matrix& x, const typename Matrix::field_type a)
@@ -73,216 +72,234 @@ class StaticMatrix
x.scalar() += a;
}
-
-
- // add scalar times matrix to matrix *******************************************************
- template<class MatrixA, class MatrixB>
- static void axpy(MatrixA& A, const typename MatrixA::field_type a, const MatrixB& B)
+ // add transformed matrix A += T1^t*B*T2 ******************************************************
+ template<class MatrixA, class MatrixB, class TransformationMatrix, bool AisScalar, bool BisScalar, bool TisScalar>
+ struct TransformMatrixHelper
{
- for(size_t i=0; i<B.N(); ++i)
+ static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
- typename MatrixB::row_type::ConstIterator it=B[i].begin();
- typename MatrixB::row_type::ConstIterator end=B[i].end();
- for(; it!=end; ++it)
- A[i][it.index()] += a * (*it);
+ 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()];
+ }
+ }
+ }
+ }
}
- }
-// template<class MatrixA, class MatrixB>
-// static void axpy(MatrixA& x, const typename MatrixA::field_type a, const MatrixB& y)
-// {
-// for(int i=0; i<MatrixA::rows; ++i)
-// for(int j=0; j<MatrixA::cols; ++j)
-// x[i][j] += a * y[i][j];
-// }
-
- template <typename FieldType, int n, class MatrixB>
- static void axpy(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a, const MatrixB& y)
- {
- x.scalar() += a * y[0][0];
- }
+// 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 <typename FieldType, int n, class MatrixB>
- static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const typename Dune::DiagonalMatrix<FieldType,n>::field_type a, const MatrixB& y)
+ 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>
{
- for(int i=0; i<n; ++i)
- x.diagonal()[i] += a * y[i][i];
- }
+ typedef Dune::FieldMatrix<K1, m, m> MatrixA;
+ typedef Dune::FieldMatrix<K3, n, n> MatrixB;
+ typedef Dune::FieldMatrix<K2, n, m> TransformationMatrix;
- template <typename FieldType, int n, typename Scalar>
- static void axpy(Dune::FieldVector<FieldType,n>& x, const Scalar a, const Dune::FieldVector<FieldType,n>& y)
- {
- x.axpy(a, y);
- }
-
-
-
-
-
- // add transformed matrix X = A^T*Y*B ******************************************************
- template<class MatrixA, class MatrixB, class TransformationMatrix>
- static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
- {
- for (size_t i=0; i<A.N(); ++i)
+ static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
- typename MatrixA::row_type::ConstIterator jIt=A[i].begin();
- typename MatrixA::row_type::ConstIterator jEnd=A[i].end();
- for(; jIt!=jEnd; ++jIt)
+ Dune::FieldMatrix<K1, m, n> T1transposedB;
+ T1transposedB = 0;
+ for (size_t i=0; i<MatrixA::rows; ++i)
{
- for (size_t k=0; k<B.N(); ++k)
+ for (size_t k=0; k<MatrixB::rows; ++k)
{
if (T1[k][i]!=0)
- {
- typename MatrixB::row_type::ConstIterator lIt=B[k].begin();
- typename MatrixB::row_type::ConstIterator lEnd=B[k].end();
- for(; lIt!=lEnd; ++lIt)
- A[i][jIt.index()] += T1[k][i] * B[k][lIt.index()] * T2[lIt.index()][jIt.index()];
- }
+ 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];
}
}
}
- }
- template<class K1, class K2, class K3, int n, int m>
- static void addTransformedMatrix(
- typename Dune::FieldMatrix<K1, m, m>& A,
- const typename Dune::FieldMatrix<K2, n, m>& T1,
- const typename Dune::FieldMatrix<K3, n, n>& B,
- const typename Dune::FieldMatrix<K2, n, m>& T2)
- {
- typename Dune::FieldMatrix<K1, m, n> T1transposedB;
- T1transposedB = 0;
- for (size_t i=0; i<A.N(); ++i)
+// 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)
{
- for (size_t k=0; k<B.N(); ++k)
- {
- if (T1[k][i]!=0)
- for (size_t l=0; l<B.M(); ++l)
- T1transposedB[i][l] += T1[k][i] * B[k][l];
- }
+ Arithmetic::addProduct(A, T1*T2, B);
}
- for (size_t k=0; k<T2.N(); ++k)
+ };
+
+ 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 l=0; l<T2.M(); ++l)
+ for (size_t k=0; k<TransformationMatrix::rows; ++k)
{
- if (T2[k][l]!=0)
- for (size_t i=0; i<A.N(); ++i)
- A[i][l] += T1transposedB[i][k] * T2[k][l];
+ 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 MatrixA, class MatrixB>
- static void addTransformedMatrix(MatrixA& X, const double& A, const MatrixB& Y, const double& B)
+ template<class FieldType, int n, class ScalarB, class TransformationMatrix, bool AisScalar, bool TisScalar>
+ struct TransformMatrixHelper<Dune::ScaledIdentityMatrix<FieldType,n>, ScalarB, TransformationMatrix, AisScalar, true, TisScalar>
{
- axpy(X, A*B, Y);
- }
+ 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>
- static void addTransformedMatrix(MatrixA& X, const TransformationMatrix& A, const Dune::DiagonalMatrix<FieldType,n>& Y, const TransformationMatrix& B)
+ struct TransformMatrixHelper<MatrixA, Dune::DiagonalMatrix<FieldType,n>, TransformationMatrix, false, false, false>
{
- for (int i=0; i<MatrixA::rows; ++i)
- for (int j=0; j<MatrixA::cols; ++j)
- for (int k=0; k<n; k++)
- X[i][j] += A[k][i] * Y.diagonal()[k] * B[k][j];
- }
+ 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>
- static void addTransformedMatrix(MatrixA& X, const TransformationMatrix& A, const Dune::ScaledIdentityMatrix<FieldType,n>& Y, const TransformationMatrix& B)
+ struct TransformMatrixHelper<MatrixA, Dune::ScaledIdentityMatrix<FieldType,n>, TransformationMatrix, false, false, false>
{
- for (int i=0; i<MatrixA::rows; ++i)
- for (int j=0; j<MatrixA::cols; ++j)
- for (int k=0; k<n; k++)
- X[i][j] += A[k][i] * Y.scalar() * B[k][j];
- }
+ 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>
- static void addTransformedMatrix(Dune::DiagonalMatrix<FieldType,n>& X, const TransformationMatrix& A, const Dune::DiagonalMatrix<FieldType,n>& Y, const TransformationMatrix& B)
+ struct TransformMatrixHelper<Dune::DiagonalMatrix<FieldType,n>, Dune::DiagonalMatrix<FieldType,n>, TransformationMatrix, false, false, false>
{
- for (int i=0; i<n; ++i)
+ 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++)
- X.diagonal()[i] += A[k][i] * Y.diagonal()[k] * B[k][i];
- }
+ {
+ 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>
- static void addTransformedMatrix(Dune::ScaledIdentityMatrix<FieldType,n>& X, const TransformationMatrix& A, const Dune::ScaledIdentityMatrix<FieldType,n>& Y, const TransformationMatrix& B)
- {
- for (int k=0; k<n; k++)
- X.scalar() += A[k][0] * Y.scalar() * B[k][0];
- }
-
-
-
- // compute v^T*A*v for an edge vector e = e_i - e_j with i!=j ******************************
- template<class Matrix>
- static typename Matrix::field_type simplexEdgeDiagonal(const Matrix& A, int i, int j)
- {
- return A[i][i] - A[i][j] - A[j][i] + A[j][j];
- }
-
- template <typename FieldType, int n>
- static FieldType simplexEdgeDiagonal(const Dune::DiagonalMatrix<FieldType,n>& A, int i, int j)
- {
- return A.diagonal(i) + A.diagonal(j);
- }
-
- template <typename FieldType, int n>
- static FieldType simplexEdgeDiagonal(const Dune::ScaledIdentityMatrix<FieldType,n>& A, int i, int j)
- {
- return 2*A.scalar();
- }
-
-
-
- // compute i-th row of A*v for an edge vector e = e_i - e_j with i!=j **********************
- template<class Matrix>
- static typename Matrix::field_type simplexEdgeResidual(const Matrix& A, int i, int j)
- {
- return A[i][i] - A[i][j];
- }
-
- template <typename FieldType, int n>
- static FieldType simplexEdgeResidual(const Dune::DiagonalMatrix<FieldType,n>& A, int i, int j)
+ struct TransformMatrixHelper<Dune::ScaledIdentityMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n>, TransformationMatrix, false, false, false>
{
- return A.diagonal(i);
- }
+ 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>
- static FieldType simplexEdgeResidual(const Dune::ScaledIdentityMatrix<FieldType,n>& A, int i, int j)
- {
- return A.scalar();
- }
-
-
- // transform matrix A = X^T * B * Y ********************************************************
- template<class MatrixA, class MatrixB, class MatrixT>
- static void transformMatrix(MatrixA& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
+ struct TransformMatrixHelper<Dune::ScaledIdentityMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n>, false, false, false>
{
- A = 0.0;
- for (int k=0; k<B.N(); ++k)
- for (int l=0; l<B.M(); ++l)
- for (int i=0; i<X.M(); ++i)
- for (int j=0; j<Y.M(); ++j)
- A[i][j] = X[k][i] * B[k][l] * Y[l][j];
- }
+ 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>
- static void transformMatrix(MatrixA& A, const double& X, const MatrixB& B, const double& Y)
+ template<class MatrixA, class MatrixB, class TransformationMatrix>
+ void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
- A = 0.0;
- axpy(A, X*Y, B);
+ TransformMatrixHelper<MatrixA,MatrixB,TransformationMatrix,
+ Arithmetic::ScalarTraits<MatrixA>::isScalar, Arithmetic::ScalarTraits<MatrixB>::isScalar, Arithmetic::ScalarTraits<TransformationMatrix>::isScalar>
+ ::addTransformedMatrix(A,T1,B,T2);
}
- template<class MatrixA, class MatrixT>
- static void transformMatrix(MatrixA& A, const MatrixT& X, const double& B, const MatrixT& Y)
+ template<class MatrixA, class MatrixB, class TransformationMatrix>
+ void transformMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{
- A = X*B*Y;
+ A = 0;
+ TransformMatrixHelper<MatrixA,MatrixB,TransformationMatrix,
+ Arithmetic::ScalarTraits<MatrixA>::isScalar, Arithmetic::ScalarTraits<MatrixB>::isScalar, Arithmetic::ScalarTraits<TransformationMatrix>::isScalar>
+ ::addTransformedMatrix(A,T1,B,T2);
}
- template<class MatrixBlockA, class MatrixB, class MatrixT>
- static void transformMatrix(Dune::BCRSMatrix<MatrixBlockA>& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
+ 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, X, B, Y);
+ transformMatrixPattern(A, T1, B, T2);
A = 0.0;
for (int k=0; k<B.N(); ++k)
{
@@ -292,19 +309,19 @@ class StaticMatrix
{
int l = BklIt.index();
- typename MatrixT::row_type::ConstIterator XkiIt = X[k].begin();
- typename MatrixT::row_type::ConstIterator XkiEnd = X[k].end();
- for (; XkiIt!=XkiEnd; ++XkiIt)
+ typename TransformationMatrix::row_type::ConstIterator T1kiIt = T1[k].begin();
+ typename TransformationMatrix::row_type::ConstIterator T1kiEnd = T1[k].end();
+ for (; T1kiIt!=T1kiEnd; ++T1kiIt)
{
- int i = XkiIt.index();
+ int i = T1kiIt.index();
- typename MatrixT::row_type::ConstIterator YljIt = Y[l].begin();
- typename MatrixT::row_type::ConstIterator YljEnd = Y[l].end();
- for (; YljIt!=YljEnd; ++YljIt)
+ typename TransformationMatrix::row_type::ConstIterator T2ljIt = T2[l].begin();
+ typename TransformationMatrix::row_type::ConstIterator T2ljEnd = T2[l].end();
+ for (; T2ljIt!=T2ljEnd; ++T2ljIt)
{
- int j = YljIt.index();
+ int j = T2ljIt.index();
MatrixBlockA Aij;
- transformMatrix(Aij,*XkiIt, *BklIt, *YljIt);
+ transformMatrix(Aij,*T1kiIt, *BklIt, *T2ljIt);
A[i][j] += Aij;
}
}
@@ -312,10 +329,10 @@ class StaticMatrix
}
}
- template<class MatrixBlockA, class MatrixB, class MatrixT>
- static void transformMatrixPattern(Dune::BCRSMatrix<MatrixBlockA>& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
+ 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(X.M(), Y.M());
+ Dune::MatrixIndexSet indices(T1.M(), T2.M());
for (int k=0; k<B.N(); ++k)
{
typename MatrixB::row_type::ConstIterator BklIt = B[k].begin();
@@ -324,17 +341,17 @@ class StaticMatrix
{
int l = BklIt.index();
- typename MatrixT::row_type::ConstIterator XkiIt = X[k].begin();
- typename MatrixT::row_type::ConstIterator XkiEnd = X[k].end();
- for (; XkiIt!=XkiEnd; ++XkiIt)
+ typename TransformationMatrix::row_type::ConstIterator T1kiIt = T1[k].begin();
+ typename TransformationMatrix::row_type::ConstIterator T1kiEnd = T1[k].end();
+ for (; T1kiIt!=T1kiEnd; ++T1kiIt)
{
- int i = XkiIt.index();
+ int i = T1kiIt.index();
- typename MatrixT::row_type::ConstIterator YljIt = Y[l].begin();
- typename MatrixT::row_type::ConstIterator YljEnd = Y[l].end();
- for (; YljIt!=YljEnd; ++YljIt)
+ typename TransformationMatrix::row_type::ConstIterator T2ljIt = T2[l].begin();
+ typename TransformationMatrix::row_type::ConstIterator T2ljEnd = T2[l].end();
+ for (; T2ljIt!=T2ljEnd; ++T2ljIt)
{
- int j = YljIt.index();
+ int j = T2ljIt.index();
indices.add(i, j);
}
}
@@ -344,6 +361,48 @@ class StaticMatrix
}
+ // compute v^T*A*v for an edge vector e = e_i - e_j with i!=j ******************************
+ template<class Matrix>
+ static typename Matrix::field_type simplexEdgeDiagonal(const Matrix& A, int i, int j)
+ {
+ return A[i][i] - A[i][j] - A[j][i] + A[j][j];
+ }
+
+ template <typename FieldType, int n>
+ static FieldType simplexEdgeDiagonal(const Dune::DiagonalMatrix<FieldType,n>& A, int i, int j)
+ {
+ return A.diagonal(i) + A.diagonal(j);
+ }
+
+ template <typename FieldType, int n>
+ static FieldType simplexEdgeDiagonal(const Dune::ScaledIdentityMatrix<FieldType,n>& A, int i, int j)
+ {
+ return 2*A.scalar();
+ }
+
+
+
+ // compute i-th row of A*v for an edge vector e = e_i - e_j with i!=j **********************
+ template<class Matrix>
+ static typename Matrix::field_type simplexEdgeResidual(const Matrix& A, int i, int j)
+ {
+ return A[i][i] - A[i][j];
+ }
+
+ template <typename FieldType, int n>
+ static FieldType simplexEdgeResidual(const Dune::DiagonalMatrix<FieldType,n>& A, int i, int j)
+ {
+ return A.diagonal(i);
+ }
+
+ template <typename FieldType, int n>
+ static FieldType simplexEdgeResidual(const Dune::ScaledIdentityMatrix<FieldType,n>& A, int i, int j)
+ {
+ return A.scalar();
+ }
+
+
+
// compute i-th row of A*v for an edge vector e = e_i - e_j with i!=j **********************
template <class K>
@@ -458,7 +517,10 @@ class StaticMatrix
}
}
-};
+
+
+
+}
#endif
diff --git a/dune/solvers/iterationsteps/truncatedblockgsstep.hh b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
index bd9b1433a99896959c62f4411f01a10fe3a0304a..c0c2083810079cf36e42e1aa97508c4be03e28dc 100644
--- a/dune/solvers/iterationsteps/truncatedblockgsstep.hh
+++ b/dune/solvers/iterationsteps/truncatedblockgsstep.hh
@@ -186,9 +186,7 @@ public:
{
size_t col = it.index();
if (col == row)
- {
- StaticMatrix::axpy(Aii, 1.0, *it);
- }
+ Arithmetic::addProduct(Aii, 1.0, *it);
else
it->mmv(x[col],r);
}