Forked from
agnumpde / dune-solvers
88 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
genericvectortools.hh 10.17 KiB
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef GENERIC_VECTOR_TOOL_HH
#define GENERIC_VECTOR_TOOL_HH
/** \file
\brief Various tools for working with istl vectors of arbitrary nesting depth
*/
#include<iostream>
#include<bitset>
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/indices.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/istl/scaledidmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/matrix-vector/genericvectortools.hh>
#include <dune/matrix-vector/resize.hh>
#include <dune/solvers/operators/sumoperator.hh>
/** \brief Various tools for working with istl vectors of arbitrary nesting depth
*/
struct GenericVector
{
//! Write vector to given stream
template <class VectorType>
DUNE_DEPRECATED_MSG("Please use Dune::MatrixVector::Generic::writeBinary instead.")
static void writeBinary(std::ostream& s, const VectorType& v)
{
Dune::MatrixVector::Generic::writeBinary(s, v);
}
//! Read vector from a given stream
template <class VectorType>
DUNE_DEPRECATED_MSG("Please use Dune::MatrixVector::Generic::readBinary instead.")
static void readBinary(std::istream& s, VectorType& v)
{
Dune::MatrixVector::Generic::readBinary(s, v);
}
//! Resize vector recursivly to size of given vector/matrix
template <class VectorTypeA, class VectorTypeB>
DUNE_DEPRECATED_MSG("Please use Dune::MatrixVector::resize instead.")
static void resize(VectorTypeA& a, const VectorTypeB& b)
{
Dune::MatrixVector::resize(a, b);
}
//! Set vector to zero at indices that are true in bitvector recursivly
template <class VectorType, class BitVectorType>
DUNE_DEPRECATED_MSG("Please use Dune::MatrixVector::Generic::truncate instead.")
static void truncate(VectorType& v, const BitVectorType& tr)
{
Dune::MatrixVector::Generic::truncate(v, tr);
}
//! weave two vector blocks into each other e.g. [[u1...un][w1...wn]] --> [[u1 w1]...[un wn]]
template <class LVectorType, class RVectorType>
static void interlace(const LVectorType& lvec DUNE_UNUSED, RVectorType& rvec DUNE_UNUSED)
{ DUNE_THROW(Dune::NotImplemented,"GenericVector::interlace not implemented for given VectorTypes");
}
/** \brief weave two vector blocks into each other e.g. [[u1...un][w1...wn]] --> [[u1 w1]...[un wn]]
*
* interlace: \f$\left((\mathbb R^k)^m\right)^n\longrightarrow(\mathbb R^k\cdot n)^m\f$, i.e. for example
*
* \f[
* \begin{pmatrix}
* \begin{pmatrix}
* \begin{pmatrix} a_{11}\\ a_{12} \end{pmatrix}\\
* \begin{pmatrix} a_{21}\\ a_{22} \end{pmatrix}\\
* \begin{pmatrix} a_{31}\\ a_{32} \end{pmatrix}
* \end{pmatrix}\\
* \begin{pmatrix}
* \begin{pmatrix} b_{11}\\ b_{12} \end{pmatrix}\\
* \begin{pmatrix} b_{21}\\ b_{22} \end{pmatrix}\\
* \begin{pmatrix} b_{31}\\ b_{32} \end{pmatrix}
* \end{pmatrix}\\
* \begin{pmatrix}
* \begin{pmatrix} c_{11}\\ c_{12} \end{pmatrix}\\
* \begin{pmatrix} c_{21}\\ c_{22} \end{pmatrix}\\
* \begin{pmatrix} c_{31}\\ c_{32} \end{pmatrix}
* \end{pmatrix}\\
* \begin{pmatrix}
* \begin{pmatrix} d_{11}\\ d_{12} \end{pmatrix}\\
* \begin{pmatrix} d_{21}\\ d_{22} \end{pmatrix}\\
* \begin{pmatrix} d_{31}\\ d_{32} \end{pmatrix}
* \end{pmatrix}
* \end{pmatrix}
*
* \xrightarrow{interlace}
*
* \begin{pmatrix}
* \begin{pmatrix}
* a_{11}\\ a_{12}\\ b_{11}\\ b_{12}\\ c_{11}\\ c_{12}\\ d_{11}\\ d_{12}
* \end{pmatrix}\\
* \begin{pmatrix}
* a_{21}\\ a_{22}\\ b_{21}\\ b_{22}\\ c_{21}\\ c_{22}\\ d_{21}\\ d_{22}
* \end{pmatrix}\\
* \begin{pmatrix}
* a_{31}\\ a_{32}\\ b_{31}\\ b_{32}\\ c_{31}\\ c_{32}\\ d_{31}\\ d_{32}
* \end{pmatrix}
* \end{pmatrix}
* \f]
*/
template <class LVectorBlock, class RVectorBlock>
static void interlace(const Dune::BlockVector<Dune::BlockVector<LVectorBlock> >& lvec, Dune::BlockVector<RVectorBlock>& rvec)
{
static_assert(RVectorBlock::dimension % LVectorBlock::dimension == 0, "Block sizes don't match.");
rvec.resize(lvec[0].size());
rvec = 0.0;
size_t offset = LVectorBlock::dimension;
for (size_t i=0; i<lvec[0].size(); ++i)
for (size_t k=0; k<lvec.size(); ++k)
for (size_t j=0; j<offset; ++j)
rvec[i][j+k*offset] = lvec[k][i][j];
}
//! unweave two vectors previously interlaced e.g. [[u1 w1]...[un wn]] --> [[u1...un][w1...wn]]
template <class LVectorType, class RVectorType>
static void deinterlace(const LVectorType& lvec DUNE_UNUSED, RVectorType& rvec DUNE_UNUSED)
{
DUNE_THROW(Dune::NotImplemented,"GenericVector::deinterlace not implemented for given VectorTypes");
}
template <class LVectorBlock, class RVectorBlock> static void deinterlace(const Dune::BlockVector<LVectorBlock>& lvec, Dune::BlockVector<Dune::BlockVector<RVectorBlock> >& rvec)
{
static_assert(LVectorBlock::dimension % RVectorBlock::dimension == 0, "Block sizes don't match.");
rvec.resize(LVectorBlock::dimension / RVectorBlock::dimension);
for (size_t k=0; k<rvec.size(); ++k)
rvec[k].resize(lvec.size());
rvec = 0.0;
size_t offset = RVectorBlock::dimension;
for (size_t i=0; i<lvec.size(); ++i)
for (size_t k=0; k<rvec.size(); ++k)
for (size_t j=0; j<offset; ++j)
rvec[k][i][j] = lvec[i][j+k*offset];
}
template <class VectorType>
static std::size_t isnan(const VectorType& x)
{
std::size_t r = 0;
typename VectorType::ConstIterator x_it=x.begin();
typename VectorType::ConstIterator x_end=x.end();
for (;x_it!=x_end; ++x_it)
r += isnan(*x_it);
return r;
}
template <class field_type, int n>
static std::size_t isnan(const Dune::FieldVector<field_type,n>& x)
{
std::size_t r=0;
for (std::size_t i=0; i<n; ++i)
r += (size_t)(std::fpclassify(x[i])==FP_NAN);
return r;
}
};
struct GenericMatrix
{
//! Set matrix to zero at indices that are true in bitvector recursivly
template <class MatrixType, class BitVectorTypeR, class BitVectorTypeC>
static void truncate(MatrixType& mat, const BitVectorTypeR& trows, const BitVectorTypeC& tcols, bool setTruncatedDiagonalOne)
{
assert(trows.size()==mat.N());
assert(tcols.size()==mat.M());
typename MatrixType::RowIterator row_it = mat.begin();
typename MatrixType::RowIterator row_end = mat.end();
for(; row_it!=row_end; ++row_it)
{
size_t row = row_it.index();
typename MatrixType::row_type::iterator col_it = row_it->begin();
typename MatrixType::row_type::iterator col_end = row_it->end();
for (; col_it!=col_end; ++col_it)
{
size_t col = col_it.index();
GenericMatrix::truncate(*col_it, trows[row], tcols[col], setTruncatedDiagonalOne and (col==row));
}
}
}
template <class field_type, int n, int m, class BitVectorTypeR, class BitVectorTypeC>
static void truncate(Dune::FieldMatrix<field_type,n,m>& mat, const BitVectorTypeR& trows, const BitVectorTypeC& tcols, bool setTruncatedDiagonalOne)
{// static_assert(BitVectorTypeR.size()==n and BitVectorTypeC.size()==m,"BitVector length doesn't match matrix size.");
assert(trows.size()==n and tcols.size()==m);
for(size_t row=0; row<n; ++row)
{
if (trows[row])
{
mat[row] = 0.0;
if (setTruncatedDiagonalOne)
mat[row][row] = 1.0;
continue;
}
for (size_t col=0; col<m; ++col)
{
if (tcols[col])
mat[row][col] = (setTruncatedDiagonalOne and col==row)?1.0:0.0;
}
}
}
template <class field_type, int n, class BitVectorType>
static void truncate(Dune::DiagonalMatrix<field_type,n>& mat, const BitVectorType& trows, const BitVectorType& tcols, bool setTruncatedDiagonalOne)
{
for(size_t row=0; row<n; ++row)
{
if (trows[row])
{
mat.diagonal(row) = setTruncatedDiagonalOne?1.0:0.0;
continue;
}
if (tcols[row])
mat.diagonal(row) = setTruncatedDiagonalOne?1.0:0.0;
}
}
template <class MatrixType>
static std::size_t isnan(const MatrixType& A)
{
std::size_t r=0;
typename MatrixType::ConstRowIterator row_it=A.begin();
typename MatrixType::ConstRowIterator row_end=A.end();
for (;row_it!=row_end; ++row_it)
{
typename MatrixType::ConstColIterator col_it=row_it->begin();
typename MatrixType::ConstColIterator col_end=row_it->end();
for (;col_it!=col_end; ++col_it)
r += isnan(*col_it);
}
return r;
}
template <class field_type, int n, int m>
static std::size_t isnan(const Dune::FieldMatrix<field_type,n,m>& x)
{
std::size_t r=0;
for (std::size_t i=0; i<n; ++i)
for (std::size_t j=0; j<m; ++j)
r += (size_t)(std::fpclassify(x[i][j])==FP_NAN);
if (r!=0)
std::cout << x << std::endl;
return r;
}
template <class field_type, int n>
static std::size_t isnan(const Dune::DiagonalMatrix<field_type,n>& x)
{
std::size_t r=0;
for (std::size_t i=0; i<n; ++i)
r += (size_t)(std::fpclassify(x.diagonal(i))==FP_NAN);
if (r!=0)
std::cout << x << std::endl;
return r;
}
template <class field_type, int n>
static std::size_t isnan(const Dune::ScaledIdentityMatrix<field_type,n>& x)
{
return n*(size_t)(std::fpclassify(x.scalar())==FP_NAN);
}
};
#endif