hack: make dependency on ag-common disappear by copying staticmatrixtools.hh from ag-common to here

[[Imported from SVN: r3174]]

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

dune/solvers/common/Makefile.am

 SUBDIRS =
 
 commondir = $(includedir)/dune/dune-solvers/common
-common_HEADERS = canignore.hh
+common_HEADERS = canignore.hh staticmatrixtools.hh
 
 include $(top_srcdir)/am/global-rules

dune/solvers/common/staticmatrixtools.hh

+#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
+

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