Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • jakub.both/dune-istl
  • eduardo.bueno/dune-istl
  • tkoch/dune-istl
  • pipping/dune-istl
  • andreas.brostrom/dune-istl
  • stephan.hilb/dune-istl
  • max.kahnt/dune-istl
  • patrick.jaap/dune-istl
  • tobias.meyer.andersen/dune-istl
  • andreas.thune/dune-istl
  • dominic/dune-istl
  • ansgar/dune-istl
  • exadune/dune-istl
  • lars.lubkoll/dune-istl
  • govind.sahai/dune-istl
  • maikel.nadolski/dune-istl
  • markus.blatt/dune-istl
  • core/dune-istl
  • lisa_julia.nebel/dune-istl
  • michael.sghaier/dune-istl
  • liesel.schumacher/dune-istl
  • Xinyun.Li/dune-istl
  • lukas.renelt/dune-istl
  • simon.praetorius/dune-istl
  • rene.milk/dune-istl
  • lasse.hinrichsen/dune-istl
  • nils.dreier/dune-istl
  • claus-justus.heine/dune-istl
  • felix.mueller/dune-istl
  • kilian.weishaupt/dune-istl
  • lorenzo.cerrone/dune-istl
  • jakob.torben/dune-istl
  • liam.keegan/dune-istl
  • alexander.mueller/dune-istl
34 results
Show changes
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1687 additions and 866 deletions
doc/doxygen/mainpage.txt
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
/** \mainpage dune-istl Automatic Documentation
\section intro Introduction
......
doc/doxygen/modules.txt
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
/**
@defgroup ISTL Iterative Solvers Template Library (ISTL)
......
doc/istl.bib
View file @ bc1f69fc
@Comment SPDX-FileCopyrightInfo: Copyright DUNE Project contributors, see file LICENSE.md in module root
@Comment SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
@String{ap = "Ann. der Physik"}
@String{as = "Acta Stereol."}
@String{awr = "Adv. Water Res."}
......@@ -156,4 +159,4 @@
OPTnote = {},
OPTannote = {},
publisher = {Springer-Verlag Wien}
}
\ No newline at end of file
}
doc/istl.tex
View file @ bc1f69fc
% SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
% SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
\documentclass[11pt]{article}
\usepackage{multicol}
\usepackage{ifthen}
%\usepackage{multitoc}
%\usepackage{german}
%\usepackage{bibgerm}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{color}
\usepackage{hyperref}
\usepackage[dvips]{epsfig}
\usepackage{subfigure}
\usepackage[dvips]{graphicx}
\usepackage[a4paper,body={148mm,240mm,nohead}]{geometry}
\usepackage[ansinew]{inputenc}
\usepackage{listings}
\lstset{language=C++, basicstyle=\ttfamily,
stringstyle=\ttfamily, commentstyle=\it, extendedchars=true}
\newif\ifpdf
\ifnum\ifx\pdfoutput\undefined0\else\pdfoutput\fi<1
\pdffalse % we are not running PDFLaTeX
\else
\pdftrue % we are running PDFLaTeX
\fi
\ifpdf
\usepackage[pdftex]{graphicx}
\else
\usepackage{graphicx}
\fi
\ifpdf
\DeclareGraphicsExtensions{.pdf, .jpg, .tif}
\else
\DeclareGraphicsExtensions{.eps, .jpg}
\fi
\newcommand{\C}{\mathbb{C}}
\newcommand{\R}{\mathbb{R}}
......@@ -668,7 +649,7 @@ contains the computed update to
the current guess, i.~e. $\vec v = M^{-1} \vec
d$ where $M$ is the approximate inverse of the operator $A$
characterizing the preconditioner. The method \lstinline!void post(X& x)!
should be called after all computations to give the precondtioner the
should be called after all computations to give the preconditioner the
chance to clean allocated resources.
See Table \ref{tab:precond} for a list of available
......@@ -831,7 +812,7 @@ now lack of computational efficiency due to the generic implementation.
% bibtex bibliography
\bibliographystyle{plain}
\bibliography{istl.bib}
\bibliography{istl}
% some links
% http://www.netlib.org/blas/blast-forum/
......
dune-istl.pc.in
View file @ bc1f69fc
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
prefix=@prefix@
exec_prefix=@exec_prefix@
libdir=@libdir@
......@@ -8,8 +11,8 @@ DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: Dune (Distributed and Unified Numerics Environment) istl module
URL: http://dune-project.org/
Description: @DESCRIPTION@
URL: @URL@
Requires: ${DEPENDENCIES}
Libs:
Cflags: -I${includedir}
dune.module
View file @ bc1f69fc
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
Module: dune-istl
Version: 2.7-git
Version: 2.11-git
Author: The Dune Core developers
Maintainer: dune-devel@lists.dune-project.org
Depends: dune-common (>= 2.6)
Description: Iterative solver template library which provides generic sparse matrix/vector classes and a variety of solvers based on these classes
URL: https://gitlab.dune-project.org/core/dune-istl
Python-Requires:
Depends: dune-common (>= 2.11)
Whitespace-Hook: Yes
dune/CMakeLists.txt
View file @ bc1f69fc
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
add_subdirectory(istl)
# if Python bindings are enabled, include necessary sub directories.
if( DUNE_ENABLE_PYTHONBINDINGS )
add_subdirectory("python")
endif()
dune/istl/CMakeLists.txt
View file @ bc1f69fc
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
add_subdirectory("common")
add_subdirectory("eigenvalue")
add_subdirectory("paamg")
add_subdirectory(test)
......@@ -6,12 +10,16 @@ add_subdirectory(test)
install(FILES
allocator.hh
basearray.hh
bccsmatrix.hh
bccsmatrixinitializer.hh
bcrsmatrix.hh
bdmatrix.hh
blocklevel.hh
btdmatrix.hh
bvector.hh
cholmod.hh
colcompmatrix.hh
dilu.hh
foreach.hh
gsetc.hh
ildl.hh
ilu.hh
......@@ -31,7 +39,6 @@ install(FILES
operators.hh
overlappingschwarz.hh
owneroverlapcopy.hh
pardiso.hh
preconditioner.hh
preconditioners.hh
repartition.hh
......@@ -40,6 +47,8 @@ install(FILES
schwarz.hh
solvercategory.hh
solver.hh
solverfactory.hh
solverregistry.hh
solvers.hh
solvertype.hh
spqr.hh
......
dune/istl/allocator.hh
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_ISTL_ALLOCATOR_HH
#define DUNE_ISTL_ALLOCATOR_HH
#include <dune/common/typetraits.hh>
#include <memory>
#include <type_traits>
#include <dune/common/typetraits.hh>
namespace Dune {
......@@ -18,7 +22,7 @@ namespace Dune {
};
template<typename T>
struct DefaultAllocatorTraits<T, void_t<typename T::allocator_type> >
struct DefaultAllocatorTraits<T, std::void_t<typename T::allocator_type> >
{
using type = typename T::allocator_type;
};
......@@ -30,7 +34,7 @@ namespace Dune {
using AllocatorType = typename AllocatorTraits<T>::type;
template<typename T, typename X>
using ReboundAllocatorType = typename AllocatorTraits<T>::type::template rebind<X>::other;
using ReboundAllocatorType = typename std::allocator_traits<typename AllocatorTraits<T>::type>::template rebind_alloc<X>;
} // end namespace Dune
......
dune/istl/basearray.hh
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_BASEARRAY_HH
#define DUNE_ISTL_BASEARRAY_HH
#include "assert.h"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <memory>
#include <algorithm>
#include <type_traits>
#include "istlexception.hh"
#include <dune/common/iteratorfacades.hh>
......@@ -30,7 +33,7 @@ namespace Imp {
- random access
This container has *NO* memory management at all,
copy constuctor, assignment and destructor are all default.
copy constructor, assignment and destructor are all default.
The constructor is made protected to emphasize that objects
are only usable in derived classes.
......@@ -39,13 +42,9 @@ namespace Imp {
Setting the compile time switch DUNE_ISTL_WITH_CHECKING
enables error checking.
\todo There shouldn't be an allocator argument here, because the array is 'unmanaged'.
And indeed, of the allocator, only its size_type is used. Hence, the signature
of this class should be changed to <class B, int stype>
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
template<class B, class ST=std::size_t >
class base_array_unmanaged
{
public:
......@@ -55,11 +54,8 @@ namespace Imp {
//! export the type representing the components
typedef B member_type;
//! export the allocator type
typedef A allocator_type;
//! the type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! the type used for references
using reference = B&;
......@@ -102,17 +98,27 @@ namespace Imp {
friend class RealIterator<ValueType>;
//! constructor
RealIterator ()
: p(0), i(0)
{}
RealIterator () = default;
RealIterator (const B* _p, B* _i) : p(_p), i(_i)
{ }
RealIterator (const B* _p, B* _i)
: p(_p), i(_i)
{}
RealIterator(const RealIterator<ValueType>& it)
: p(it.p), i(it.i)
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator (const RealIterator<T_>& other)
: p(other.p), i(other.i)
{}
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator& operator= (const RealIterator<T_>& other)
{
p = other.p;
i = other.i;
return *this;
}
//! return index
size_type index () const
{
......@@ -168,8 +174,8 @@ namespace Imp {
i+=d;
}
const B* p;
B* i;
const B* p = nullptr;
B* i = nullptr;
};
//! iterator type for sequential access
......@@ -289,7 +295,7 @@ namespace Imp {
- find returning iterator
This container has *NO* memory management at all,
copy constuctor, assignment and destructor are all default.
copy constructor, assignment and destructor are all default.
The constructor is made protected to emphasize that objects
are only usably in derived classes.
......@@ -300,7 +306,7 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
template<class B, class ST=std::size_t >
class compressed_base_array_unmanaged
{
public:
......@@ -310,11 +316,8 @@ namespace Imp {
//! export the type representing the components
typedef B member_type;
//! export the allocator type
typedef A allocator_type;
//! The type used for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! the type used for references
using reference = B&;
......@@ -357,22 +360,29 @@ namespace Imp {
friend class RealIterator<ValueType>;
//! constructor
RealIterator ()
: p(0), j(0), i(0)
{}
RealIterator () = default;
//! constructor
RealIterator (B* _p, size_type* _j, size_type _i)
: p(_p), j(_j), i(_i)
{ }
{}
/**
* @brief Copy constructor from mutable iterator
*/
RealIterator(const RealIterator<ValueType>& it)
: p(it.p), j(it.j), i(it.i)
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator (const RealIterator<T_>& other)
: p(other.p), j(other.j), i(other.i)
{}
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator& operator= (const RealIterator<T_>& other)
{
p = other.p;
j = other.j;
i = other.i;
return *this;
}
//! equality
bool equals (const RealIterator<ValueType>& it) const
......@@ -432,9 +442,9 @@ namespace Imp {
return p[i];
}
B* p;
size_type* j;
size_type i;
B* p = nullptr;
size_type* j = nullptr;
size_type i = 0;
};
/** @brief The iterator type. */
......
dune/istl/bccsmatrix.hh 0 → 100644
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_BCCSMATRIX_HH
#define DUNE_ISTL_BCCSMATRIX_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
namespace Dune::ISTL::Impl
{
/**
* @brief A block matrix with compressed-column storage
*
* @tparam B The type of a matrix entry (may be a matrix itself)
* @tparam I The type used for row and column indices.
* When using BCCSMatrix to set up a SuiteSparse solver, this type
* should be set to the index type used by the solver. That way,
* the solver can use the row/column information without any copying.
*
* Currently this class is mainly a vehicle to get matrices into the solvers
* of the SuiteSparse package. The class tries to follow the dune-istl
* matrix interface, but it doesn't implement the interface entirely.
*/
template<class B, class I = typename std::allocator<B>::size_type>
class BCCSMatrix
{
public:
using Index = I;
using size_type = std::size_t;
/** \brief Default constructor
*/
BCCSMatrix()
: N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
{}
/** @brief Destructor */
~BCCSMatrix()
{
if(N_+M_+Nnz_!=0)
free();
}
/** \brief Set matrix size */
void setSize(size_type rows, size_type columns)
{
N_ = rows;
M_ = columns;
}
/**
* @brief Get the number of rows.
* @return The number of rows.
*/
size_type N() const
{
return N_;
}
/** \brief Return the number of nonzero block entries */
size_type nonzeroes() const
{
return Nnz_;
}
/**
* @brief Get the number of columns.
* @return The number of columns.
*/
size_type M() const
{
return M_;
}
/** \brief Direct access to the array of matrix entries
*
* This is meant to be handed directly to SuiteSparse solvers.
* Note that the 'Index' type is controllable via the second
* class template argument.
*/
B* getValues() const
{
return values;
}
/** \brief Direct access to the array of row indices
*
* This is meant to be handed directly to SuiteSparse solvers.
* Note that the 'Index' type is controllable via the second
* class template argument.
*/
Index* getRowIndex() const
{
return rowindex;
}
/** \brief Direct access to the array of column entry points
*
* This is meant to be handed directly to SuiteSparse solvers.
* Note that the 'Index' type is controllable via the second
* class template argument.
*/
Index* getColStart() const
{
return colstart;
}
/** \brief Assignment operator */
BCCSMatrix& operator=(const BCCSMatrix& mat)
{
if(N_+M_+Nnz_!=0)
free();
N_=mat.N_;
M_=mat.M_;
Nnz_= mat.Nnz_;
if(M_>0) {
colstart=new size_type[M_+1];
for(size_type i=0; i<=M_; ++i)
colstart[i]=mat.colstart[i];
}
if(Nnz_>0) {
values = new B[Nnz_];
rowindex = new size_type[Nnz_];
for(size_type i=0; i<Nnz_; ++i)
values[i]=mat.values[i];
for(size_type i=0; i<Nnz_; ++i)
rowindex[i]=mat.rowindex[i];
}
return *this;
}
/** @brief free allocated space. */
virtual void free()
{
delete[] values;
delete[] rowindex;
delete[] colstart;
N_ = 0;
M_ = 0;
Nnz_ = 0;
}
public:
size_type N_, M_, Nnz_;
B* values;
Index* rowindex;
Index* colstart;
};
}
#endif
dune/istl/colcompmatrix.hh → dune/istl/bccsmatrixinitializer.hh
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_COLCOMPMATRIX_HH
#define DUNE_ISTL_COLCOMPMATRIX_HH
#include "bcrsmatrix.hh"
#include "bvector.hh"
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/unused.hh>
#ifndef DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
#define DUNE_ISTL_BCCSMATRIX_INITIALIZER_HH
#include <limits>
#include <set>
namespace Dune
{
/**
* @brief Provides access to an iterator over all matrix rows.
*
* @tparam M The type of the matrix.
*/
template<class M>
class MatrixRowSet
{
public:
//! @brief The type of the matrix.
typedef M Matrix;
//! @brief The matrix row iterator type.
typedef typename Matrix::ConstRowIterator const_iterator;
#include <dune/common/typetraits.hh>
#include <dune/common/scalarmatrixview.hh>
/**
* @brief Construct an row set over all matrix rows.
* @param m The matrix for which we manage the rows.
*/
MatrixRowSet(const Matrix& m)
: m_(m)
{}
#include <dune/istl/bccsmatrix.hh>
//! @brief Get the row iterator at the first row.
const_iterator begin() const
{
return m_.begin();
}
//! @brief Get the row iterator at the end of all rows.
const_iterator end() const
{
return m_.end();
}
private:
const Matrix& m_;
};
namespace Dune
{
template<class I, class S, class D>
class OverlappingSchwarzInitializer;
}
namespace Dune::ISTL::Impl
{
/**
* @brief Provides access to an iterator over an arbitrary subset
* of matrix rows.
......@@ -137,327 +110,178 @@ namespace Dune
};
/**
* @brief Utility class for converting an ISTL Matrix into a column-compressed matrix
* @tparam M the matrix type
*/
template<class M>
struct ColCompMatrix
{};
/**
* @brief Inititializer for the ColCompMatrix
* @brief Initializer for a BCCSMatrix,
* as needed by OverlappingSchwarz
* @tparam M the matrix type
*/
template<class M>
struct ColCompMatrixInitializer
{};
template<class M, class X, class TM, class TD, class T1>
class SeqOverlappingSchwarz;
template<class T, bool flag>
struct SeqOverlappingSchwarzAssemblerHelper;
/**
* @brief Converter for BCRSMatrix to column-compressed Matrix.
* specialization for BCRSMatrix
* @tparam M The matrix type
* @tparam I The type used for row and column indices
*/
template<class B, class TA, int n, int m>
class ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
template<class M, class I = typename M::size_type>
class BCCSMatrixInitializer
{
friend struct ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<B,n,m>,TA> >;
template<class IList, class S, class D>
friend class Dune::OverlappingSchwarzInitializer;
public:
/** @brief The type of the matrix to convert. */
typedef BCRSMatrix<FieldMatrix<B,n,m>,TA> Matrix;
friend struct SeqOverlappingSchwarzAssemblerHelper<ColCompMatrix<Matrix>, true>;
using Matrix = M;
using Index = I;
typedef Dune::ISTL::Impl::BCCSMatrix<typename Matrix::field_type, I> OutputMatrix;
typedef typename Matrix::size_type size_type;
/**
* @brief Constructor that initializes the data.
* @param mat The matrix to convert.
/** \brief Constructor for dense matrix-valued matrices
*/
explicit ColCompMatrix(const Matrix& mat);
BCCSMatrixInitializer(OutputMatrix& mat_)
: mat(&mat_), cols(mat_.M())
{
if constexpr (Dune::IsNumber<typename M::block_type>::value)
{
n = m = 1;
}
else
{
// WARNING: This assumes that all blocks are dense and identical
n = M::block_type::rows;
m = M::block_type::cols;
}
ColCompMatrix();
mat->Nnz_=0;
}
/** @brief Destructor */
virtual ~ColCompMatrix();
BCCSMatrixInitializer()
: mat(0), cols(0), n(0), m(0)
{}
/**
* @brief Get the number of rows.
* @return The number of rows.
*/
size_type N() const
virtual ~BCCSMatrixInitializer()
{}
template<typename Iter>
void addRowNnz(const Iter& row) const
{
return N_;
mat->Nnz_+=row->getsize();
}
size_type nnz() const
template<typename Iter, typename FullMatrixIndex>
void addRowNnz(const Iter& row, const std::set<FullMatrixIndex>& indices) const
{
return Nnz_/n/m;
auto siter =indices.begin();
for (auto entry=row->begin(); entry!=row->end(); ++entry)
{
for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
if(siter==indices.end())
break;
if(*siter==entry.index())
// index is in subdomain
++mat->Nnz_;
}
}
/**
* @brief Get the number of columns.
* @return The number of columns.
*/
size_type M() const
template<typename Iter, typename SubMatrixIndex>
void addRowNnz(const Iter& row, const std::vector<SubMatrixIndex>& indices) const
{
return M_;
for (auto entry=row->begin(); entry!=row->end(); ++entry)
if (indices[entry.index()]!=std::numeric_limits<SubMatrixIndex>::max())
++mat->Nnz_;
}
B* getValues() const
void allocate()
{
return values;
allocateMatrixStorage();
allocateMarker();
}
int* getRowIndex() const
template<typename Iter, typename CIter>
void countEntries([[maybe_unused]] const Iter& row, const CIter& col) const
{
return rowindex;
countEntries(col.index());
}
int* getColStart() const
void countEntries(size_type colindex) const
{
return colstart;
for(size_type i=0; i < m; ++i)
{
assert(colindex*m+i<cols);
marker[colindex*m+i]+=n;
}
}
ColCompMatrix& operator=(const Matrix& mat);
ColCompMatrix& operator=(const ColCompMatrix& mat);
/**
* @brief Initialize data from a given set of matrix rows and columns
* @param mat the matrix with the values
* @param mrs The set of row (and column) indices to remove
*/
virtual void setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs);
/** @brief free allocated space. */
virtual void free();
/** @brief Initialize data from given matrix. */
virtual void setMatrix(const Matrix& mat);
public:
int N_, M_, Nnz_;
B* values;
int* rowindex;
int* colstart;
};
template<class T, class A, int n, int m>
class ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >
{
template<class I, class S, class D>
friend class OverlappingSchwarzInitializer;
public:
typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
typedef Dune::ColCompMatrix<Matrix> ColCompMatrix;
typedef typename Matrix::row_type::const_iterator CIter;
typedef typename Matrix::size_type size_type;
ColCompMatrixInitializer(ColCompMatrix& lum);
ColCompMatrixInitializer();
virtual ~ColCompMatrixInitializer();
template<typename Iter>
void addRowNnz(const Iter& row) const;
template<typename Iter, typename FullMatrixIndex>
void addRowNnz(const Iter& row, const std::set<FullMatrixIndex>& indices) const;
template<typename Iter, typename SubMatrixIndex>
void addRowNnz(const Iter& row, const std::vector<SubMatrixIndex>& indices) const;
void allocate();
template<typename Iter>
void countEntries(const Iter& row, const CIter& col) const;
void countEntries(size_type colidx) const;
void calcColstart() const;
template<typename Iter>
void copyValue(const Iter& row, const CIter& col) const;
void copyValue(const CIter& col, size_type rowindex, size_type colidx) const;
virtual void createMatrix() const;
protected:
void allocateMatrixStorage() const;
void allocateMarker();
ColCompMatrix* mat;
size_type cols;
mutable size_type *marker;
};
template<class T, class A, int n, int m>
ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::ColCompMatrixInitializer(ColCompMatrix& mat_)
: mat(&mat_), cols(mat_.M()), marker(0)
{
mat->Nnz_=0;
}
template<class T, class A, int n, int m>
ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::ColCompMatrixInitializer()
: mat(0), cols(0), marker(0)
{}
template<class T, class A, int n, int m>
ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::~ColCompMatrixInitializer()
{
if(marker)
delete[] marker;
}
template<class T, class A, int n, int m>
template<typename Iter>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row) const
{
mat->Nnz_+=row->getsize();
}
template<class T, class A, int n, int m>
template<typename Iter, typename FullMatrixIndex>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row,
const std::set<FullMatrixIndex>& indices) const
{
typedef typename Iter::value_type::const_iterator RIter;
typedef typename std::set<FullMatrixIndex>::const_iterator MIter;
MIter siter =indices.begin();
for(RIter entry=row->begin(); entry!=row->end(); ++entry)
void calcColstart() const
{
for(; siter!=indices.end() && *siter<entry.index(); ++siter) ;
if(siter==indices.end())
break;
if(*siter==entry.index())
// index is in subdomain
++mat->Nnz_;
mat->colstart[0]=0;
for(size_type i=0; i < cols; ++i) {
assert(i<cols);
mat->colstart[i+1]=mat->colstart[i]+marker[i];
marker[i]=mat->colstart[i];
}
}
}
template<class T, class A, int n, int m>
template<typename Iter, typename SubMatrixIndex>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row,
const std::vector<SubMatrixIndex>& indices) const
{
using RIter = typename Iter::value_type::const_iterator;
for(RIter entry=row->begin(); entry!=row->end(); ++entry)
if (indices[entry.index()]!=std::numeric_limits<SubMatrixIndex>::max())
++mat->Nnz_;
}
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocate()
{
allocateMatrixStorage();
allocateMarker();
}
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage() const
{
mat->Nnz_*=n*m;
// initialize data
mat->values=new T[mat->Nnz_];
mat->rowindex=new int[mat->Nnz_];
mat->colstart=new int[cols+1];
}
template<typename Iter, typename CIter>
void copyValue(const Iter& row, const CIter& col) const
{
copyValue(col, row.index(), col.index());
}
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMarker()
{
marker = new typename Matrix::size_type[cols];
template<typename CIter>
void copyValue(const CIter& col, size_type rowindex, size_type colindex) const
{
for(size_type i=0; i<n; i++) {
for(size_type j=0; j<m; j++) {
assert(colindex*m+j<cols-1 || (size_type)marker[colindex*m+j]<(size_type)mat->colstart[colindex*m+j+1]);
assert((size_type)marker[colindex*m+j]<mat->Nnz_);
mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
mat->values[marker[colindex*m+j]] = Dune::Impl::asMatrix(*col)[i][j];
++marker[colindex*m+j]; // index for next entry in column
}
}
}
for(size_type i=0; i < cols; ++i)
marker[i]=0;
}
virtual void createMatrix() const
{
marker.clear();
}
template<class T, class A, int n, int m>
template<typename Iter>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(const Iter& row, const CIter& col) const
{
DUNE_UNUSED_PARAMETER(row);
countEntries(col.index());
}
protected:
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(size_type colindex) const
{
for(size_type i=0; i < m; ++i)
void allocateMatrixStorage() const
{
assert(colindex*m+i<cols);
marker[colindex*m+i]+=n;
mat->Nnz_*=n*m;
// initialize data
mat->values=new typename M::field_type[mat->Nnz_];
mat->rowindex=new I[mat->Nnz_];
mat->colstart=new I[cols+1];
}
}
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::calcColstart() const
{
mat->colstart[0]=0;
for(size_type i=0; i < cols; ++i) {
assert(i<cols);
mat->colstart[i+1]=mat->colstart[i]+marker[i];
marker[i]=mat->colstart[i];
void allocateMarker()
{
marker.resize(cols);
std::fill(marker.begin(), marker.end(), 0);
}
}
template<class T, class A, int n, int m>
template<typename Iter>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const Iter& row, const CIter& col) const
{
copyValue(col, row.index(), col.index());
}
OutputMatrix* mat;
size_type cols;
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const CIter& col, size_type rowindex, size_type colindex) const
{
for(size_type i=0; i<n; i++) {
for(size_type j=0; j<m; j++) {
assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
assert((int)marker[colindex*m+j]<mat->Nnz_);
mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
mat->values[marker[colindex*m+j]]=(*col)[i][j];
++marker[colindex*m+j]; // index for next entry in column
}
}
}
// Number of rows/columns of the matrix entries
// (assumed to be scalars or dense matrices)
size_type n, m;
template<class T, class A, int n, int m>
void ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::createMatrix() const
{
delete[] marker;
marker=0;
}
mutable std::vector<size_type> marker;
};
template<class F, class MRS>
void copyToColCompMatrix(F& initializer, const MRS& mrs)
template<class F, class Matrix>
void copyToBCCSMatrix(F& initializer, const Matrix& matrix)
{
typedef typename MRS::const_iterator Iter;
typedef typename std::iterator_traits<Iter>::value_type::const_iterator CIter;
for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
for (auto row=matrix.begin(); row!= matrix.end(); ++row)
initializer.addRowNnz(row);
initializer.allocate();
for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
for (auto row=matrix.begin(); row!= matrix.end(); ++row) {
for(CIter col=row->begin(); col != row->end(); ++col)
for (auto col=row->begin(); col != row->end(); ++col)
initializer.countEntries(row, col);
}
initializer.calcColstart();
for(Iter row=mrs.begin(); row!= mrs.end(); ++row) {
for(CIter col=row->begin(); col != row->end(); ++col) {
for (auto row=matrix.begin(); row!= matrix.end(); ++row) {
for (auto col=row->begin(); col != row->end(); ++col) {
initializer.copyValue(row, col);
}
......@@ -466,7 +290,7 @@ namespace Dune
}
template<class F, class M,class S>
void copyToColCompMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
void copyToBCCSMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
{
typedef MatrixRowSubset<M,S> MRS;
typedef typename MRS::RowIndexSet SIS;
......@@ -511,100 +335,5 @@ namespace Dune
initializer.createMatrix();
}
#ifndef DOXYGEN
template<class B, class TA, int n, int m>
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::ColCompMatrix()
: N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
{}
template<class B, class TA, int n, int m>
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
::ColCompMatrix(const Matrix& mat)
: N_(n*mat.N()), M_(m*mat.M()), Nnz_(n*m*mat.nonzeroes())
{}
template<class B, class TA, int n, int m>
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const Matrix& mat)
{
if(N_+M_+Nnz_!=0)
free();
setMatrix(mat);
return *this;
}
template<class B, class TA, int n, int m>
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const ColCompMatrix& mat)
{
if(N_+M_+Nnz_!=0)
free();
N_=mat.N_;
M_=mat.M_;
Nnz_= mat.Nnz_;
if(M_>0) {
colstart=new int[M_+1];
for(int i=0; i<=M_; ++i)
colstart[i]=mat.colstart[i];
}
if(Nnz_>0) {
values = new B[Nnz_];
rowindex = new int[Nnz_];
for(int i=0; i<Nnz_; ++i)
values[i]=mat.values[i];
for(int i=0; i<Nnz_; ++i)
rowindex[i]=mat.rowindex[i];
}
return *this;
}
template<class B, class TA, int n, int m>
void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
::setMatrix(const Matrix& mat)
{
N_=n*mat.N();
M_=m*mat.M();
ColCompMatrixInitializer<Matrix> initializer(*this);
copyToColCompMatrix(initializer, MatrixRowSet<Matrix>(mat));
}
template<class B, class TA, int n, int m>
void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
::setMatrix(const Matrix& mat, const std::set<std::size_t>& mrs)
{
if(N_+M_+Nnz_!=0)
free();
N_=mrs.size()*n;
M_=mrs.size()*m;
ColCompMatrixInitializer<Matrix> initializer(*this);
copyToColCompMatrix(initializer, MatrixRowSubset<Matrix,std::set<std::size_t> >(mat,mrs));
}
template<class B, class TA, int n, int m>
ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::~ColCompMatrix()
{
if(N_+M_+Nnz_!=0)
free();
}
template<class B, class TA, int n, int m>
void ColCompMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::free()
{
delete[] values;
delete[] rowindex;
delete[] colstart;
N_ = 0;
M_ = 0;
Nnz_ = 0;
}
#endif // DOXYGEN
}
#endif
dune/istl/bcrsmatrix.hh
View file @ bc1f69fc
This diff is collapsed.
dune/istl/bdmatrix.hh
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_BDMATRIX_HH
......@@ -6,8 +8,10 @@
#include <memory>
#include <dune/common/rangeutilities.hh>
#include <dune/common/scalarmatrixview.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/blocklevel.hh>
/** \file
\author Oliver Sander
......@@ -46,9 +50,6 @@ namespace Dune {
//! The type for the index access and the size
typedef typename A::size_type size_type;
//! increment block level counter
static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
/** \brief Default constructor */
BDMatrix() : BCRSMatrix<B,A>() {}
......@@ -106,17 +107,25 @@ namespace Dune {
return *this;
}
/** \brief Solve the system Ax=b in O(n) time
*
* \exception ISTLError if the matrix is singular
*
*/
template <class V>
void solve (V& x, const V& rhs) const {
for (size_type i=0; i<this->N(); i++)
{
auto&& xv = Impl::asVector(x[i]);
auto&& rhsv = Impl::asVector(rhs[i]);
Impl::asMatrix((*this)[i][i]).solve(xv,rhsv);
}
}
/** \brief Inverts the matrix */
void invert() {
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (size_type i=0; i<this->N(); i++)
(*this)[i][i] = 1.0 / id((*this)[i][i]);
},
[&](auto id) {
for (size_type i=0; i<this->N(); i++)
id((*this)[i][i]).invert();
});
for (size_type i=0; i<this->N(); i++)
Impl::asMatrix((*this)[i][i]).invert();
}
private:
......@@ -134,6 +143,13 @@ namespace Dune {
void addindex (size_type row, size_type col) {}
void endindices () {}
};
template<typename B, typename A>
struct FieldTraits< BDMatrix<B, A> >
{
using field_type = typename BDMatrix<B, A>::field_type;
using real_type = typename FieldTraits<field_type>::real_type;
};
/** @}*/
} // end namespace Dune
......
dune/istl/blocklevel.hh 0 → 100644
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_BLOCKLEVEL_HH
#define DUNE_ISTL_BLOCKLEVEL_HH
#include <algorithm>
#include <type_traits>
#include <dune/common/indices.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/hybridutilities.hh>
/*!
* \file
* \brief Helper functions for determining the vector/matrix block level
*/
// forward declaration
namespace Dune {
template<typename... Args>
class MultiTypeBlockVector;
template<typename FirstRow, typename... Args>
class MultiTypeBlockMatrix;
} // end namespace Dune
namespace Dune { namespace Impl {
// forward declaration
template<typename T> struct MaxBlockLevel;
template<typename T> struct MinBlockLevel;
//! recursively determine block level of a MultiTypeBlockMatrix
template<typename M, template<typename B> typename BlockLevel, typename Op>
constexpr std::size_t blockLevelMultiTypeBlockMatrix(const Op& op)
{
// inialize with zeroth diagonal block
using namespace Dune::Indices;
using Block00 = typename std::decay_t<decltype(std::declval<M>()[_0][_0])>;
std::size_t blockLevel = BlockLevel<Block00>::value() + 1;
// iterate over all blocks to determine min/max block level
using namespace Dune::Hybrid;
forEach(integralRange(index_constant<M::N()>()), [&](auto&& i) {
using namespace Dune::Hybrid; // needed for icc, see issue #31
forEach(integralRange(index_constant<M::M()>()), [&](auto&& j) {
using Block = typename std::decay_t<decltype(std::declval<M>()[i][j])>;
blockLevel = op(blockLevel, BlockLevel<Block>::value() + 1);
});
});
return blockLevel;
}
//! recursively determine block level of a MultiTypeBlockVector
template<typename V, template<typename B> typename BlockLevel, typename Op>
constexpr std::size_t blockLevelMultiTypeBlockVector(const Op& op)
{
// inialize with zeroth block
using namespace Dune::Indices;
using Block0 = typename std::decay_t<decltype(std::declval<V>()[_0])>;
std::size_t blockLevel = BlockLevel<Block0>::value() + 1;
// iterate over all blocks to determine min/max block level
using namespace Dune::Hybrid;
forEach(integralRange(index_constant<V::size()>()), [&](auto&& i) {
using Block = typename std::decay_t<decltype(std::declval<V>()[i])>;
blockLevel = op(blockLevel, BlockLevel<Block>::value() + 1);
});
return blockLevel;
}
template<typename T>
struct MaxBlockLevel
{
static constexpr std::size_t value(){
if constexpr (IsNumber<T>::value)
return 0;
else
return MaxBlockLevel<typename T::block_type>::value() + 1;
}
};
template<typename T>
struct MinBlockLevel
{
// the default implementation assumes minBlockLevel == maxBlockLevel
static constexpr std::size_t value()
{ return MaxBlockLevel<T>::value(); }
};
// max block level for MultiTypeBlockMatrix
template<typename FirstRow, typename... Args>
struct MaxBlockLevel<Dune::MultiTypeBlockMatrix<FirstRow, Args...>>
{
static constexpr std::size_t value()
{
using M = MultiTypeBlockMatrix<FirstRow, Args...>;
constexpr auto max = [](const auto& a, const auto& b){ return std::max(a,b); };
return blockLevelMultiTypeBlockMatrix<M, MaxBlockLevel>(max);
}
};
// min block level for MultiTypeBlockMatrix
template<typename FirstRow, typename... Args>
struct MinBlockLevel<Dune::MultiTypeBlockMatrix<FirstRow, Args...>>
{
static constexpr std::size_t value()
{
using M = MultiTypeBlockMatrix<FirstRow, Args...>;
constexpr auto min = [](const auto& a, const auto& b){ return std::min(a,b); };
return blockLevelMultiTypeBlockMatrix<M, MinBlockLevel>(min);
}
};
// max block level for MultiTypeBlockVector
template<typename... Args>
struct MaxBlockLevel<Dune::MultiTypeBlockVector<Args...>>
{
static constexpr std::size_t value()
{
using V = MultiTypeBlockVector<Args...>;
constexpr auto max = [](const auto& a, const auto& b){ return std::max(a,b); };
return blockLevelMultiTypeBlockVector<V, MaxBlockLevel>(max);
}
};
// min block level for MultiTypeBlockVector
template<typename... Args>
struct MinBlockLevel<Dune::MultiTypeBlockVector<Args...>>
{
static constexpr std::size_t value()
{
using V = MultiTypeBlockVector<Args...>;
constexpr auto min = [](const auto& a, const auto& b){ return std::min(a,b); };
return blockLevelMultiTypeBlockVector<V, MinBlockLevel>(min);
}
};
// special case: empty MultiTypeBlockVector
template<>
struct MaxBlockLevel<Dune::MultiTypeBlockVector<>>
{
static constexpr std::size_t value()
{ return 0; };
};
// special case: empty MultiTypeBlockVector
template<>
struct MinBlockLevel<Dune::MultiTypeBlockVector<>>
{
static constexpr std::size_t value()
{ return 0; };
};
}} // end namespace Dune::Impl
namespace Dune {
//! Determine the maximum block level of a possibly nested vector/matrix type
template<typename T>
constexpr std::size_t maxBlockLevel()
{ return Impl::MaxBlockLevel<T>::value(); }
//! Determine the minimum block level of a possibly nested vector/matrix type
template<typename T>
constexpr std::size_t minBlockLevel()
{ return Impl::MinBlockLevel<T>::value(); }
//! Determine if a vector/matrix has a uniquely determinable block level
template<typename T>
constexpr bool hasUniqueBlockLevel()
{ return maxBlockLevel<T>() == minBlockLevel<T>(); }
//! Determine the block level of a possibly nested vector/matrix type
template<typename T>
constexpr std::size_t blockLevel()
{
static_assert(hasUniqueBlockLevel<T>(), "Block level cannot be uniquely determined!");
return Impl::MaxBlockLevel<T>::value();
}
} // end namespace Dune
#endif
dune/istl/btdmatrix.hh
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_BTDMATRIX_HH
#define DUNE_ISTL_BTDMATRIX_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/blocklevel.hh>
/** \file
\author Oliver Sander
......@@ -43,9 +48,6 @@ namespace Dune {
//! The type for the index access and the size
typedef typename A::size_type size_type;
//! increment block level counter
static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
/** \brief Default constructor */
BTDMatrix() : BCRSMatrix<B,A>() {}
......@@ -120,13 +122,9 @@ namespace Dune {
// special handling for 1x1 matrices. The generic algorithm doesn't work for them
if (this->N()==1) {
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
x[0] = id(rhs[0]) / id((*this)[0][0]);
},
[&](auto id) {
id(*this)[0][0].solve(x[0],rhs[0]);
});
auto&& x0 = Impl::asVector(x[0]);
auto&& rhs0 = Impl::asVector(rhs[0]);
Impl::asMatrix((*this)[0][0]).solve(x0, rhs0);
return;
}
......@@ -138,95 +136,48 @@ namespace Dune {
/* Modify the coefficients. */
block_type a_00_inv = (*this)[0][0];
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
a_00_inv = 1.0 / id(a_00_inv);
},
[&](auto id) {
id(a_00_inv).invert();
});
//c[0] /= (*this)[0][0]; /* Division by zero risk. */
block_type c_0_tmp = c[0];
Hybrid::ifElse(IsNumber<B>(), /* Division by zero risk. */
[&](auto id) {
c[0] = a_00_inv * id(c_0_tmp);
},
[&](auto id) {
FMatrixHelp::multMatrix(id(a_00_inv), id(c_0_tmp), id(c[0]));
});
Impl::asMatrix(a_00_inv).invert();
//c[0] /= (*this)[0][0]; /* Division by zero risk. */
block_type tmp = a_00_inv;
Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[0]));
c[0] = tmp;
// d = a^{-1} d /* Division by zero would imply a singular matrix. */
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
d[0] *= id(a_00_inv);
},
[&](auto id) {
auto d_0_tmp = d[0];
id(a_00_inv).mv(d_0_tmp,d[0]);
});
auto d_0_tmp = d[0];
auto&& d_0 = Impl::asVector(d[0]);
Impl::asMatrix(a_00_inv).mv(Impl::asVector(d_0_tmp),d_0);
for (unsigned int i = 1; i < this->N(); i++) {
// id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
block_type tmp;
Hybrid::ifElse(IsNumber<B>(),
[&](auto metaId) {
tmp = metaId((*this)[i][i-1]) * metaId(c[i-1]);
},
[&](auto metaId) {
FMatrixHelp::multMatrix(metaId((*this)[i][i-1]), metaId(c[i-1]), metaId(tmp));
});
tmp = (*this)[i][i-1];
Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[i-1]));
block_type id = (*this)[i][i];
id -= tmp;
Hybrid::ifElse(IsNumber<B>(), /* Division by zero risk. */
[&](auto metaId) {
id = 1.0 / metaId(id);
},
[&](auto metaId) {
metaId(id).invert();
});
Impl::asMatrix(id).invert(); /* Division by zero risk. */
if (i<c.size()) {
// c[i] *= id
Hybrid::ifElse(IsNumber<B>(), /* Last value calculated is redundant. */
[&](auto metaId) {
c[i] *= metaId(id);
},
[&](auto metaId) {
tmp = c[i];
FMatrixHelp::multMatrix(metaId(id), metaId(tmp), metaId(c[i]));
});
Impl::asMatrix(c[i]).leftmultiply(Impl::asMatrix(id)); /* Last value calculated is redundant. */
}
// d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
Hybrid::ifElse(IsNumber<B>(),
[&](auto metaId) {
d[i] -= (*this)[i][i-1] * metaId(d[i-1]);
d[i] *= metaId(id);
},
[&](auto metaId) {
metaId((*this)[i][i-1]).mmv(d[i-1], d[i]);
auto tmpVec = d[i];
metaId(id).mv(tmpVec, d[i]);
});
auto&& d_i = Impl::asVector(d[i]);
Impl::asMatrix((*this)[i][i-1]).mmv(Impl::asVector(d[i-1]), d_i);
auto tmpVec = d[i];
Impl::asMatrix(id).mv(Impl::asVector(tmpVec), d_i);
}
/* Now back substitute. */
x[this->N() - 1] = d[this->N() - 1];
Hybrid::ifElse(IsNumber<B>(),
[&](auto metaId) {
for (int i = this->N() - 2; i >= 0; i--)
x[i] = d[i] - c[i] * metaId(x[i+1]);
},
[&](auto metaId) {
for (int i = this->N() - 2; i >= 0; i--) {
//x[i] = d[i] - c[i] * x[i + 1];
x[i] = d[i];
metaId(c[i]).mmv(x[i+1], x[i]);
}
});
for (int i = this->N() - 2; i >= 0; i--) {
//x[i] = d[i] - c[i] * x[i + 1];
x[i] = d[i];
auto&& x_i = Impl::asVector(x[i]);
Impl::asMatrix(c[i]).mmv(Impl::asVector(x[i+1]), x_i);
}
}
......@@ -245,6 +196,14 @@ namespace Dune {
void addindex (size_type row, size_type col) {}
void endindices () {}
};
template<typename B, typename A>
struct FieldTraits< BTDMatrix<B, A> >
{
using field_type = typename BTDMatrix<B, A>::field_type;
using real_type = typename FieldTraits<field_type>::real_type;
};
/** @}*/
} // end namespace Dune
......
dune/istl/bvector.hh
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
......@@ -13,13 +15,15 @@
#include <utility>
#include <vector>
#include <dune/common/deprecated.hh>
#include <dune/common/dotproduct.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/unused.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/istl/blocklevel.hh>
#include "basearray.hh"
#include "istlexception.hh"
......@@ -49,11 +53,6 @@ namespace Imp {
{
public:
using field_type = B;
static constexpr unsigned int blockLevel()
{
return 0;
}
};
template <class B>
......@@ -61,11 +60,6 @@ namespace Imp {
{
public:
using field_type = typename B::field_type;
static constexpr unsigned int blockLevel()
{
return B::blocklevel;
}
};
/** \brief Define some derived types transparently for number types and dune-istl matrix/vector types
......@@ -86,8 +80,8 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
class block_vector_unmanaged : public base_array_unmanaged<B,A>
template<class B, class ST=std::size_t >
class block_vector_unmanaged : public base_array_unmanaged<B,ST>
{
public:
......@@ -97,17 +91,14 @@ namespace Imp {
//! export the type representing the components
typedef B block_type;
//! export the allocator type
typedef A allocator_type;
//! The size type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! make iterators available as types
typedef typename base_array_unmanaged<B,A>::iterator Iterator;
typedef typename base_array_unmanaged<B,ST>::iterator Iterator;
//! make iterators available as types
typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
typedef typename base_array_unmanaged<B,ST>::const_iterator ConstIterator;
//! for STL compatibility
typedef B value_type;
......@@ -169,15 +160,8 @@ namespace Imp {
#ifdef DUNE_ISTL_WITH_CHECKING
if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
(*this)[i] += a*id(y[i]);
},
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
id((*this)[i]).axpy(a,id(y[i]));
});
for (size_type i=0; i<this->n; ++i)
Impl::asVector((*this)[i]).axpy(a,Impl::asVector(y[i]));
return *this;
}
......@@ -190,8 +174,8 @@ namespace Imp {
* @param y other (compatible) vector
* @return
*/
template<class OtherB, class OtherA>
auto operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
template<class OtherB, class OtherST>
auto operator* (const block_vector_unmanaged<OtherB,OtherST>& y) const
{
typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
PromotedType sum(0);
......@@ -211,8 +195,8 @@ namespace Imp {
* @param y other (compatible) vector
* @return
*/
template<class OtherB, class OtherA>
auto dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
template<class OtherB, class OtherST>
auto dot(const block_vector_unmanaged<OtherB,OtherST>& y) const
{
typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
PromotedType sum(0);
......@@ -220,16 +204,8 @@ namespace Imp {
if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
#endif
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
using namespace std;
for (size_type i=0; i<this->n; ++i)
sum += Dune::dot(id((*this)[i]),y[i]);
},
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
sum += id((*this)[i]).dot(y[i]);
});
for (size_type i=0; i<this->n; ++i)
sum += Impl::asVector((*this)[i]).dot(Impl::asVector(y[i]));
return sum;
}
......@@ -240,16 +216,8 @@ namespace Imp {
typename FieldTraits<field_type>::real_type one_norm () const
{
typename FieldTraits<field_type>::real_type sum=0;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
using namespace std;
for (size_type i=0; i<this->n; ++i)
sum += abs(id((*this)[i]));
},
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
sum += id((*this)[i]).one_norm();
});
for (size_type i=0; i<this->n; ++i)
sum += Impl::asVector((*this)[i]).one_norm();
return sum;
}
......@@ -257,15 +225,8 @@ namespace Imp {
typename FieldTraits<field_type>::real_type one_norm_real () const
{
typename FieldTraits<field_type>::real_type sum=0;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
sum += fvmeta::absreal((*this)[i]);
},
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
sum += id((*this)[i]).one_norm_real();
});
for (size_type i=0; i<this->n; ++i)
sum += Impl::asVector((*this)[i]).one_norm_real();
return sum;
}
......@@ -280,15 +241,8 @@ namespace Imp {
typename FieldTraits<field_type>::real_type two_norm2 () const
{
typename FieldTraits<field_type>::real_type sum=0;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
sum += fvmeta::abs2(id((*this)[i]));
},
[&](auto id) {
for (size_type i=0; i<this->n; ++i)
sum += id((*this)[i]).two_norm2();
});
for (size_type i=0; i<this->n; ++i)
sum += Impl::asVector((*this)[i]).two_norm2();
return sum;
}
......@@ -300,20 +254,10 @@ namespace Imp {
using std::max;
real_type norm = 0;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (auto const &x : *this) {
using std::abs;
real_type const a = abs(x);
norm = max(a, norm);
}
},
[&](auto id) {
for (auto const &x : *this) {
real_type const a = x.infinity_norm();
norm = max(a, norm);
}
});
for (auto const &xi : *this) {
real_type const a = Impl::asVector(xi).infinity_norm();
norm = max(a, norm);
}
return norm;
}
......@@ -325,19 +269,10 @@ namespace Imp {
using std::max;
real_type norm = 0;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (auto const &x : *this) {
real_type const a = fvmeta::absreal(x);
norm = max(a, norm);
}
},
[&](auto id) {
for (auto const &x : *this) {
real_type const a = x.infinity_norm_real();
norm = max(a, norm);
}
});
for (auto const &xi : *this) {
real_type const a = Impl::asVector(xi).infinity_norm_real();
norm = max(a, norm);
}
return norm;
}
......@@ -352,23 +287,11 @@ namespace Imp {
real_type norm = 0;
real_type isNaN = 1;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (auto const &x : *this) {
using std::abs;
real_type const a = abs(x);
norm = max(a, norm);
isNaN += a;
}
},
[&](auto id) {
for (auto const &x : *this) {
real_type const a = id(x).infinity_norm();
norm = max(a, norm);
isNaN += a;
}
});
for (auto const &xi : *this) {
real_type const a = Impl::asVector(xi).infinity_norm();
norm = max(a, norm);
isNaN += a;
}
return norm * (isNaN / isNaN);
}
......@@ -382,21 +305,11 @@ namespace Imp {
real_type norm = 0;
real_type isNaN = 1;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
for (auto const &x : *this) {
real_type const a = fvmeta::absreal(x);
norm = max(a, norm);
isNaN += a;
}
},
[&](auto id) {
for (auto const &x : *this) {
real_type const a = id(x).infinity_norm_real();
norm = max(a, norm);
isNaN += a;
}
});
for (auto const &xi : *this) {
real_type const a = Impl::asVector(xi).infinity_norm_real();
norm = max(a, norm);
isNaN += a;
}
return norm * (isNaN / isNaN);
}
......@@ -414,21 +327,15 @@ namespace Imp {
{
size_type d=0;
Hybrid::ifElse(IsNumber<B>(),
[&](auto id) {
d = this->n;
},
[&](auto id) {
for (size_type i=0; i<this->n; i++)
d += (*this)[i].dim();
});
for (size_type i=0; i<this->n; i++)
d += Impl::asVector((*this)[i]).dim();
return d;
}
protected:
//! make constructor protected, so only derived classes can be instantiated
block_vector_unmanaged () : base_array_unmanaged<B,A>()
block_vector_unmanaged () : base_array_unmanaged<B,ST>()
{ }
};
......@@ -481,7 +388,7 @@ namespace Imp {
enables error checking.
*/
template<class B, class A=std::allocator<B> >
class BlockVector : public Imp::block_vector_unmanaged<B,A>
class BlockVector : public Imp::block_vector_unmanaged<B,typename A::size_type>
{
public:
......@@ -499,14 +406,11 @@ namespace Imp {
//! The type for the index access
typedef typename A::size_type size_type;
//! increment block level counter
static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
//===== constructors and such
......@@ -563,37 +467,11 @@ namespace Imp {
*/
void reserve(size_type capacity)
{
DUNE_UNUSED const auto &guard =
[[maybe_unused]] const auto &guard =
Imp::makeScopeGuard([this]{ syncBaseArray(); });
storage_.reserve(capacity);
}
/**
* @brief Reserve space.
*
* After calling this method the vector can hold up to
* capacity values. If the specified capacity is smaller
* than the current capacity and bigger than the current size
* space will be freed.
*
* If the template parameter copyOldValues is true the values will
* be copied. If it is false the old values are lost.
*
* @param capacity The maximum number of elements the vector
* needs to hold.
* @param copyOldValues Ignored, values are always copied.
*
* \deprecated This method is deprecated, since the extra copyOldValues
* parameter is unusual, and the previous implementation did
* not always reduce memory anyway.
*/
DUNE_DEPRECATED_MSG("Use the overload without the second parameter, "
"values are always copied")
void reserve(size_type capacity, bool copyOldValues)
{
reserve(capacity);
}
/**
* @brief Get the capacity of the vector.
*
......@@ -617,38 +495,11 @@ namespace Imp {
*/
void resize(size_type size)
{
DUNE_UNUSED const auto &guard =
[[maybe_unused]] const auto &guard =
Imp::makeScopeGuard([this]{ syncBaseArray(); });
storage_.resize(size);
}
/**
* @brief Resize the vector.
*
* After calling this method BlockVector::N() will return size
* If the capacity of the vector is smaller than the specified
* size then reserve(size) will be called.
*
* If the template parameter copyOldValues is true the values
* will be copied if the capacity changes. If it is false
* the old values are lost.
* @param size The new size of the vector.
* @param copyOldValues Ignored, values are always copied.
*
* \deprecated This method is deprecated, since the extra copyOldValues
* parameter is unusual and conflicts with the usual meaning
* (default value for newly created elements).
*/
DUNE_DEPRECATED_MSG("Use the overload without the second parameter, "
"values are always copied")
void resize(size_type size, bool copyOldValues)
{
resize(size);
}
//! copy constructor
BlockVector(const BlockVector &a)
noexcept(noexcept(std::declval<BlockVector>().storage_ = a.storage_))
......@@ -668,7 +519,7 @@ namespace Imp {
BlockVector& operator= (const BlockVector& a)
noexcept(noexcept(std::declval<BlockVector>().storage_ = a.storage_))
{
DUNE_UNUSED const auto &guard =
[[maybe_unused]] const auto &guard =
Imp::makeScopeGuard([this]{ syncBaseArray(); });
storage_ = a.storage_;
return *this;
......@@ -687,7 +538,7 @@ namespace Imp {
noexcept(noexcept(
std::declval<BlockVector&>().storage_.swap(other.storage_)))
{
DUNE_UNUSED const auto &guard = Imp::makeScopeGuard([&]{
[[maybe_unused]] const auto &guard = Imp::makeScopeGuard([&]{
syncBaseArray();
other.syncBaseArray();
});
......@@ -698,7 +549,7 @@ namespace Imp {
BlockVector& operator= (const field_type& k)
{
// forward to operator= in base class
(static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
(static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
return *this;
}
......@@ -765,7 +616,7 @@ namespace Imp {
#else
template<class B, class A=std::allocator<B> >
#endif
class BlockVectorWindow : public Imp::block_vector_unmanaged<B,A>
class BlockVectorWindow : public Imp::block_vector_unmanaged<B,typename A::size_type>
{
public:
......@@ -783,19 +634,16 @@ namespace Imp {
//! The type for the index access
typedef typename A::size_type size_type;
//! increment block level counter
static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
//===== constructors and such
//! makes empty array
BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
BlockVectorWindow () : Imp::block_vector_unmanaged<B,size_type>()
{ }
//! make array from given pointer and size
......@@ -831,7 +679,7 @@ namespace Imp {
//! assign from scalar
BlockVectorWindow& operator= (const field_type& k)
{
(static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
(static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
return *this;
}
......@@ -890,8 +738,8 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,A>
template<class B, class ST=std::size_t >
class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,ST>
{
public:
......@@ -903,17 +751,14 @@ namespace Imp {
//! export the type representing the components
typedef B block_type;
//! export the allocator type
typedef A allocator_type;
//! make iterators available as types
typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
typedef typename compressed_base_array_unmanaged<B,ST>::iterator Iterator;
//! make iterators available as types
typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
typedef typename compressed_base_array_unmanaged<B,ST>::const_iterator ConstIterator;
//! The type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//===== assignment from scalar
......@@ -956,7 +801,8 @@ namespace Imp {
#ifdef DUNE_ISTL_WITH_CHECKING
if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
#endif
for (size_type i=0; i<y.n; ++i) (this->operator[](y.j[i])).axpy(a,y.p[i]);
for (size_type i=0; i<y.n; ++i)
Impl::asVector((*this)[y.j[i]]).axpy(a,Impl::asVector(y.p[i]));
return *this;
}
......@@ -1109,7 +955,7 @@ namespace Imp {
protected:
//! make constructor protected, so only derived classes can be instantiated
compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,ST>()
{ }
//! return true if index sets coincide
......@@ -1143,8 +989,8 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,A>
template<class B, class ST=std::size_t >
class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,ST>
{
public:
......@@ -1156,25 +1002,19 @@ namespace Imp {
//! export the type representing the components
typedef B block_type;
//! export the allocator type
typedef A allocator_type;
//! The type for the index access
typedef typename A::size_type size_type;
//! increment block level counter
static constexpr unsigned int blocklevel = Imp::BlockTraits<B>::blockLevel()+1;
typedef ST size_type;
//! make iterators available as types
typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
typedef typename compressed_block_vector_unmanaged<B,ST>::Iterator Iterator;
//! make iterators available as types
typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
typedef typename compressed_block_vector_unmanaged<B,ST>::ConstIterator ConstIterator;
//===== constructors and such
//! makes empty array
CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,ST>()
{ }
//! make array from given pointers and size
......@@ -1213,7 +1053,7 @@ namespace Imp {
//! assign from scalar
CompressedBlockVectorWindow& operator= (const field_type& k)
{
(static_cast<compressed_block_vector_unmanaged<B,A>&>(*this)) = k;
(static_cast<compressed_block_vector_unmanaged<B,ST>&>(*this)) = k;
return *this;
}
......
dune/istl/cholmod.hh
View file @ bc1f69fc
This diff is collapsed.
dune/istl/common/CMakeLists.txt 0 → 100644
View file @ bc1f69fc
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#install headers
install(FILES
counter.hh
registry.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/istl/common)
dune/istl/common/counter.hh 0 → 100644
View file @ bc1f69fc
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_ISTL_COMMON_COUNTER_HH
#define DUNE_ISTL_COMMON_COUNTER_HH
#include <cassert>
#include <typeinfo>
#include <iostream>
#include <memory>
#include <tuple>
#include <utility>
#include <dune/common/typeutilities.hh>
constexpr std::size_t maxcount = 100;
#define DUNE_GET_COUNTER(Tag) \
(counterFunc(Dune::PriorityTag<maxcount>{}, Tag{}, Dune::CounterImpl::ADLTag{}))
#define DUNE_INC_COUNTER(Tag) \
namespace { \
namespace CounterImpl { \
constexpr std::size_t \
counterFunc(Dune::PriorityTag<DUNE_GET_COUNTER(Tag)+1> p, Tag, ADLTag) \
{ \
return p.value; \
} \
} \
} \
static_assert(true, "unfudge indentation")
namespace Dune {
namespace {
namespace CounterImpl {
struct ADLTag {};
template<class Tag>
constexpr std::size_t counterFunc(Dune::PriorityTag<0>, Tag, ADLTag)
{
return 0;
}
} // end namespace CounterImpl
} // end empty namespace
} // end namespace Dune
#endif // DUNE_ISTL_COMMON_COUNTER_HH