Commit ff4baac0 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files
parents f217d1eb 2fd03136
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_CG_INCLUDE
......@@ -16,7 +16,6 @@
#include <cmath>
#include <cassert>
#include <iostream>
#include <boost/mpl/bool.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/itl/itl_fwd.hpp>
......@@ -29,15 +28,14 @@
#include <boost/numeric/mtl/operation/conj.hpp>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/operation/lazy.hpp>
#include <boost/numeric/mtl/utility/static_assert.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
namespace itl {
/// Conjugate Gradients without preconditioning
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Iteration& iter)
{
mtl::vampir_trace<7001> tracer;
......@@ -48,20 +46,20 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Scalar rho(0), rho_1(0), alpha(0), alpha_1(0);
Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
r = b - A*x;
rho = dot(r, r);
while (! iter.finished(Real(sqrt(abs(rho))))) {
++iter;
if (iter.first())
p = r;
else
p = r + (rho / rho_1) * p;
else
p = r + (rho / rho_1) * p;
// q = A * p; alpha = rho / dot(p, q);
(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
alpha= rho / alpha_1;
x += alpha * p;
rho_1 = rho;
(lazy(r) -= alpha * q) || (lazy(rho) = lazy_unary_dot(r));
......@@ -71,9 +69,9 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
}
/// Conjugate Gradients
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Preconditioner, typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& L, Iteration& iter)
{
using pc::is_identity;
......@@ -88,7 +86,7 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
Scalar rho(0), rho_1(0), rr, alpha(0), alpha_1;
Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
r = b - A*x;
rr = dot(r, r);
while (! iter.finished(Real(sqrt(abs(rr))))) {
......@@ -97,12 +95,12 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
if (iter.first())
p = z;
else
else
p = z + (rho / rho_1) * p;
(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
alpha= rho / alpha_1;
x += alpha * p;
rho_1 = rho;
(lazy(r) -= alpha * q) || (lazy(rr) = lazy_unary_dot(r));
......@@ -111,9 +109,9 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
}
/// Conjugate Gradients with ignored right preconditioner to unify interface
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
typename Preconditioner, typename RightPreconditioner, typename Iteration >
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
const Preconditioner& L, const RightPreconditioner&, Iteration& iter)
{
return cg(A, x, b, L, iter);
......@@ -121,7 +119,7 @@ int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
/// Solver class for CG method; right preconditioner ignored (prints warning if not identity)
/** Methods inherited from \ref base_solver. **/
template < typename LinearOperator, typename Preconditioner,
template < typename LinearOperator, typename Preconditioner,
typename RightPreconditioner>
class cg_solver
: public base_solver< cg_solver<LinearOperator, Preconditioner, RightPreconditioner>, LinearOperator >
......@@ -129,16 +127,16 @@ class cg_solver
typedef base_solver< cg_solver<LinearOperator, Preconditioner, RightPreconditioner>, LinearOperator > base;
public:
/// Construct solver from a linear operator; generate (left) preconditioner from it
explicit cg_solver(const LinearOperator& A) : base(A), L(A)
explicit cg_solver(const LinearOperator& A) : base(A), L(A)
{
// MTL_STATIC_ASSERT((!pc::static_is_identity<RightPreconditioner>::value),
// static_assert((!pc::static_is_identity<RightPreconditioner>::value),
// "Right preconditioner must be identity!");
if (!pc::static_is_identity<RightPreconditioner>::value)
std::cerr << "Right Preconditioner ignored!" << std::endl;
}
/// Construct solver from a linear operator and (left) preconditioner
cg_solver(const LinearOperator& A, const Preconditioner& L) : base(A), L(L)
cg_solver(const LinearOperator& A, const Preconditioner& L) : base(A), L(L)
{
if (!pc::static_is_identity<RightPreconditioner>::value)
std::cerr << "Right Preconditioner ignored!" << std::endl;
......
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG, www.simunova.com.
// 2010 SimuNova UG, www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also tools/license/license.mtl.txt in the distribution.
#ifndef ITL_PC_SOLVER_INCLUDE
......@@ -15,14 +15,10 @@
namespace itl {
#ifdef MTL_WITH_TEMPLATE_ALIAS
template <typename Matrix, template <typename, typename, typename> class Solver,
template <typename> class Left, template <typename> class Right>
using pc_solver= Solver<Matrix, Left<Matrix>, Right<Matrix> >;
#endif // MTL_WITH_TEMPLATE_ALIAS
} // namespace itl
#endif // ITL_PC_SOLVER_INCLUDE
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG, www.simunova.com.
// 2010 SimuNova UG, www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also tools/license/license.mtl.txt in the distribution.
#ifndef ITL_REPEATING_SOLVER_INCLUDE
#define ITL_REPEATING_SOLVER_INCLUDE
#include <boost/mpl/if.hpp>
#include <boost/static_assert.hpp>
#include <type_traits>
#include <boost/numeric/itl/itl_fwd.hpp>
namespace itl {
......@@ -25,14 +25,14 @@ namespace itl {
template <typename Solver, unsigned N, bool Stored>
class repeating_solver
{
typedef typename boost::mpl::if_c<Stored, Solver, const Solver&>::type solver_type;
typedef typename std::conditional<Stored, Solver, const Solver&>::type solver_type;
public:
explicit repeating_solver(const Solver& s) : s(s) {}
template <typename Matrix>
explicit repeating_solver(const Matrix& A) : s(A)
{
BOOST_STATIC_ASSERT((Stored)); // if matrix is passed class must own solver
static_assert(Stored, ""); // if matrix is passed class must own solver
}
/// Perform N iterations of the referred solver
......
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG, www.simunova.com.
// 2010 SimuNova UG, www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also tools/license/license.mtl.txt in the distribution.
#ifndef ITL_PC_CONCAT_INCLUDE
#define ITL_PC_CONCAT_INCLUDE
#include <boost/mpl/if.hpp>
#include <type_traits>
#include <boost/numeric/mtl/operation/resource.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
#include <boost/numeric/itl/pc/solver.hpp>
......@@ -24,14 +25,14 @@ namespace itl { namespace pc {
template <typename PC1, typename PC2, typename Matrix, bool Store1= true, bool Store2= true>
class concat
{
typedef typename boost::mpl::if_c<Store1, PC1, const PC1&>::type pc1_type;
typedef typename boost::mpl::if_c<Store2, PC2, const PC2&>::type pc2_type;
typedef typename std::conditional<Store1, PC1, const PC1&>::type pc1_type;
typedef typename std::conditional<Store2, PC2, const PC2&>::type pc2_type;
public:
/// Construct both preconditioners from matrix \p A
explicit concat(const Matrix& A) : A(A), pc1(A), pc2(A)
{
BOOST_STATIC_ASSERT((Store1 && Store2));
static_assert(Store1 && Store2,"");
}
/// Both preconditioners are already constructed and passed as arguments
......
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_IC_0_INCLUDE
#define ITL_PC_IC_0_INCLUDE
#include <boost/mpl/bool.hpp>
#include <type_traits>
#include <boost/numeric/linear_algebra/identity.hpp>
#include <boost/numeric/linear_algebra/inverse.hpp>
......@@ -99,7 +99,7 @@ class ic_0
void adjoint_solve(const VectorIn& x, VectorOut& y) const
{
mtl::vampir_trace<5044> tracer;
solve(x, y);
solve(x, y);
}
......@@ -109,46 +109,46 @@ class ic_0
protected:
template <typename VectorOut, typename Solver> friend struct ic_0_evaluator;
// Dummy type to perform factorization in initializer list not in
// Dummy type to perform factorization in initializer list not in
struct factorizer
{
factorizer(const Matrix &A, U_type& U)
{ factorize(A, U, mtl::traits::is_sparse<Matrix>(), boost::is_same<Value, typename mtl::Collection<Matrix>::value_type>()); }
{ factorize(A, U, mtl::traits::is_sparse<Matrix>(), std::is_same<Value, typename mtl::Collection<Matrix>::value_type>()); }
template <typename T>
void factorize(const Matrix&, U_type&, boost::mpl::false_, T)
void factorize(const Matrix&, U_type&, std::false_type, T)
{ MTL_THROW_IF(true, mtl::logic_error("IC(0) is not suited for dense matrices")); }
// When we change the value_type then the factorization is still performed with that of A
template <typename UF>
void factorize(const Matrix& A, UF& U, boost::mpl::true_, boost::mpl::false_)
void factorize(const Matrix& A, UF& U, std::true_type, std::false_type)
{
typedef mtl::mat::compressed2D<typename mtl::Collection<Matrix>::value_type, para> tmp_type;
tmp_type U_tmp;
factorize(A, U_tmp, boost::mpl::true_(), boost::mpl::true_());
factorize(A, U_tmp, std::true_type(), std::true_type());
U= U_tmp;
}
// Factorization adapted from Saad
// Undefined (runtime) behavior if matrix is not symmetric
// Undefined (runtime) behavior if matrix is not symmetric
// UF is type for the factorization
template <typename UF>
void factorize(const Matrix& A, UF& U, boost::mpl::true_, boost::mpl::true_)
template <typename UF>
void factorize(const Matrix& A, UF& U, std::true_type, std::true_type)
{
using namespace mtl; using namespace mtl::tag; using mtl::traits::range_generator;
using namespace mtl; using namespace mtl::tag; using mtl::traits::range_generator;
using math::reciprocal; using mtl::mat::upper;
mtl::vampir_trace<5035> tracer;
// For the factorization we take still the value_type of A and later we copy it maybe to another value_type
typedef typename mtl::Collection<Matrix>::value_type value_type;
typedef typename range_generator<row, UF>::type cur_type;
typedef typename range_generator<nz, cur_type>::type icur_type;
typedef typename range_generator<row, UF>::type cur_type;
typedef typename range_generator<nz, cur_type>::type icur_type;
MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
U= upper(A);
typename mtl::traits::col<UF>::type col(U);
typename mtl::traits::value<UF>::type value(U);
typename mtl::traits::value<UF>::type value(U);
cur_type kc= begin<row>(U), kend= end<row>(U);
for (size_type k= 0; kc != kend; ++kc, ++k) {
......@@ -159,7 +159,7 @@ class ic_0
// U[k][k]= 1.0 / sqrt(U[k][k]);
value_type inv_dia= reciprocal(sqrt(value(*ic)));
value(*ic, inv_dia);
// icur_type jbegin=
// icur_type jbegin=
++ic;
for (; ic != iend; ++ic) {
// U[k][i] *= U[k][k]
......@@ -174,8 +174,8 @@ class ic_0
icur_type jc= begin<nz>(irow), jend= end<nz>(irow);
while (col(*jc) <= k) ++jc;
while (col(*--jend) > i) ;
++jend;
++jend;
for (; jc != jend; ++jc) {
size_type j= col(*jc);
U.lvalue(j, i)-= d * U[k][j];
......@@ -191,7 +191,7 @@ class ic_0
L_type L;
lower_solver_t lower_solver;
upper_solver_t upper_solver;
};
};
#if 0
template <typename Matrix, typename Value, typename Vector>
......@@ -204,9 +204,9 @@ struct ic_0_solver
template <typename VectorOut>
void assign_to(VectorOut& y) const
{ P.solve(x, y); }
{ P.solve(x, y); }
const ic_0<Matrix, Value>& P;
const ic_0<Matrix, Value>& P;
const Vector& x;
};
#endif
......@@ -219,7 +219,7 @@ struct ic_0_evaluator
typedef typename mtl::Collection<VectorOut>::value_type out_value_type;
ic_0_evaluator(VectorOut& y, const Solver& s)
ic_0_evaluator(VectorOut& y, const Solver& s)
: y(y), s(s), U(s.P.U), y0(s.P.solve_lower(s.x, y)) { MTL_DEBUG_ARG(lr= 99999999); }
......
// Software License for MTL
//
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_ILU_0_INCLUDE
#define ITL_PC_ILU_0_INCLUDE
#include <boost/mpl/bool.hpp>
#include <type_traits>
#include <boost/numeric/linear_algebra/identity.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
......@@ -30,7 +31,7 @@
namespace itl { namespace pc {
// Dummy type to perform factorization in initializer list not in
// Dummy type to perform factorization in initializer list not in
struct ilu_0_factorizer
{
template <typename Matrix, typename L_type, typename U_type>
......@@ -38,14 +39,14 @@ struct ilu_0_factorizer
{ factorize(A, L, U, mtl::traits::is_sparse<Matrix>()); }
template <typename Matrix, typename L_type, typename U_type>
void factorize(const Matrix&, L_type&, U_type&, boost::mpl::false_)
void factorize(const Matrix&, L_type&, U_type&, std::false_type)
{ MTL_THROW_IF(true, mtl::logic_error("ILU is not intended for dense matrices")); }
template <typename Matrix, typename L_type, typename U_type>
void factorize(const Matrix& A, L_type& L, U_type& U, boost::mpl::true_)
void factorize(const Matrix& A, L_type& L, U_type& U, std::true_type)
{
using namespace mtl; using namespace mtl::tag; using mtl::traits::range_generator;
using math::reciprocal;
using namespace mtl; using namespace mtl::tag; using mtl::traits::range_generator;
using math::reciprocal;
MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
mtl::vampir_trace<5038> tracer;
......@@ -55,10 +56,10 @@ struct ilu_0_factorizer
typedef mtl::mat::compressed2D<value_type, para> LU_type;
LU_type LU(A);
typedef typename range_generator<row, LU_type>::type cur_type;
typedef typename range_generator<nz, cur_type>::type icur_type;
typedef typename range_generator<row, LU_type>::type cur_type;
typedef typename range_generator<nz, cur_type>::type icur_type;
typename mtl::traits::col<LU_type>::type col(LU);
typename mtl::traits::value<LU_type>::type value(LU);
typename mtl::traits::value<LU_type>::type value(LU);
mtl::vec::dense_vector<value_type, mtl::vec::parameters<> > inv_dia(num_rows(A));
cur_type ic= begin<row>(LU), iend= end<row>(LU);
for (size_type i= 0; ic != iend; ++ic, ++i) {
......@@ -72,14 +73,14 @@ struct ilu_0_factorizer
for (icur_type jc= kc + 1; jc != kend; ++jc)
value(*jc, value(*jc) - aik * LU[k][col(*jc)]);
// std::cout << "LU after eliminating A[" << i << "][" << k << "] =\n" << LU;
// std::cout << "LU after eliminating A[" << i << "][" << k << "] =\n" << LU;
}
inv_dia[i]= reciprocal(LU[i][i]);
}
invert_diagonal(LU);
invert_diagonal(LU);
L= strict_lower(LU);
U= upper(LU);
}
}
};
template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
......
// Software License for MTL
//
// Copyright (c) 2007 The Trustees of Indiana University.
//
// Copyright (c) 2007 The Trustees of Indiana University.
// 2008 Dresden University of Technology and the Trustees of Indiana University.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
//
//
// This file is part of the Matrix Template Library
//
//
// See also license.mtl.txt in the distribution.
#ifndef ITL_PC_ILUT_INCLUDE
#define ITL_PC_ILUT_INCLUDE
#include <type_traits>
#include <boost/numeric/mtl/vector/sparse_vector.hpp>
#include <boost/numeric/mtl/operation/invert_diagonal.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>
......@@ -29,26 +31,26 @@ struct ilut_factorizer
// column-major matrices are copied first
template <typename Matrix, typename Para, typename L_type, typename U_type, bool B>
void factorize(const Matrix &A, const Para& p, L_type& L, U_type& U, boost::mpl::bool_<B>)
void factorize(const Matrix &A, const Para& p, L_type& L, U_type& U, std::integral_constant<bool, B>)
{
typedef typename mtl::Collection<Matrix>::value_type value_type;
typedef typename mtl::Collection<Matrix>::size_type size_type;
typedef mtl::mat::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;
typedef mtl::mat::compressed2D<value_type, para> LU_type;
LU_type LU(A);
factorize(LU, p, L, U, boost::mpl::true_());
factorize(LU, p, L, U, std::true_type());
}
// According Yousef Saad: ILUT, NLAA, Vol 1(4), 387-402 (1994)
#if 0
template <typename Value, typename MPara, typename Para, typename L_type, typename U_type>
factorize(const mtl::mat::compressed2D<Value, MPara>& A, const Para& p, L_type& L, U_type& U, boost::mpl::true_)
factorize(const mtl::mat::compressed2D<Value, MPara>& A, const Para& p, L_type& L, U_type& U, std::true_type)
#endif
template <typename Matrix, typename Para, typename L_type, typename U_type>
void factorize(const Matrix& A, const Para& p, L_type& L, U_type& U, boost::mpl::true_)
void factorize(const Matrix& A, const Para& p, L_type& L, U_type& U, std::true_type)
{
{
mtl::vampir_trace<5049> tracer;
using std::abs; using mtl::traits::range_generator; using mtl::begin; using mtl::end;
using namespace mtl::tag;
......@@ -56,22 +58,22 @@ struct ilut_factorizer
typedef typename mtl::Collection<Matrix>::value_type value_type;
typedef typename mtl::Collection<Matrix>::size_type size_type;
typedef typename range_generator<row, Matrix>::type cur_type;
typedef typename range_generator<nz, cur_type>::type icur_type;
typedef typename range_generator<row, Matrix>::type cur_type;
typedef typename range_generator<nz, cur_type>::type icur_type;
typename mtl::traits::col<Matrix>::type col(A);
typename mtl::traits::const_value<Matrix>::type value(A);
typename mtl::traits::const_value<Matrix>::type value(A);
size_type n= num_rows(A);
L.change_dim(n, n);
L.change_dim(n, n);
U.change_dim(n, n);
{
mtl::mat::inserter<L_type> L_ins(L, p.first);
mtl::mat::inserter<U_type> U_ins(U, p.first + 1); // plus one for diagonal
mtl::sparse_vector<value_type> vec(n); // corr. row in paper
cur_type ic= begin<row>(A); // , iend= end<row>(A);
for (size_type i= 0; i < n; ++i, ++ic) {
for (icur_type kc= begin<nz>(ic), kend= end<nz>(ic); kc != kend; ++kc) // row= A[i][*]
vec.insert(col(*kc), value(*kc));
// std::cerr << "vec_" << i << " = " << vec << std::endl;
......@@ -96,7 +98,7 @@ struct ilut_factorizer
// std::cerr << "vec_" << i << " = " << vec << std::endl;
vec.sort_on_data();
// std::cerr << "vec_" << i << " sorted on data = " << vec << std::endl;
// std::cout << "vec at " << i << " is " << vec << '\n';
// mtl::vampir_trace<9904> tracer2;
bool diag_found= false;
......@@ -111,7 +113,7 @@ struct ilut_factorizer
U_ins[i][k] << v;
} else // i > k
if (cntl++ < p.first)
L_ins[i][k] << v;
L_ins[i][k] << v;
}
if (!diag_found) std::cerr << "Deleted diagonal!!!!\n";
vec.make_empty();
......@@ -128,7 +130,7 @@ class ilut
{
typedef ilu<Matrix, ilut_factorizer, Value> base;
public:
ilut(const Matrix& A, std::size_t p, typename mtl::Collection<Matrix>::value_type tau)
ilut(const Matrix& A, std::size_t p, typename mtl::Collection<Matrix>::value_type tau)
: base(A, std::make_pair(p, tau)) {}
};
......