From 73fa481cf640b5790ab6b14df6322f59c1afa313 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 22 Jan 2010 08:46:23 +0000
Subject: [PATCH] hack: make dependency on ag-common disappear by copying
staticmatrixtools.hh from ag-common to here
[[Imported from SVN: r3174]]
---
dune.module | 2 +-
dune/solvers/common/Makefile.am | 2 +-
dune/solvers/common/staticmatrixtools.hh | 320 ++++++++++++++++++
.../genericmultigridtransfer.hh | 2 +-
4 files changed, 323 insertions(+), 3 deletions(-)
create mode 100644 dune/solvers/common/staticmatrixtools.hh
diff --git a/dune.module b/dune.module
index dd6e2a32..ae30173e 100644
--- a/dune.module
+++ b/dune.module
@@ -7,4 +7,4 @@ Module: dune-solvers
Version: 0.9
Maintainer: sander@mi.fu-berlin.de
#depending on
-Depends: dune-common dune-grid dune-istl dune-localfunctions ag-common
+Depends: dune-common dune-grid dune-istl dune-localfunctions
diff --git a/dune/solvers/common/Makefile.am b/dune/solvers/common/Makefile.am
index fb1b08c5..ccc8e9ca 100644
--- a/dune/solvers/common/Makefile.am
+++ b/dune/solvers/common/Makefile.am
@@ -1,6 +1,6 @@
SUBDIRS =
commondir = $(includedir)/dune/dune-solvers/common
-common_HEADERS = canignore.hh
+common_HEADERS = canignore.hh staticmatrixtools.hh
include $(top_srcdir)/am/global-rules
diff --git a/dune/solvers/common/staticmatrixtools.hh b/dune/solvers/common/staticmatrixtools.hh
new file mode 100644
index 00000000..af26041b
--- /dev/null
+++ b/dune/solvers/common/staticmatrixtools.hh
@@ -0,0 +1,320 @@
+#ifndef STATIC_MATRIX_TOOL_HH
+#define STATIC_MATRIX_TOOL_HH
+
+#include "dune/common/fmatrix.hh"
+#include "dune/istl/scaledidmatrix.hh"
+#include "dune/istl/diagonalmatrix.hh"
+#include "dune/istl/bcrsmatrix.hh"
+#include "dune/istl/matrixindexset.hh"
+
+class StaticMatrix
+{
+ public:
+
+ // type promotion (smallest matrix that can hold the sum of two matrices *******************
+ template<class MatrixA, class MatrixB>
+ struct Promote
+ {
+ typedef Dune::FieldMatrix<typename MatrixA::field_type, MatrixA::rows, MatrixA::cols> Type;
+ };
+
+ template<class Matrix>
+ struct Promote<Matrix, Matrix>
+ {
+ typedef Matrix Type;
+ };
+
+ template<typename FieldType, int n>
+ struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n> >
+ {
+ typedef Dune::DiagonalMatrix<FieldType,n> Type;
+ };
+
+ template<typename FieldType, int n>
+ struct Promote<Dune::ScaledIdentityMatrix<FieldType,n>, Dune::DiagonalMatrix<FieldType,n> >
+ {
+ typedef Dune::DiagonalMatrix<FieldType,n> Type;
+ };
+
+
+
+ // add scalar to matrix diagonal ***********************************************************
+ template<class Matrix>
+ static void addToDiagonal(Matrix& x, const typename Matrix::field_type a)
+ {
+ for(int i=0; i<Matrix::rows; ++i)
+ x[i][i] += a;
+ }
+
+ template <typename FieldType, int n>
+ static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a)
+ {
+ 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)
+ {
+ for(int i=0; i<B.N(); ++i)
+ {
+ 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);
+ }
+ }
+
+// 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 FieldType a, const MatrixB& y)
+ {
+ x.scalar() += a * y[0][0];
+ }
+
+ template <typename FieldType, int n, class MatrixB>
+ static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const FieldType a, const MatrixB& y)
+ {
+ for(int i=0; i<n; ++i)
+ x.diagonal()[i] += a * y[i][i];
+ }
+
+
+
+ // 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 (int i=0; i<A.N(); ++i)
+ {
+ typename MatrixA::row_type::ConstIterator jIt=A[i].begin();
+ typename MatrixA::row_type::ConstIterator jEnd=A[i].end();
+ for(; jIt!=jEnd; ++jIt)
+ {
+ 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)
+ A[i][jIt.index()] += T1[k][i] * B[k][lIt.index()] * T2[lIt.index()][jIt.index()];
+ }
+ }
+ }
+ }
+
+ template<class MatrixA, class MatrixB>
+ static void addTransformedMatrix(MatrixA& X, const double& A, const MatrixB& Y, const double& B)
+ {
+ axpy(X, A*B, Y);
+ }
+
+ 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)
+ {
+ 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];
+ }
+
+ 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)
+ {
+ 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];
+ }
+
+ 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)
+ {
+ for (int i=0; i<n; ++i)
+ for (int k=0; k<n; k++)
+ X.diagonal()[i] += A[k][i] * Y.diagonal()[k] * B[k][i];
+ }
+
+ 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::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::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)
+ {
+ 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];
+ }
+
+ template<class MatrixA, class MatrixB>
+ static void transformMatrix(MatrixA& A, const double& X, const MatrixB& B, const double& Y)
+ {
+ A = 0.0;
+ axpy(A, X*Y, B);
+ }
+
+ template<class MatrixA, class MatrixT>
+ static void transformMatrix(MatrixA& A, const MatrixT& X, const double& B, const MatrixT& Y)
+ {
+ A = X*B*Y;
+ }
+
+ template<class MatrixBlockA, class MatrixB, class MatrixT>
+ static void transformMatrix(Dune::BCRSMatrix<MatrixBlockA>& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
+ {
+ typedef typename Dune::BCRSMatrix<MatrixBlockA> MatrixA;
+
+ transformMatrixPattern(A, X, B, Y);
+ 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 MatrixT::row_type::ConstIterator XkiIt = X[k].begin();
+ typename MatrixT::row_type::ConstIterator XkiEnd = X[k].end();
+ for (; XkiIt!=XkiEnd; ++XkiIt)
+ {
+ int i = XkiIt.index();
+
+ typename MatrixT::row_type::ConstIterator YljIt = Y[l].begin();
+ typename MatrixT::row_type::ConstIterator YljEnd = Y[l].end();
+ for (; YljIt!=YljEnd; ++YljIt)
+ {
+ int j = YljIt.index();
+ MatrixBlockA Aij;
+ transformMatrix(Aij,*XkiIt, *BklIt, *YljIt);
+ A[i][j] += Aij;
+ }
+ }
+ }
+ }
+ }
+
+ template<class MatrixBlockA, class MatrixB, class MatrixT>
+ static void transformMatrixPattern(Dune::BCRSMatrix<MatrixBlockA>& A, const MatrixT& X, const MatrixB& B, const MatrixT& Y)
+ {
+ typedef typename Dune::BCRSMatrix<MatrixBlockA> MatrixA;
+
+ Dune::MatrixIndexSet indices(X.M(), Y.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 MatrixT::row_type::ConstIterator XkiIt = X[k].begin();
+ typename MatrixT::row_type::ConstIterator XkiEnd = X[k].end();
+ for (; XkiIt!=XkiEnd; ++XkiIt)
+ {
+ int i = XkiIt.index();
+
+ typename MatrixT::row_type::ConstIterator YljIt = Y[l].begin();
+ typename MatrixT::row_type::ConstIterator YljEnd = Y[l].end();
+ for (; YljIt!=YljEnd; ++YljIt)
+ {
+ int j = YljIt.index();
+ indices.add(i, j);
+ }
+ }
+ }
+ }
+ indices.exportIdx(A);
+ }
+
+
+
+ // compute i-th row of A*v for an edge vector e = e_i - e_j with i!=j **********************
+ template <class K>
+ static Dune::FieldMatrix<K,1,1>& toMatrix(K& x)
+ {
+ return *(reinterpret_cast< Dune::FieldMatrix<K,1,1>* > (&x));
+ }
+
+ template <class K, int n, int m>
+ static Dune::FieldMatrix<K,n,m>& toMatrix(Dune::FieldMatrix<K,n,m>& x)
+ {
+ return x;
+ }
+
+ template <class K, int n>
+ static Dune::ScaledIdentityMatrix<K,n>& toMatrix(Dune::ScaledIdentityMatrix<K,n>& x)
+ {
+ return x;
+ }
+
+ template <class K, int n>
+ static Dune::DiagonalMatrix<K,n>& toMatrix(Dune::DiagonalMatrix<K,n>& x)
+ {
+ return x;
+ }
+
+ template <class K>
+ static Dune::FieldMatrix<K,1,1>& toMatrix(Dune::ScaledIdentityMatrix<K,1>& x)
+ {
+ return *(reinterpret_cast< Dune::FieldMatrix<K,1,1>* > (&x));
+ }
+
+ template <class K>
+ static Dune::FieldMatrix<K,1,1>& toMatrix(Dune::DiagonalMatrix<K,1>& x)
+ {
+ return *(reinterpret_cast< Dune::FieldMatrix<K,1,1>* > (&x));
+ }
+
+
+};
+
+#endif
+
diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index d9ab4195..686989c8 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh
@@ -10,7 +10,7 @@
#include <dune/localfunctions/lagrange/pqkfactory.hh>
-#include "dune/ag-common/staticmatrixtools.hh"
+#include "dune/solvers/common/staticmatrixtools.hh"
/** \brief Restriction and prolongation operator for standard multigrid
--
GitLab