staticmatrixtools.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef STATIC_MATRIX_TOOL_HH
#define STATIC_MATRIX_TOOL_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"

#include <dune/matrix-vector/addtodiagonal.hh>
#include <dune/matrix-vector/axpy.hh>
#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"

namespace StaticMatrix
{
        // type promotion (smallest matrix that can hold the sum of two matrices *******************
        template<class MatrixA, class MatrixB>
        struct Promote
        {
          typedef typename Dune::MatrixVector::Promote<MatrixA, MatrixB>::Type Type;
        };

        // add scalar to matrix diagonal ***********************************************************
        template<class Matrix>
        static void addToDiagonal(Matrix& x, const typename Matrix::field_type a)
        {
          Dune::MatrixVector::addToDiagonal(x, a);
        }

        // add transformed matrix A += T1^t*B*T2 ******************************************************
        template<class MatrixA, class MatrixB, class TransformationMatrix>
        void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& 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)
        {
            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::MatrixVector::transformMatrixPattern(A,T1,B,T2);
        }


        // 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();
        }




        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));
        }

        /** \brief Compute an LDL^T decomposition
         *
         * The methods computes a decomposition A=LDL^T of a given dense
         * symmetric matrix A such that L is lower triangular with all
         * diagonal elements equal to 1 and D is diagonal. If A is positive
         * definite then A=(LD^0.5)(LD^0.5)^T is the Cholesky decomposition.
         * However, the LDL^T decomposition does also work for indefinite
         * symmetric matrices and is more stable than the Cholesky decomposition
         * since no square roots are required.
         *
         * The method does only change the nontrivial entries of the given matrix
         * L and D, i.e., it does not set the trivial 0 and 1 entries.
         * Thus one can store both in a single matrix M and use
         * M as argument for L and D.
         *
         * The method can furthermore work in-place, i.e., it is safe to
         * use A as argument for L and D. In this case the entries of A
         * below and on the diagonal are changed to those to those of
         * L and D, respectively.
         *
         * \param A Matrix to be decomposed. Only the lower triangle is used.
         * \param L Matrix to store the lower triangle. Only entries below the diagonal are written.
         * \param D Matrix to store the diagonal. Only diagonal entries are written.
         */
        template<class SymmetricMatrix, class LowerTriangularMatrix, class DiagonalMatrix>
        static void ldlt(const SymmetricMatrix& A, LowerTriangularMatrix& L, DiagonalMatrix& D)
        {
          Dune::MatrixVector::ldlt(A,L,D);
        }

        /** \brief Solve linear system using a LDL^T decomposition.
         *
         * The methods solves a linear system Mx=b where A is given
         * by a decomposition A=LDL^T. The method does only use
         * the values of L and D below and on the diagonal, respectively.
         * The 1 entries on the diagonal of L are not required.
         * If L and D are stored in a single matrix it is safe
         * the use this matrix as argument for both.
         *
         * Note that the solution vector must already have the correct size.
         *
         * \param L Matrix containing the lower triangle of the decomposition
         * \param D Matrix containing the diagonal of the decomposition
         * \param b Right hand side on the linear system
         * \param x Vector to store the solution of the linear system
         */
        template<class LowerTriangularMatrix, class DiagonalMatrix, class RhsVector, class SolVector>
        static void ldltSolve(const LowerTriangularMatrix& L, const DiagonalMatrix& D, const RhsVector& b, SolVector& x)
        {
          Dune::MatrixVector::ldltSolve(L,D,b,x);
        }
}

#endif