Commit f6447413 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

removed a lot obsolete code

parent 9e2dadcb
#pragma once
namespace AMDiS
{
namespace Impl
{
// Type-Trait needs to be specialized by linear algebra libraries
template <class Vector>
struct BaseVector { using type = Vector; };
template <class Matrix>
struct BaseMatrix { using type = Matrix; };
}
template <class Vector>
using BaseVector_t = typename Impl::BaseVector<Vector>::type;
template <class Matrix>
using BaseMatrix_t = typename Impl::BaseMatrix<Matrix>::type;
template <class T>
struct Triplet
{
std::size_t row;
std::size_t col;
T value;
};
namespace Impl
{
template <class Vector>
class BaseWrapper
{
public:
using value_type = typename Vector::value_type;
using size_type = typename Vector::size_type;
explicit BaseWrapper(Vector& vec) : vec(vec) {}
void resize(size_type s)
{
vec.resize(s);
}
size_type size() const
{
return vec.size();
}
value_type& operator[](size_type i) { return vec[i]; }
value_type const& operator[](size_type i) const { return vec[i]; }
protected:
Vector& vec;
};
} // end namespace Impl
template <class Vector>
Impl::BaseWrapper<std::remove_reference_t<Vector>> wrapper(Vector&& vec)
{
return {std::forward<Vector>(vec)};
}
} // end namespace AMDiS
......@@ -12,7 +12,6 @@
*/
#include <amdis/common/Utility.hpp>
#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
#include <amdis/linear_algebra/PreconditionerInterface.hpp>
#include <amdis/linear_algebra/RunnerInterface.hpp>
......@@ -25,12 +24,11 @@ namespace AMDiS
class LinearSolverInterface
{
protected:
using RunnerBase = RunnerInterface<Matrix, BaseVector_t<VectorX>,
BaseVector_t<VectorB>>;
using RunnerBase = RunnerInterface<Matrix, VectorX, VectorB>;
public:
/// Destructor.
virtual ~LinearSolverInterface() {}
virtual ~LinearSolverInterface() = default;
/// Public method to call in order to solve a linear system Ax = b.
/**
......
#pragma once
#include <dune/common/ftraits.hh>
#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
namespace Dune
{
template <class Elt, class Parameters>
struct FieldTraits<mtl::compressed2D<Elt,Parameters>>
{
using field_type = typename FieldTraits<Elt>::field_type;
using real_type = typename FieldTraits<Elt>::real_type;
};
template <class Value, class Parameters>
struct FieldTraits<mtl::dense2D<Value,Parameters>>
{
using field_type = typename FieldTraits<Value>::field_type;
using real_type = typename FieldTraits<Value>::real_type;
};
template <class Value, class Parameters>
struct FieldTraits<mtl::dense_vector<Value,Parameters>>
{
using field_type = typename FieldTraits<Value>::field_type;
using real_type = typename FieldTraits<Value>::real_type;
};
} // end namespace Dune
......@@ -7,7 +7,7 @@ namespace AMDiS
struct PreconditionerInterface
{
/// Virtual destructor.
virtual ~PreconditionerInterface() {}
virtual ~PreconditionerInterface() = default;
/// Is called a the beginning of a solution procedure
virtual void init(Matrix const& fullMatrix) = 0;
......
......@@ -14,7 +14,7 @@ namespace AMDiS
public:
/// virtual destructor
virtual ~RunnerInterface() {}
virtual ~RunnerInterface() = default;
/// Is called at the beginning of a solution procedure
virtual void init(Matrix const& A) = 0;
......
#pragma once
#include "itl/block_diagonal.hpp"
#include <amdis/CreatorMap.hpp>
#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
namespace AMDiS
{
template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
class DefaultCreators<PreconditionerInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>>
{
using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
using PreconBase = PreconditionerInterface<Matrix, Vector>;
template <template<class> class ITLPrecon>
using PreconCreator
= typename Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>::Creator;
using Map = CreatorMap<PreconBase>;
public:
static void init()
{
auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
Map::addCreator("diag", pc_diag);
Map::addCreator("jacobi", pc_diag);
auto pc_id = new PreconCreator<IdentityPreconditioner>;
Map::addCreator("identity", pc_id);
Map::addCreator("no", pc_id);
Map::addCreator("default", pc_id);
}
};
} // end namespace AMDiS
#pragma once
#include <amdis/CreatorMap.hpp>
#include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
namespace AMDiS
{
/// Adds default creators for linear solvers based on `BlockMTLMatrix`.
/**
* Adds creators for block matrix aware solvers.
* - *cg*: conjugate gradient method, \see cg_solver_type
* - *cgs*: stabilized conjugate gradient mtheod, \see cgs_solver_type
* - *bcgs*: stabilized bi-conjugate gradient method, \see bicgstab_type
* - *tfqmr*: Transposed-Free Quasi-Minimal Residual method, \see tfqmr_solver_type
* - *bcgsl*: stabilized BiCG(ell) method, \see bicgstab_ell_type
* - *gmres*: Generalized minimal residula method, \see gmres_type
* - *fgmres*: Flexible GMRES, \see fgmres_type
* - *minres*: Minimal residul method, \see minres_solver_type
* - *gcr*: Generalized conjugate residual method, \see gcr_type
* - *umfpack*: external UMFPACK solver, \see UmfpackRunner
**/
template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
class DefaultCreators<LinearSolverInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector, Vector>>
{
using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
template <class ITLSolver>
using SolverCreator
= typename LinearSolver<Matrix, Vector,
KrylovRunner<ITLSolver, Matrix, BaseVector_t<Vector>>>::Creator;
#ifdef HAVE_UMFPACK
using UmfpackSolverCreator
= typename LinearSolver<Matrix, Vector,
UmfpackRunner<Matrix, BaseVector_t<Vector>>>::Creator;
#endif
using Map = CreatorMap<SolverBase>;
public:
static void init()
{
auto cg = new SolverCreator<cg_solver_type>;
Map::addCreator("cg", cg);
auto cgs = new SolverCreator<cgs_solver_type>;
Map::addCreator("cgs", cgs);
auto bcgs = new SolverCreator<bicgstab_type>;
Map::addCreator("bicgstab", bcgs);
Map::addCreator("bcgs", bcgs);
auto tfqmr = new SolverCreator<tfqmr_solver_type>;
Map::addCreator("tfqmr", tfqmr);
auto bcgsl = new SolverCreator<bicgstab_ell_type>;
Map::addCreator("bicgstab_ell", bcgsl);
Map::addCreator("bcgsl", bcgsl);
auto gmres = new SolverCreator<gmres_type>;
Map::addCreator("gmres", gmres);
auto fgmres = new SolverCreator<fgmres_type>;
Map::addCreator("fgmres", fgmres);
auto minres = new SolverCreator<minres_solver_type>;
Map::addCreator("minres", minres);
auto gcr = new SolverCreator<gcr_type>;
Map::addCreator("gcr", gcr);
#ifdef HAVE_UMFPACK
auto umfpack = new UmfpackSolverCreator;
Map::addCreator("umfpack", umfpack);
Map::addCreator("direct", umfpack);
#endif
Map::addCreator("default", gmres);
}
};
} // end namespace AMDiS
/** \file BlockMTLMatrix.hpp */
#pragma once
#include <array>
#include <boost/numeric/mtl/matrices.hpp>
#include <amdis/common/Utility.hpp>
#include <amdis/common/Literals.hpp>
#include <amdis/common/Loops.hpp>
#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
namespace AMDiS
{
/// A wrapper for AMDiS::SolverMatrix to be used in MTL/ITL solvers
template <class MTLMatrix, std::size_t _N, std::size_t _M>
class BlockMTLMatrix
: public std::array<std::array<MTLMatrix, _M>, _N>
{
using Self = BlockMTLMatrix;
public:
/// The index/size - type
using size_type = typename MTLMatrix::size_type;
/// The type of the elements of the MTLMatrix
using value_type = typename MTLMatrix::value_type;
/// The underlying mtl matrix type
using BaseMatrix = MTLMatrix;
public:
/// Return the (R,C)'th matrix block
template <std::size_t R, std::size_t C>
auto& operator()(const index_t<R>, const index_t<C>)
{
static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
return std::get<C>(std::get<R>(*this));
}
/// Return the (R,C)'th matrix block
template <std::size_t R, std::size_t C>
auto const& operator()(const index_t<R>, const index_t<C>) const
{
static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
return std::get<C>(std::get<R>(*this));
}
/// Return the number of row blocks
static constexpr std::size_t N() { return _N; }
/// Return the number of column blocks
static constexpr std::size_t M() { return _M; }
/// perform blockwise multiplication A*b -> x
template <class VectorIn, class VectorOut, class Assign>
void mult(VectorIn const& b, VectorOut& x, Assign) const
{
// create iranges to access array blocks
std::array<mtl::irange, _N> r_rows;
std::array<mtl::irange, _M> r_cols;
getRanges(r_rows, r_cols);
forEach(range_<0, _N>, [&](const auto _i)
{
bool first = true;
// a reference to the i'th block of x
VectorOut x_i(x[r_rows[_i]]);
forEach(range_<0, _M>, [&](const auto _j)
{
auto const& A_ij = this->operator()(_i, _j);
if (num_rows(A_ij) > 0 && A_ij.nnz() > 0) {
// a reference to the j'th block of b
const VectorIn b_j(b[r_cols[_j]]);
if (first) {
Assign::first_update(x_i, A_ij * b_j);
first = false;
}
else {
Assign::update(x_i, A_ij * b_j);
}
}
});
});
}
/// A Multiplication operator returns a multiplication-expresssion.
/// Calls \ref mult internally.
template <class VectorIn>
mtl::vec::mat_cvec_multiplier<Self, VectorIn>
operator*(VectorIn const& v) const
{
return {*this, v};
}
/// Fill an array of irange corresponding to the row-sizes, used
/// to access sub-vectors
void getRowRanges(std::array<mtl::irange, _N>& r_rows) const
{
std::size_t start = 0;
forEach(range_<0, _N>, [&](const auto _r) {
std::size_t finish = start + num_rows((*this)(_r, 0_c));
r_rows[_r].set(start, finish);
start = finish;
});
}
/// Fill an array of irange corresponding to the column-sizes, used
/// to access sub-vectors
void getColRanges(std::array<mtl::irange, _M>& r_cols) const
{
std::size_t start = 0;
forEach(range_<0, _M>, [&](const auto _c) {
std::size_t finish = start + num_cols((*this)(0_c, _c));
r_cols[_c].set(start, finish);
start = finish;
});
}
/// Fill two arrays of irange corresponding to row and column sizes.
/// \see getRowRanges() and \see getColRanges()
void getRanges(std::array<mtl::irange, _N>& r_rows,
std::array<mtl::irange, _M>& r_cols) const
{
getRowRanges(r_rows);
getColRanges(r_cols);
}
};
namespace Impl
{
/// Specialization of Impl::MTLMatrix from \file LinearAlgebraBase.hpp
template <class MTLMatrix, std::size_t _N, std::size_t _M>
struct BaseMatrix<BlockMTLMatrix<MTLMatrix, _N, _M>>
{
using type = MTLMatrix;
};
}
/// Return the number of overall rows of a BlockMTLMatrix
template <class MTLMatrix, std::size_t _N, std::size_t _M>
inline std::size_t num_rows(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
{
std::size_t nRows = 0;
forEach(range_<0, _N>, [&](const auto _r) {
nRows += num_rows(A(_r, 0_c));
});
return nRows;
}
/// Return the number of overall columns of a BlockMTLMatrix
template <class MTLMatrix, std::size_t _N, std::size_t _M>
inline std::size_t num_cols(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
{
std::size_t nCols = 0;
forEach(range_<0, _M>, [&](const auto _c) {
nCols += num_cols(A(0_c, _c));
});
return nCols;
}
/// Return the size, i.e. rows*columns of a BlockMTLMatrix
template <class MTLMatrix, std::size_t _N, std::size_t _M>
inline std::size_t size(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
{
return num_rows(A) * num_cols(A);
}
/// Nullify a BlockMTLMatrix, i.e. nullify each block.
template <class MTLMatrix, std::size_t _N, std::size_t _M>
inline void set_to_zero(BlockMTLMatrix<MTLMatrix, _N, _M>& A)
{
forEach(range_<0, _N>, [&](const auto _r) {
forEach(range_<0, _M>, [&](const auto _c) {
set_to_zero(A(_r,_c));
});
});
}
} // end namespace AMDiS
/// \cond HIDDEN_SYMBOLS
namespace mtl
{
template <class MTLMatrix, std::size_t _N, std::size_t _M>
struct Collection<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
{
using value_type = typename MTLMatrix::value_type;
using size_type = typename MTLMatrix::size_type;
};
namespace ashape
{
template <class MTLMatrix, std::size_t _N, std::size_t _M>
struct ashape_aux<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
{
using type = nonscal;
};
} // end namespace ashape
} // end namespace mtl
/// \endcond
#pragma once
#include <array>
#include <type_traits>
#include <boost/numeric/mtl/utility/irange.hpp>
#include <amdis/common/Mpl.hpp>
#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
namespace AMDiS
{
/// A simple block-vector class
template <class MTLVector, std::size_t N>
class BlockMTLVector
: public std::array<MTLVector, N>
{
public:
/// The type of the blocks
using BaseVector = MTLVector;
/// The index/size - type
using size_type = typename MTLVector::size_type;
/// The type of the elements of the MTLVector
using value_type = typename MTLVector::value_type;
};
namespace Impl
{
/// Specialization of Impl::BaseVector from \file LinearAlgebraBase.hpp
template <class MTLVector, std::size_t N>
struct BaseVector<BlockMTLVector<MTLVector, N>>
{
using type = MTLVector;
};
}
/// Return the number of overall rows of a BlockMTLVector
template <class MTLVector, std::size_t N>
inline std::size_t num_rows(BlockMTLVector<MTLVector, N> const& vec)
{
std::size_t nRows = 0;
forEach(range_<0, N>, [&](const auto _i) {
nRows += num_rows(std::get<_i>(vec));
});
return nRows;
}
/// Return the number of overall columns of a BlockMTLVector
template <class MTLVector, std::size_t N>
inline std::size_t num_cols(BlockMTLVector<MTLVector, N> const& vec)
{
return 1;
}
/// Return the size, i.e. rows*columns of a BlockMTLVector
template <class MTLVector, std::size_t N>
inline std::size_t size(BlockMTLVector<MTLVector, N> const& vec)
{
return num_rows(vec);
}
/// Nullify a BlockMTLVector, i.e. nullify each block.
template <class MTLVector, std::size_t N>
inline void set_to_zero(BlockMTLVector<MTLVector, N>& vec)
{
forEach(range_<0, N>, [&](const auto _i) {
set_to_zero(std::get<_i>(vec));
});
}
/// A wrapper, that creates a contiguos vector corresponding to a block-vector
/// and copy the value on construction and eventually back on destruction, if
/// required.
template <class BlockVector, class Vector, class NonConstBlockVector>
class BlockVectorWrapper;
// specialization for BlockMTLVector
template <class BlockVector, class Vector, class MTLVector, std::size_t N>
class BlockVectorWrapper<BlockVector, Vector, BlockMTLVector<MTLVector, N>>
{
static_assert( std::is_same< std::remove_const_t<BlockVector>, BlockMTLVector<MTLVector, N> >::value,
"This specialization is for BlockMTLVectors only.");
public:
explicit BlockVectorWrapper(BlockVector& blockVector,
bool copyBack = !std::is_const<BlockVector>::value)
: blockVector(blockVector)
, vector_(num_rows(blockVector))
, copyBack(copyBack)
{
assignTo();
}
~BlockVectorWrapper()
{
if (copyBack)
assignFrom(bool_<!std::is_const<BlockVector>::value>);
}
/// Return a reference to the block-vector
BlockVector const& getBlockVector() const { return blockVector; }
/// Return a reference to the contiguose-vector
Vector& vector() { return vector_; }
Vector const& vector() const { return vector_; }
private:
/// copy from block-vector to vector
void assignTo()
{
std::size_t start = 0;
for (std::size_t r = 0; r < N; ++r) {
std::size_t finish = start + num_rows(blockVector[r]);
mtl::irange range(start, finish);
vector_[range] = blockVector[r];