-
Elias Pipping authored
This reverts commit 0e1afee4. Adding these functions to dune-matrix-vector directly instead
Elias Pipping authoredThis reverts commit 0e1afee4. Adding these functions to dune-matrix-vector directly instead
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
staticmatrixtools.hh 18.48 KiB
// -*- 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/diagonalmatrix.hh"
#include "dune/common/fmatrix.hh"
#include "dune/istl/scaledidmatrix.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::FieldMatrix<FieldType,n,n>, Dune::DiagonalMatrix<FieldType,n> >
{
typedef Dune::FieldMatrix<FieldType,n,n> Type;
};
template<typename FieldType, int n>
struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::FieldMatrix<FieldType,n,n> >
{
typedef Dune::FieldMatrix<FieldType,n,n> 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)
// the commented line should NOT be used as it may lead to ambiguities in which case the general method above will be used.
// an example is the line taken from massassembler.hh (line 83 at the time):
// StaticMatrix::addToDiagonal(localMatrix[i][i], values[i] * zi);
// where values[i] is of type FieldVector<FieldType,1> and zi is a double. This call then does not exactly fit this specialization (without the implicit cast of FV<FT,1>)
// and hence some wild template voodoo decides which of the templates is to be taken - in this case with a gcc 4.4.5 that was the one above leading to a wrong factor n
// in the diagonal value static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type 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(size_t 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 typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type 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 typename Dune::DiagonalMatrix<FieldType,n>::field_type a, const MatrixB& y)
{
for(int i=0; i<n; ++i)
x.diagonal()[i] += a * y[i][i];
}
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)
{
typename MatrixA::row_type::ConstIterator jIt=A[i].begin();
typename MatrixA::row_type::ConstIterator jEnd=A[i].end();
for(; jIt!=jEnd; ++jIt)
{
for (size_t k=0; k<B.N(); ++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()];
} }
}
}
}
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)
{
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];
}
}
for (size_t k=0; k<T2.N(); ++k)
{
for (size_t l=0; l<T2.M(); ++l)
{
if (T2[k][l]!=0)
for (size_t i=0; i<A.N(); ++i)
A[i][l] += T1transposedB[i][k] * T2[k][l];
}
}
}
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::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();
}
// 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)
{
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)
{
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));
}
/** \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)
{
for(unsigned int i = 0; i < A.N(); ++i)
{
D[i][i] = A[i][i];
for(unsigned int j = 0; j < i; ++j)
{
L[i][j] = A[i][j];
for(unsigned int k = 0; k < j; ++k)
L[i][j] -= L[i][k] * L[j][k] * D[k][k];
L[i][j] /= D[j][j]; }
for(unsigned int k = 0; k < i; ++k)
D[i][i] -= L[i][k]*L[i][k] * D[k][k];
}
}
/** \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)
{
for(unsigned int i = 0; i < x.size(); ++i)
{
x[i] = b[i];
for(unsigned int j = 0; j < i; ++j)
x[i] -= L[i][j] * x[j];
}
for(unsigned int i = 0; i < x.size(); ++i)
x[i] /= D[i][i];
for(int i = x.size()-1; i >=0; --i)
{
for(unsigned int j = i+1; j < x.size(); ++j)
x[i] -= L[j][i] * x[j];
}
}
};
#endif