diff --git a/cmake/modules/AmdisMacros.cmake b/cmake/modules/AmdisMacros.cmake
index 8ae3c82b0955582c19de4adfd6c56eed5216bd45..b670bfcd309a3eb030325317e7b43ffa1b1c17e8 100644
--- a/cmake/modules/AmdisMacros.cmake
+++ b/cmake/modules/AmdisMacros.cmake
@@ -24,10 +24,14 @@ endif (NOT BACKEND)
if (BACKEND STREQUAL "MTL" OR BACKEND STREQUAL "MTL4")
find_package(MTL REQUIRED)
set(HAVE_MTL TRUE)
- if (MTL_FOUND)
- message(STATUS " Found MTL, version: ${MTL_VERSION}")
- dune_register_package_flags(LIBRARIES MTL::MTL COMPILE_DEFINITIONS "ENABLE_MTL=1")
- endif (MTL_FOUND)
+ message(STATUS " Found MTL, version: ${MTL_VERSION}")
+ dune_register_package_flags(LIBRARIES MTL::MTL COMPILE_DEFINITIONS "ENABLE_MTL=1")
+
+ find_package(HYPRE)
+ if (HYPRE_FOUND)
+ set(HAVE_HYPRE TRUE)
+ dune_register_package_flags(LIBRARIES HYPRE::HYPRE COMPILE_DEFINITIONS "ENABLE_HYPRE=1")
+ endif(HYPRE_FOUND)
elseif (BACKEND STREQUAL "EIGEN" OR BACKEND STREQUAL "EIGEN3")
find_package(Eigen3 REQUIRED 3.3.5)
set(HAVE_EIGEN TRUE)
diff --git a/cmake/modules/CMakeLists.txt b/cmake/modules/CMakeLists.txt
index 6d866c15744a1740cbf0c02832a8ff3495953c7f..5474caf6870ee17a78f4d6204c5cf42c1d8ba0b0 100644
--- a/cmake/modules/CMakeLists.txt
+++ b/cmake/modules/CMakeLists.txt
@@ -2,6 +2,7 @@ install(FILES
AddAmdisExecutable.cmake
AmdisMacros.cmake
AmdisCXXFeatures.cmake
+ FindHYPRE.cmake
FindMTL.cmake
FindPETSc.cmake
PkgConfigLinkLibraries.cmake
diff --git a/cmake/modules/FindHYPRE.cmake b/cmake/modules/FindHYPRE.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..f3fbf34ad47106e00798e3397db8ba41e0add0ed
--- /dev/null
+++ b/cmake/modules/FindHYPRE.cmake
@@ -0,0 +1,75 @@
+# FindHYPRE.cmake
+#
+# Finds the HYPRE library
+#
+# This will define the following variables
+#
+# HYPRE_FOUND
+# HYPRE_INCLUDE_DIR
+# HYPRE_LIBRARIES
+#
+# and the following imported targets
+#
+# HYPRE::HYPRE
+#
+# Author: Simon Praetorius <simon.praetorius@tu-dresden.de>
+
+
+mark_as_advanced(HYPRE_FOUND HYPRE_INCLUDE_DIR HYPRE_LIBRARIES)
+
+find_path(HYPRE_INCLUDE_DIR HYPRE.h
+ PATHS
+ /opt/software/hypre
+ ${HYPRE_DIR} ${HYPRE_ROOT}
+ ENV HYPRE_DIR ENV HYPRE_ROOT
+ ENV EBROOTHYPRE
+ PATH_SUFFIXES
+ hypre include include/hypre
+ NO_DEFAULT_PATH)
+find_path(HYPRE_INCLUDE_DIR HYPRE.h
+ PATH_SUFFIXES hypre include include/hypre)
+
+find_library(HYPRE_LIBRARY HYPRE
+ PATHS
+ /opt/software/hypre
+ ${HYPRE_DIR} ${HYPRE_ROOT}
+ ENV HYPRE_DIR ENV HYPRE_ROOT
+ ENV EBROOTHYPRE
+ PATH_SUFFIXES
+ lib
+ NO_DEFAULT_PATH)
+find_library(HYPRE_LIBRARY HYPRE
+ PATH_SUFFIXES lib)
+
+if (HYPRE_LIBRARY)
+ set(HYPRE_LIBRARIES ${HYPRE_LIBRARY})
+ get_filename_component(HYPRE_LIBRARY_DIR ${HYPRE_LIBRARY} DIRECTORY)
+
+ foreach(_lib_name "HYPRE_FEI" "HYPRE_core")
+ find_library(_lib ${_lib_name} HINTS ${HYPRE_LIBRARY_DIR} NO_DEFAULT_PATH)
+ if (_lib)
+ list(APPEND HYPRE_LIBRARIES ${_lib})
+ endif (_lib)
+ unset(_lib)
+ endforeach()
+ unset(_lib_name)
+endif (HYPRE_LIBRARY)
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(HYPRE
+ REQUIRED_VARS HYPRE_INCLUDE_DIR HYPRE_LIBRARY
+)
+
+# text for feature summary
+set_package_properties("HYPRE" PROPERTIES
+ DESCRIPTION "high performance preconditioners")
+
+
+if(HYPRE_FOUND AND NOT TARGET HYPRE::HYPRE)
+ add_library(HYPRE::HYPRE INTERFACE IMPORTED)
+
+ set_target_properties(HYPRE::HYPRE PROPERTIES
+ INTERFACE_INCLUDE_DIRECTORIES "${HYPRE_INCLUDE_DIR}"
+ INTERFACE_LINK_LIBRARIES "${HYPRE_LIBRARIES}"
+ )
+endif()
diff --git a/cmake/modules/FindMTL.cmake b/cmake/modules/FindMTL.cmake
index 8c20f0752d99a13cfde4541073c4c5619016b4f0..80a84429c1100d6a6d28d04d5a5aebe17bc413f3 100644
--- a/cmake/modules/FindMTL.cmake
+++ b/cmake/modules/FindMTL.cmake
@@ -18,16 +18,15 @@
mark_as_advanced(MTL_FOUND MTL_COMPILE_DEFINITIONS MTL_INCLUDE_DIR)
find_path(MTL_INCLUDE_DIR boost/numeric/mtl/mtl.hpp
- HINTS
- ${CMAKE_CURRENT_SOURCE_DIR}/../../libs/mtl4
- ${CMAKE_CURRENT_SOURCE_DIR}/../../include/mtl4
- ${MTL_DIR}/../../include
- $ENV{MTL_DIR}/../../include
- ENV EBROOTMTL
PATHS
/opt/software/mtl
/opt/sources/mtl
/opt/development/mtl
+ ${MTL_DIR} ${MTL_ROOT}
+ ENV MTL_DIR ENV MTL_ROOT
+ ENV EBROOTMTL
+ PATH_SUFFIXES
+ include
)
set(MTL_COMPILE_DEFINITIONS MTL_ASSERT_FOR_THROW=1)
@@ -40,6 +39,10 @@ find_package_handle_standard_args(MTL
REQUIRED_VARS MTL_INCLUDE_DIR
)
+# text for feature summary
+set_package_properties("MTL" PROPERTIES
+ DESCRIPTION "Matrix Template Library")
+
if(MTL_FOUND AND NOT TARGET MTL::MTL)
add_library(MTL::MTL INTERFACE IMPORTED)
@@ -48,8 +51,12 @@ if(MTL_FOUND AND NOT TARGET MTL::MTL)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
endif (SuiteSparse_FOUND)
+ find_package(HYPRE QUIET)
+ if (HYPRE_FOUND)
+ list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_HYPRE")
+ endif (HYPRE_FOUND)
+
set_target_properties(MTL::MTL PROPERTIES
INTERFACE_INCLUDE_DIRECTORIES "${MTL_INCLUDE_DIR}"
- INTERFACE_COMPILE_DEFINITIONS "${MTL_COMPILE_DEFINITIONS}"
- )
+ INTERFACE_COMPILE_DEFINITIONS "${MTL_COMPILE_DEFINITIONS}")
endif()
diff --git a/config.h.cmake b/config.h.cmake
index ec1f0739d074008e846f16240807478bb9e0a836..9317b1e2b2a8e976ecb3ef4dd1e5123a714febe0 100644
--- a/config.h.cmake
+++ b/config.h.cmake
@@ -43,11 +43,14 @@
/* Define to true if the MTL library is available */
#cmakedefine HAVE_MTL 1
+/* Define to true if the HYPRE library is available */
+#cmakedefine HAVE_HYPRE ENABLE_HYPRE
+
/* Define to true if the Eigen3 library is available */
#cmakedefine HAVE_EIGEN 1
/* Define to true if the PETSc library is available */
-#cmakedefine HAVE_PETSC ENABLE_PETSC
+#cmakedefine HAVE_PETSC 1
/* some detected compiler features may be used in AMDiS */
#cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1
diff --git a/src/amdis/linearalgebra/mtl/CMakeLists.txt b/src/amdis/linearalgebra/mtl/CMakeLists.txt
index 1306877d8b3f4d0fcca5642dc4f6a459a3c68d23..8b79796a4bc8b8463ef0fad70e1a362783b82fd2 100644
--- a/src/amdis/linearalgebra/mtl/CMakeLists.txt
+++ b/src/amdis/linearalgebra/mtl/CMakeLists.txt
@@ -1,10 +1,12 @@
install(FILES
Constraints.hpp
+ HyprePrecon.hpp
ITL_Preconditioner.hpp
ITL_Solver.hpp
KrylovRunner.hpp
MatrixBackend.hpp
Preconditioner.hpp
+ SolverPrecon.hpp
Traits.hpp
UmfpackRunner.hpp
VectorBackend.hpp
diff --git a/src/amdis/linearalgebra/mtl/HyprePrecon.hpp b/src/amdis/linearalgebra/mtl/HyprePrecon.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb6935196b37201b75976a1bf7fd912ed9704067
--- /dev/null
+++ b/src/amdis/linearalgebra/mtl/HyprePrecon.hpp
@@ -0,0 +1,153 @@
+#pragma once
+
+#if HAVE_HYPRE
+
+#include <dune/common/unused.hh>
+#include <amdis/linearalgebra/mtl/itl/hypre.hpp>
+
+namespace AMDiS
+{
+ /// \brief Interface to the HYPRE BoomerAMG preconditioner [...]
+ /**
+ * Parameters provided by AMDiS:
+ *
+ * [PRECON]->cycle mode:
+ * 1...V-cycle
+ * 2...W-cycle
+ *
+ * [PRECON]->interpolation type:
+ * 0...classical modified interpolation
+ * 1...LS interpolation (for use with GSMG)
+ * 2...classical modified interpolation for hyperbolic PDEs
+ * 3...direct interpolation (with separation of weights)
+ * 4...multipass interpolation
+ * 5...multipass interpolation (with separation of weights)
+ * 6...extended+i interpolation
+ * 7...extended+i (if no common C neighbor) interpolation
+ * 8...standard interpolation
+ * 9...standard interpolation (with separation of weights)
+ * 10..classical block interpolation (for use with nodal systems version only)
+ * 11..classical block interpolation (for use with nodal systems version only)
+ * with diagonalized diagonal blocks
+ * 12..FF interpolation
+ * 13..FF1 interpolation
+ * 14..extended interpolation
+ *
+ * [PRECON]->info:
+ * 0...no printout (default)
+ * 1...print setup information
+ * 2...print solve information
+ * 3...print both setup and solve information
+ *
+ * [PRECON]->relax type:
+ * 0...Jacobi
+ * 1...Gauss-Seidel, sequential (very slow!)
+ * 2...Gauss-Seidel, interior points in parallel, boundary sequential (slow!)
+ * 3...hybrid Gauss-Seidel or SOR, forward solve
+ * 4...hybrid Gauss-Seidel or SOR, backward solve
+ * 5...hybrid chaotic Gauss-Seidel (works only with OpenMP)
+ * 6...hybrid symmetric Gauss-Seidel or SSOR
+ * 9...Gaussian elimination (only on coarsest level)
+ **/
+ template <class Traits>
+ class HyprePrecon
+ : public PreconditionerInterface<Traits>
+ {
+ using Self = HyprePrecon;
+ using Matrix = typename Traits::Mat;
+ using Vector = typename Traits::Vec;
+ using Super = PreconditionerInterface<Traits>;
+
+ public:
+ /// A creator to be used instead of the constructor.
+ struct Creator : CreatorInterfaceName<Super>
+ {
+ std::unique_ptr<Super> createWithString(std::string prefix) override
+ {
+ return std::make_unique<Self>(prefix);
+ }
+ };
+
+ public:
+ HyprePrecon(std::string const& prefix)
+ {
+ // read HYPRE parameters
+ int maxIter = 1, cycleMode = -1, interpolation = -1, relaxation = -1, info = 1;
+ Parameters::get(prefix + "->max iteration", maxIter);
+ Parameters::get(prefix + "->cycle mode", cycleMode);
+ Parameters::get(prefix + "->interpolation type", interpolation);
+ Parameters::get(prefix + "->relax type", relaxation);
+ Parameters::get(prefix + "->info", info);
+
+ config_.maxIter(maxIter);
+ config_.interpolationType(interpolation);
+ config_.relaxType(relaxation);
+ config_.cycle(cycleMode);
+ config_.printLevel(info);
+ }
+
+ /// Implementation of \ref PreconditionerBase::init()
+ void init(Matrix const& fullMatrix) override
+ {
+ hypreMatrix_.init(fullMatrix);
+ HYPRE_IJMatrixGetObject(hypreMatrix_, (void**)(&matrix_));
+ HYPRE_BoomerAMGCreate(&solver_);
+
+ // apply solver parameters
+ config_(solver_);
+
+ mtl::dense_vector<double> swap(1, 0.0);
+ mtl::HypreParVector x(swap);
+ HYPRE_BoomerAMGSetup(solver_, matrix_, x, x);
+
+ solverCreated_ = true;
+ }
+
+ /// Implementation of \ref PreconditionerInterface::exit()
+ void exit() override
+ {
+ if (solverCreated_)
+ HYPRE_BoomerAMGDestroy(solver_);
+ solverCreated_ = false;
+ }
+
+ /// Implementation of \ref PreconditionerBase::solve()
+ void solve(Vector const& mtlB, Vector& mtlX) const override
+ {
+ assert(solverCreated_);
+ mtl::HypreParVector x(mtlB); // use b as initial guess
+ mtl::HypreParVector b(mtlB);
+
+ DUNE_UNUSED int error = HYPRE_BoomerAMGSolve(solver_, matrix_, b, x);
+ assert(error != 0);
+
+ // write output back to MTL vector
+ mtl::convert(x.getHypreVector(), mtlX);
+ }
+
+ /// Implementation of \ref PreconditionerBase::adjoint_solve()
+ void adjoint_solve(Vector const& mtlB, Vector& mtlX) const override
+ {
+ assert(solverCreated_);
+ mtl::HypreParVector x(mtlB); // use b as initial guess
+ mtl::HypreParVector b(mtlB);
+
+ DUNE_UNUSED int error = HYPRE_BoomerAMGSolveT(solver_, matrix_, b, x);
+ assert(error != 0);
+
+ // write output back to MTL vector
+ mtl::convert(x.getHypreVector(), mtlX);
+ }
+
+ private:
+ mtl::HypreMatrix hypreMatrix_;
+ HYPRE_ParCSRMatrix matrix_;
+ HYPRE_Solver solver_;
+
+ mtl::AMGConfigurator config_;
+ bool solverCreated_ = false;
+ };
+
+} // end namespace AMDiS
+
+#endif // HAVE_HYPRE
diff --git a/src/amdis/linearalgebra/mtl/ITL_Preconditioner.hpp b/src/amdis/linearalgebra/mtl/ITL_Preconditioner.hpp
index 687d708891e2193350293fa0b189f274616c5f38..1bb2a41312f86be043bb8902deb8cb6387407629 100644
--- a/src/amdis/linearalgebra/mtl/ITL_Preconditioner.hpp
+++ b/src/amdis/linearalgebra/mtl/ITL_Preconditioner.hpp
@@ -6,10 +6,15 @@
#include <boost/numeric/itl/pc/ic_0.hpp>
#include <boost/numeric/mtl/vector/assigner.hpp>
-#include <amdis/linearalgebra/mtl/itl/masslumping.hpp>
+#include <amdis/CreatorMap.hpp>
#include <amdis/linearalgebra/mtl/Preconditioner.hpp>
+#include <amdis/linearalgebra/mtl/SolverPrecon.hpp>
#include <amdis/linearalgebra/mtl/Traits.hpp>
-#include <amdis/CreatorMap.hpp>
+#include <amdis/linearalgebra/mtl/itl/masslumping.hpp>
+
+#if HAVE_HYPRE
+#include <amdis/linearalgebra/mtl/HyprePrecon.hpp>
+#endif
namespace AMDiS
{
@@ -73,6 +78,7 @@ namespace AMDiS
template <class Matrix>
using ICPreconditioner = itl::pc::ic_0<Matrix>;
+
template <class Traits>
class DefaultCreators< PreconditionerInterface<Traits> >
{
@@ -107,20 +113,28 @@ namespace AMDiS
Map::addCreator("identity", pc_id);
Map::addCreator("no", pc_id);
+ auto pc_solver = new typename SolverPrecon<Traits>::Creator;
+ Map::addCreator("solver", pc_solver);
+
+#if HAVE_HYPRE
+ auto pc_hypre = new typename HyprePrecon<Traits>::Creator;
+ Map::addCreator("hypre", pc_hypre);
+#endif
+
Map::addCreator("default", pc_id);
}
};
- template <class Traits>
- itl::pc::solver<PreconditionerInterface<Traits>, typename Traits::Vec, false>
- solve(PreconditionerInterface<Traits> const& P, typename Traits::Vec const& vin)
+ template <class Traits, class Vector>
+ itl::pc::solver<PreconditionerInterface<Traits>, Vector, false>
+ solve(PreconditionerInterface<Traits> const& P, Vector const& vin)
{
return {P, vin};
}
- template <class Traits>
- itl::pc::solver<PreconditionerInterface<Traits>, typename Traits::Vec, true>
- adjoint_solve(PreconditionerInterface<Traits> const& P, typename Traits::Vec const& vin)
+ template <class Traits, class Vector>
+ itl::pc::solver<PreconditionerInterface<Traits>, Vector, true>
+ adjoint_solve(PreconditionerInterface<Traits> const& P, Vector const& vin)
{
return {P, vin};
}
diff --git a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
index bc2360a6f047fda0d0b05bc5cba4a41d86fcc544..3aeae41e922182c8b8895e3acf1a7dba3ef82ec2 100644
--- a/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
+++ b/src/amdis/linearalgebra/mtl/ITL_Solver.hpp
@@ -36,10 +36,10 @@ namespace AMDiS
{
public:
explicit cg_solver_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::cg(A, x, b, l, r, iter);
+ return itl::cg(A, x, b, precon, iter);
}
};
@@ -56,10 +56,10 @@ namespace AMDiS
{
public:
explicit cgs_solver_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::cgs(A, x, b, l, iter);
+ return itl::cgs(A, x, b, precon, iter);
}
};
@@ -76,10 +76,10 @@ namespace AMDiS
{
public:
explicit bicg_solver_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::bicg(A, x, b, l, iter);
+ return itl::bicg(A, x, b, precon, iter);
}
};
@@ -96,10 +96,10 @@ namespace AMDiS
{
public:
explicit bicgstab_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::bicgstab(A, x, b, l, iter);
+ return itl::bicgstab(A, x, b, precon, iter);
}
};
@@ -116,10 +116,10 @@ namespace AMDiS
{
public:
explicit bicgstab2_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::bicgstab_2(A, x, b, l, iter);
+ return itl::bicgstab_2(A, x, b, precon, iter);
}
};
@@ -135,10 +135,11 @@ namespace AMDiS
{
public:
explicit qmr_solver_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::qmr(A, x, b, l, r, iter);
+ itl::pc::identity<LinOp> id(A);
+ return itl::qmr(A, x, b, precon, id, iter);
}
};
@@ -149,16 +150,17 @@ namespace AMDiS
* Quasi-Minimal Residual method \implements ITL_Solver
*
* Solves a linear system by the Transposed-Free Quasi-Minimal Residual method
- * (TFQMR). Does not use preconditioning currently.
+ * (TFQMR).
*/
class tfqmr_solver_type
{
public:
explicit tfqmr_solver_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::tfqmr(A, x, b, l, r, iter);
+ itl::pc::identity<LinOp> id(A);
+ return itl::tfqmr(A, x, b, precon, id, iter);
}
};
@@ -179,10 +181,11 @@ namespace AMDiS
{
Parameters::get(name + "->ell", ell);
}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::bicgstab_ell(A, x, b, l, r, iter, ell);
+ itl::pc::identity<LinOp> id(A);
+ return itl::bicgstab_ell(A, x, b, precon, id, iter, ell);
}
};
@@ -205,17 +208,18 @@ namespace AMDiS
Parameters::get(name + "->restart", restart);
Parameters::get(name + "->orthogonalization", ortho);
}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
+ itl::pc::identity<LinOp> id(A);
switch (ortho)
{
default:
case GRAM_SCHMIDT:
- return itl::gmres2(A, x, b, l, r, iter, restart);
+ return itl::gmres2(A, x, b, id, precon, iter, restart);
break;
case HOUSEHOLDER:
- return itl::gmres_householder(A, x, b, l, iter, restart);
+ return itl::gmres_householder(A, x, b, precon, iter, restart);
break;
}
}
@@ -240,6 +244,8 @@ namespace AMDiS
* Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast
* algorithms for solving large nonsymmetric linear systems.
* SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
+ *
+ * NOTE: preconditioners are ignored.
*/
class idr_s_type
{
@@ -249,10 +255,11 @@ namespace AMDiS
{
Parameters::get(name + "->s", s);
}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const&, I& iter)
{
- return itl::idr_s(A, x, b, l, r, iter, s);
+ itl::pc::identity<LinOp> id(A);
+ return itl::idr_s(A, x, b, id, id, iter, s);
}
};
@@ -269,10 +276,10 @@ namespace AMDiS
{
public:
explicit minres_solver_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::minres(A, x, b, l, r, iter);
+ return itl::minres(A, x, b, precon, iter);
}
};
@@ -295,10 +302,11 @@ namespace AMDiS
{
Parameters::get(name + "->restart", restart);
}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::gcr(A, x, b, l, r, iter, restart);
+ itl::pc::identity<LinOp> id(A);
+ return itl::gcr(A, x, b, id, precon, iter, restart);
}
};
@@ -322,17 +330,18 @@ namespace AMDiS
Parameters::get(name + "->restart", restart);
Parameters::get(name + "->orthogonalization", orthogonalization);
}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
+ itl::pc::identity<LinOp> id(A);
switch (orthogonalization)
{
default:
case GRAM_SCHMIDT:
- return itl::fgmres(A, x, b, l, r, iter, restart);
+ return itl::fgmres(A, x, b, id, precon, iter, restart);
break;
case HOUSEHOLDER:
- return itl::fgmres_householder(A, x, b, r, iter, restart);
+ return itl::fgmres_householder(A, x, b, precon, iter, restart);
break;
}
}
@@ -357,10 +366,10 @@ namespace AMDiS
{
public:
explicit preonly_type(std::string const& /*name*/) {}
- template <class LinOp, class X, class B, class L, class R, class I>
- int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& /*r*/, I& iter)
+ template <class LinOp, class X, class B, class P, class I>
+ int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
- return itl::preonly(A, x, b, l, iter);
+ return itl::preonly(A, x, b, precon, iter);
}
};
diff --git a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
index a17f1feb5ab9af5b100b1d032063f53a0fd844ea..3087971bc11c93c533a3999d3ede8554a624424d 100644
--- a/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/KrylovRunner.hpp
@@ -29,81 +29,80 @@ namespace AMDiS
using Vector = typename Traits::Vec;
using Comm = typename Traits::Comm;
using PreconBase = PreconditionerInterface<Traits>;
+ using FloatType = typename Traits::CoefficientType;
public:
/// Constructor.
explicit KrylovRunner(std::string const& prefix)
: itlSolver_(prefix)
+ , prefix_(prefix)
{
- initPrecon(prefix);
-
- Parameters::get(prefix + "->max iteration", maxIter_);
- Parameters::get(prefix + "->print cycle", printCycle_);
+ Parameters::get(prefix_ + "->absolute tolerance", aTol_);
+ Parameters::get(prefix_ + "->relative tolerance", rTol_);
+ Parameters::get(prefix_ + "->max iteration", maxIter_);
+ Parameters::get(prefix_ + "->print cycle", printCycle_);
}
/// Implementation of \ref RunnerInterface::init()
void init(Matrix const& A, Comm& comm) override
{
DUNE_UNUSED_PARAMETER(comm);
- L_->init(A);
- R_->init(A);
+ initPrecon(prefix_);
+ P_->init(A);
}
/// Implementation of \ref RunnerInterface::exit()
void exit() override
{
- L_->exit();
- R_->exit();
+ P_->exit();
}
/// Implementation of \ref RunnerInterface::solve()
int solve(Matrix const& A, Vector& x, Vector const& b, SolverInfo& solverInfo) override
{
- test_exit(bool(L_), "There is no left preconditioner");
- test_exit(bool(R_), "There is no right preconditioner");
-
- auto r = Vector(b);
- if (two_norm(x) != 0)
- r -= A * x;
+ test_exit(bool(P_), "There is no preconditioner");
// print information about the solution process
- using T = typename Dune::FieldTraits<typename Vector::value_type>::real_type;
- itl::cyclic_iteration<T> iter(two_norm(r),
- maxIter_, solverInfo.relTolerance(), solverInfo.absTolerance(), printCycle_);
+ itl::cyclic_iteration<FloatType> iter(two_norm(b - A*x),
+ maxIter_, rTol_, aTol_, printCycle_);
iter.set_quite(solverInfo.info() == 0);
iter.suppress_resume(solverInfo.info() == 0);
- int error = itlSolver_(A, x, b, *L_, *R_, iter);
+ int error = itlSolver_(A, x, b, *P_, iter);
solverInfo.setAbsResidual(iter.resid());
solverInfo.setRelResidual(iter.relresid());
return error;
}
- private:
+ protected:
/// Create left/right preconditioners from parameters given in the init-file
void initPrecon(std::string const& prefix)
{
// Creator for the left preconditioner
- std::string leftPreconName = "", rightPreconName = "";
- Parameters::get(prefix + "->left precon", leftPreconName);
- Parameters::get(prefix + "->right precon", rightPreconName);
-
- auto leftCreator
- = CreatorMap<PreconBase>::getCreator(leftPreconName, prefix + "->left precon");
- auto rightCreator
- = CreatorMap<PreconBase>::getCreator(rightPreconName, prefix + "->right precon");
+ std::string preconName = "default";
+ Parameters::get(prefix + "->precon", preconName);
- L_ = leftCreator->create();
- R_ = rightCreator->create();
+ auto creator
+ = named(CreatorMap<PreconBase>::getCreator(preconName, prefix + "->precon"));
+ P_ = creator->createWithString(prefix + "->precon");
}
private:
/// The itl solver object that performes the actual solve
ITLSolver itlSolver_;
- /// Pointer to the left and right preconditioner
- std::shared_ptr<PreconBase> L_, R_;
+ /// The parameter prefix
+ std::string prefix_;
+
+ /// Pointer to the left/right preconditioner
+ std::shared_ptr<PreconBase> P_;
+
+ /// Absolute solver tolerance ||b - A*x|| < aTol
+ FloatType aTol_ = 0;
+
+ /// Relative solver tolerance ||b - A*x||/||b - A*x0|| < rTol
+ FloatType rTol_ = 1.e-6;
/// The maximal number of iterations
std::size_t maxIter_ = 1000;
diff --git a/src/amdis/linearalgebra/mtl/Preconditioner.hpp b/src/amdis/linearalgebra/mtl/Preconditioner.hpp
index f8bd006fe559c67f128e12a238c07d8e868e90bd..906b300f0d9b675b395f0df22eb1165ff7cf6fe0 100644
--- a/src/amdis/linearalgebra/mtl/Preconditioner.hpp
+++ b/src/amdis/linearalgebra/mtl/Preconditioner.hpp
@@ -20,9 +20,9 @@ namespace AMDiS
public:
/// A creator to be used instead of the constructor.
- struct Creator : CreatorInterface<Super>
+ struct Creator : CreatorInterfaceName<Super>
{
- std::unique_ptr<Super> create() override
+ std::unique_ptr<Super> createWithString(std::string) override
{
return std::make_unique<Self>();
}
@@ -42,17 +42,17 @@ namespace AMDiS
}
/// Implementation of \ref PreconditionerBase::solve()
- void solve(Vector const& vin, Vector& vout) const override
+ void solve(Vector const& x, Vector& y) const override
{
test_exit_dbg(bool(precon_), "No preconditioner initialized!");
- precon_->solve(vin, vout);
+ precon_->solve(x, y);
}
/// Implementation of \ref PreconditionerBase::adjoint_solve()
- void adjoint_solve(Vector const& vin, Vector& vout) const override
+ void adjoint_solve(Vector const& x, Vector& y) const override
{
test_exit_dbg(bool(precon_), "No preconditioner initialized!");
- precon_->adjoint_solve(vin, vout);
+ precon_->adjoint_solve(x, y);
}
private:
diff --git a/src/amdis/linearalgebra/mtl/SolverPrecon.hpp b/src/amdis/linearalgebra/mtl/SolverPrecon.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..68b6fefadd9be890d707f5d4257507d9d6f6c558
--- /dev/null
+++ b/src/amdis/linearalgebra/mtl/SolverPrecon.hpp
@@ -0,0 +1,92 @@
+#pragma once
+
+#include <memory>
+
+#include <amdis/CreatorMap.hpp>
+#include <amdis/Initfile.hpp>
+#include <amdis/Output.hpp>
+#include <amdis/linearalgebra/LinearSolverInterface.hpp>
+#include <amdis/linearalgebra/PreconditionerInterface.hpp>
+
+namespace AMDiS
+{
+ /**
+ * \ingroup Solver
+ *
+ * \brief Use a LinearSolver as Preconditioner
+ */
+ template <class Traits>
+ class SolverPrecon
+ : public PreconditionerInterface<Traits>
+ {
+ using Self = SolverPrecon;
+ using Matrix = typename Traits::Mat;
+ using Vector = typename Traits::Vec;
+ using Super = PreconditionerInterface<Traits>;
+
+ public:
+ /// A creator to be used instead of the constructor.
+ struct Creator : CreatorInterfaceName<Super>
+ {
+ std::unique_ptr<Super> createWithString(std::string prefix) override
+ {
+ return std::make_unique<Self>(std::move(prefix));
+ }
+ };
+
+ public:
+ SolverPrecon(std::string const& prefix)
+ : solverInfo_(prefix + "->solver")
+ , comm_(Environment::comm())
+ {
+ std::string solverName = "default";
+ Parameters::get(prefix + "->solver", solverName);
+
+ auto solverCreator = named(CreatorMap<LinearSolverInterface<Traits>>::getCreator(solverName, prefix + "->solver"));
+ solver_ = solverCreator->createWithString(prefix + "->solver");
+ }
+
+ /// Implementation of \ref PreconditionerBase::init()
+ void init(Matrix const& fullMatrix) override
+ {
+ matrix_ = &fullMatrix;
+ }
+
+ /// Implementation of \ref PreconditionerInterface::exit()
+ void exit() override
+ {
+ matrix_ = nullptr;
+ }
+
+ /// Implementation of \ref PreconditionerBase::solve()
+ void solve(Vector const& x, Vector& y) const override
+ {
+ test_exit_dbg(bool(solver_), "No solver initialized!");
+ test_exit_dbg(bool(matrix_), "No matrix initialized!");
+
+ y.checked_change_resource(x);
+ test_exit(size(y) == num_cols(*matrix_), "incompatible size");
+ solver_->solve(*matrix_, y, x, comm_, solverInfo_);
+ }
+
+ /// Implementation of \ref PreconditionerBase::adjoint_solve()
+ void adjoint_solve(Vector const& x, Vector& y) const override
+ {
+ error_exit("Not Implemented.");
+ test_exit_dbg(bool(solver_), "No solver initialized!");
+ test_exit_dbg(bool(matrix_), "No matrix initialized!");
+
+ y.checked_change_resource(x);
+ test_exit(size(y) == num_rows(*matrix_), "incompatible size");
+ // solver_->adjointSolve(*matrix_, y, x, comm_, solverInfo_);
+ }
+
+ private:
+ mutable SolverInfo solverInfo_;
+ mutable typename Traits::Comm comm_;
+
+ std::shared_ptr<LinearSolverInterface<Traits>> solver_;
+ Matrix const* matrix_ = nullptr;
+ };
+
+} // namespace AMDiS
diff --git a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
index 643cfbaf865402ec4b316768547ac191c1ccc4c5..4703b6ca4670888e87d5998c1705635657955d20 100644
--- a/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
+++ b/src/amdis/linearalgebra/mtl/UmfpackRunner.hpp
@@ -72,11 +72,7 @@ namespace AMDiS
umfpack_error_exit("solve", e.code);
}
- auto r = Vector(b);
- if (two_norm(x) != 0)
- r -= A * x;
-
- solverInfo.setAbsResidual(two_norm(r));
+ solverInfo.setAbsResidual(two_norm(b - A*x));
solverInfo.setError(code);
return code;
diff --git a/src/amdis/linearalgebra/mtl/itl/CMakeLists.txt b/src/amdis/linearalgebra/mtl/itl/CMakeLists.txt
index ffe212d96d8944dd67b60cd142df139215a5cdac..4df82e7b979899cc047b0bbe00d9fa3883efc111 100644
--- a/src/amdis/linearalgebra/mtl/itl/CMakeLists.txt
+++ b/src/amdis/linearalgebra/mtl/itl/CMakeLists.txt
@@ -1,5 +1,4 @@
install(FILES
- block_diagonal.hpp
details.hpp
fgmres_householder.hpp
fgmres.hpp
@@ -10,5 +9,4 @@ install(FILES
masslumping.hpp
minres.hpp
preonly.hpp
- umfpack2_solve.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/mtl/itl)
diff --git a/src/amdis/linearalgebra/mtl/itl/block_diagonal.hpp b/src/amdis/linearalgebra/mtl/itl/block_diagonal.hpp
deleted file mode 100644
index b295a389e6e8a5dc83ecc3ea125a6675efd64613..0000000000000000000000000000000000000000
--- a/src/amdis/linearalgebra/mtl/itl/block_diagonal.hpp
+++ /dev/null
@@ -1,81 +0,0 @@
-#pragma once
-
-#include <boost/numeric/itl/pc/diagonal.hpp>
-
-#include <amdis/linearalgebra/mtl/BlockMTLMatrix.hpp>
-
-namespace itl
-{
- namespace pc
- {
- /// Diagonal Preconditioner
- template <class MTLMatrix, std::size_t _N, std::size_t _M, class Value>
- class diagonal<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>, Value>
- {
- using Self = diagonal;
- using MatrixType = AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>;
-
- public:
- using value_type = Value;
- using size_type = typename mtl::Collection<MatrixType>::size_type;
-
- /// Constructor takes matrix reference
- explicit diagonal(const MatrixType& A) : inv_diag(num_rows(A))
- {
- MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
- using math::reciprocal;
-
- std::array<mtl::irange, _N> r_rows;
- A.getRowRanges(r_rows);
-
- for (std::size_t i = 0; i < _N; i++)
- inv_diag[r_rows[i]] = mtl::mat::diagonal(A[i][i]);
-
- for (size_type i = 0; i < num_rows(A); ++i)
- inv_diag[i]= reciprocal(inv_diag[i]);
- }
-
- /// Member function solve, better use free function solve
- template <typename Vector>
- Vector solve(const Vector& x) const
- {
- Vector y(resource(x));
- solve(x, y);
- return y;
- }
-
- template <typename VectorIn, typename VectorOut>
- void solve(const VectorIn& x, VectorOut& y) const
- {
- y.checked_change_resource(x);
- MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
- for (size_type i= 0; i < size(inv_diag); ++i)
- y[i]= inv_diag[i] * x[i];
- }
-
- /// Member function for solving adjoint problem, better use free function adjoint_solve
- template <typename Vector>
- Vector adjoint_solve(const Vector& x) const
- {
- Vector y(resource(x));
- adjoint_solve(x, y);
- return y;
- }
-
- template <typename VectorIn, typename VectorOut>
- void adjoint_solve(const VectorIn& x, VectorOut& y) const
- {
- using mtl::conj;
- y.checked_change_resource(x);
- MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
- for (size_type i= 0; i < size(inv_diag); ++i)
- y[i]= conj(inv_diag[i]) * x[i];
- }
-
- protected:
- mtl::vec::dense_vector<value_type> inv_diag;
- };
-
-
- }
-} // namespace itl::pc
diff --git a/src/amdis/linearalgebra/mtl/itl/details.hpp b/src/amdis/linearalgebra/mtl/itl/details.hpp
index 310b19fcdf518aa7dedcbe5f28e8185fe13cfb16..878ecada7320f7256bb98dc6cc767c3806e04d7e 100644
--- a/src/amdis/linearalgebra/mtl/itl/details.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/details.hpp
@@ -1,37 +1,11 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius
-
-
-#ifndef ITL_DETAIL_INCLUDE
-#define ITL_DETAIL_INCLUDE
+#pragma once
#include <boost/math/special_functions/sign.hpp>
namespace itl
{
-
namespace details
{
-
/// Compute the Givens rotation matrix parameters for a and b.
// template<typename T>
// void rotmat(const T& a, const T& b , T& c, T& s)
@@ -40,12 +14,12 @@ namespace itl
//
// const T zero = math::zero(T());
// if (a == zero) {
- // c = 0.0;
- // s = 1.0;
+ // c = 0.0;
+ // s = 1.0;
// } else {
- // double temp = abs(a) / sqrt( conj(a)*a + conj(b)*b );
- // c = temp;
- // s = temp * (b / a);
+ // double temp = abs(a) / sqrt( conj(a)*a + conj(b)*b );
+ // c = temp;
+ // s = temp * (b / a);
// }
// }
@@ -72,7 +46,6 @@ namespace itl
s = temp * c;
}
}
- }
-}
-#endif // ITL_DETAIL_INCLUDE
+ } // end namespace details
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/fgmres.hpp b/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
index 21e4216b4b43060d870edc37107dd68587756a3c..ebf3a37c836ff6325258685fde89c8ab56238d5e 100644
--- a/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/fgmres.hpp
@@ -1,30 +1,7 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius
-
-
-#ifndef ITL_FGMRES_INCLUDE
-#define ITL_FGMRES_INCLUDE
+#pragma once
#include <algorithm>
+
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
@@ -38,7 +15,6 @@
namespace itl
{
-
/// Flexible Generalized Minimal Residual method (without restart)
/// Cite: Youcef Saad, A Flexible Inner-Outer Preconditioned GMRES Algorithm, 1993
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
@@ -193,10 +169,4 @@ namespace itl
return iter;
}
-
-
-} // namespace itl
-
-#endif // ITL_FGMRES_INCLUDE
-
-
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp b/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
index 2f066ee7472f46b3f0e32b010cee9be92dd9a6a6..92f33adff72e9f97c2e6ddaf823d43ac660dbf5e 100644
--- a/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/fgmres_householder.hpp
@@ -1,30 +1,7 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius (adopted from previous implementation)
-
-
-#ifndef ITL_FGMRES_HOUSEHOLDER_INCLUDE
-#define ITL_FGMRES_HOUSEHOLDER_INCLUDE
+#pragma once
#include <algorithm>
+
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
@@ -38,7 +15,6 @@
namespace itl
{
-
/// Flexible Generalized Minimal Residual method (without restart) using householder othogonalization
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
regardless on whether the termination criterion is reached or not. **/
@@ -207,10 +183,4 @@ namespace itl
return iter;
}
-
-
-} // namespace itl
-
-#endif // ITL_FGMRES_HOUSEHOLDER_INCLUDE
-
-
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/gcr.hpp b/src/amdis/linearalgebra/mtl/itl/gcr.hpp
index 4e56b148979ed7310270a48930cd0bbe4863ae2d..dccecde44d0410de2f1569e103008c3a25749f58 100644
--- a/src/amdis/linearalgebra/mtl/itl/gcr.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gcr.hpp
@@ -1,34 +1,9 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius
-
-
-#ifndef ITL_GCR_INCLUDE
-#define ITL_GCR_INCLUDE
+#pragma once
#include <boost/numeric/mtl/concept/collection.hpp>
namespace itl
{
-
template <typename Matrix, typename Vector, typename Preconditioner, typename Iteration>
int gcr_full(const Matrix& A, Vector& x, const Vector& b, const Preconditioner& P, Iteration& iter)
{
@@ -117,7 +92,5 @@ Cite:
return iter;
}
-} // namespace itl;
-
-#endif // ITL_GMR_INCLUDE
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/gmres2.hpp b/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
index c65c7afe97c7dabd7393d48ad365b4dba2ea09bd..4ecbe564f444475e4562d6e78f5f18ba0ad821d1 100644
--- a/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gmres2.hpp
@@ -1,30 +1,7 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius (adopted from previous implementation)
-
-
-#ifndef ITL_GMRES2_INCLUDE
-#define ITL_GMRES2_INCLUDE
+#pragma once
#include <algorithm>
+
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
@@ -38,7 +15,6 @@
namespace itl
{
-
/// Generalized Minimal Residual method (without restart)
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
regardless on whether the termination criterion is reached or not. **/
@@ -196,10 +172,4 @@ namespace itl
return iter;
}
-
-
-} // namespace itl
-
-#endif // ITL_GMRES2_INCLUDE
-
-
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp b/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
index eea759692e35db89045145db006892faef69998a..b2565ca32103ad3477736563921a4405b735d049 100644
--- a/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/gmres_householder.hpp
@@ -1,30 +1,7 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius (adopted from previous implementation)
-
-
-#ifndef ITL_GMRES_HOUSEHOLDER_INCLUDE
-#define ITL_GMRES_HOUSEHOLDER_INCLUDE
+#pragma once
#include <algorithm>
+
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
@@ -38,7 +15,6 @@
namespace itl
{
-
/// Generalized Minimal Residual method (without restart) using householder othogonalization.
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
regardless on whether the termination criterion is reached or not. **/
@@ -212,10 +188,4 @@ namespace itl
return iter;
}
-
-
-} // namespace itl
-
-#endif // ITL_GMRES_HOUSEHOLDER_INCLUDE
-
-
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/hypre.hpp b/src/amdis/linearalgebra/mtl/itl/hypre.hpp
index b877658206fd6482e3d23cd8191ec0310e99123c..0c4540708dda77e139552147051abeca23908d13 100644
--- a/src/amdis/linearalgebra/mtl/itl/hypre.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/hypre.hpp
@@ -1,38 +1,13 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
+#pragma once
-// Written by Andreas Naumann
+#if MTL_HAS_HYPRE && HAVE_MPI
-
-#ifndef AMDIS_HYPRE_SOLVER_HPP
-#define AMDIS_HYPRE_SOLVER_HPP
-
-#ifdef MTL_HAS_HYPRE
-
-#include "mpi.h"
-#include "HYPRE.h"
-#include "HYPRE_parcsr_ls.h"
+#include <mpi.h>
+#include <HYPRE.h>
+#include <HYPRE_parcsr_ls.h>
namespace mtl
{
-
class HypreMatrix
{
public:
@@ -278,8 +253,6 @@ namespace mtl
int _maxIter;
};
-} //namespace mtl
+} // end namespace mtl
#endif // MTL_HAS_HYPRE
-
-#endif // HYPRE_SOLVER_HPP
diff --git a/src/amdis/linearalgebra/mtl/itl/masslumping.hpp b/src/amdis/linearalgebra/mtl/itl/masslumping.hpp
index 26d6b684adfbc2013db1bae540108b9c360ad7d5..3754cc87a3bc56f7d17cceacb40e7aa54f251fd0 100644
--- a/src/amdis/linearalgebra/mtl/itl/masslumping.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/masslumping.hpp
@@ -1,17 +1,4 @@
-// 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_MASSLUMPING_INCLUDE
-#define ITL_PC_MASSLUMPING_INCLUDE
+#pragma once
#include <boost/numeric/linear_algebra/inverse.hpp>
@@ -28,7 +15,6 @@ namespace itl
{
namespace pc
{
-
/// Diagonal Preconditioner
template <typename Matrix, typename Value= typename mtl::Collection<Matrix>::value_type>
class masslumping
@@ -41,7 +27,6 @@ namespace itl
/// Constructor takes matrix reference
explicit masslumping(const Matrix& A) : inv_diag(num_rows(A))
{
- // mtl::vampir_trace<5050> tracer;
MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
using namespace mtl;
using namespace mtl::tag;
@@ -127,8 +112,5 @@ namespace itl
return solver<masslumping<Matrix>, Vector, true>(P, x);
}
-
- }
-} // namespace itl::pc
-
-#endif // ITL_PC_MASSLUMPING_INCLUDE
+ } // end namespace pc
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/minres.hpp b/src/amdis/linearalgebra/mtl/itl/minres.hpp
index 1c1319eea3cfbf7a878f72e66ea8f3ec5a56e382..8d18154117732989dbe22ba2851e946a91e5b76d 100644
--- a/src/amdis/linearalgebra/mtl/itl/minres.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/minres.hpp
@@ -1,55 +1,28 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Thomas Witkowski
-
-
-#ifndef AMDIS_ITL_MINRES_INCLUDE
-#define AMDIS_ITL_MINRES_INCLUDE
+#pragma once
#include <boost/numeric/mtl/concept/collection.hpp>
namespace itl
{
-
/// Minimal Residual method
template <typename Matrix, typename Vector,
- typename LeftPreconditioner, typename RightPreconditioner,
- typename Iteration>
+ typename LeftPreconditioner, typename Iteration>
int minres(const Matrix& A, Vector& x, const Vector& b,
- const LeftPreconditioner& L, const RightPreconditioner& /*R*/,
- Iteration& iter)
+ const LeftPreconditioner& P, Iteration& iter)
{
- using std::abs;
+ using std::abs; using std::sqrt;
using math::reciprocal;
typedef typename mtl::Collection<Vector>::value_type Scalar;
if (size(b) == 0)
throw mtl::logic_error("empty rhs vector");
- Scalar zero= math::zero(b[0]), one= math::one(b[0]);
- Vector v0(size(x), zero), v1(b - A * x), v2(v1), z1(solve(L, v1)), z2(size(x), zero);
+ Scalar zero= math::zero(b[0]), one= math::one(b[0]);
+ Vector v0(size(x), zero), v1(b - A * x), v2(v1), z1(solve(P, v1)), z2(size(x), zero);
Vector w0(size(x), zero), w1(size(x), zero), w2(size(x), zero);
Scalar s0(zero), s1(zero), c0(one), c1(one), gamma0(one);
- Scalar gamma1(sqrt(dot(z1, v1))), gamma2(zero), eta(gamma1);
+ Scalar gamma1(sqrt(abs(dot(z1, v1)))), gamma2(zero), eta(gamma1);
Scalar sigma1(one), alpha0(zero), alpha1(zero), alpha2(zero), alpha3(zero);
while (!iter.finished(abs(eta)))
@@ -59,9 +32,9 @@ namespace itl
sigma1 = dot(v2, z1);
v2 += -(sigma1 / gamma1) * v1 - (gamma1 / gamma0) * v0;
- z2 = solve(L, v2);
+ z2 = solve(P, v2);
- gamma2 = sqrt(dot(z2, v2));
+ gamma2 = sqrt(abs(dot(z2, v2)));
alpha0 = c1 * sigma1 - c0 * s1 * gamma1;
alpha1 = sqrt(alpha0 * alpha0 + gamma2 * gamma2);
alpha2 = s1 * sigma1 + c0 * c1 * gamma1;
@@ -93,6 +66,4 @@ namespace itl
return iter;
}
-} // namespace itl;
-
-#endif // AMDIS_ITL_MINRES_INCLUDE
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/preonly.hpp b/src/amdis/linearalgebra/mtl/itl/preonly.hpp
index 98f1b6c56d9ac113c3b74598e7e0b5305eb3c2f0..0a926f1bd5fe8a535e31eef3ae1b68106458e0f2 100644
--- a/src/amdis/linearalgebra/mtl/itl/preonly.hpp
+++ b/src/amdis/linearalgebra/mtl/itl/preonly.hpp
@@ -1,40 +1,15 @@
-/******************************************************************************
- *
- * AMDiS - Adaptive multidimensional simulations
- *
- * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
- * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
- *
- * Authors:
- * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
- *
- * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
- * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- *
- *
- * This file is part of AMDiS
- *
- * See also license.opensource.txt in the distribution.
- *
- ******************************************************************************/
-
-// Written by Simon Praetorius
-
-
-#ifndef ITL_PREONLY_INCLUDE
-#define ITL_PREONLY_INCLUDE
+#pragma once
+#include <cassert>
#include <boost/numeric/mtl/concept/collection.hpp>
namespace itl
{
-
/// Solver that simply applies a preconditioner to the rhs vector
template <typename Matrix, typename Vector, typename Preconditioner, typename Iteration>
int preonly(const Matrix& A, Vector& x, const Vector& b, const Preconditioner& P, Iteration& iter)
{
- if (size(b) == 0)
- throw mtl::logic_error("empty rhs vector");
+ assert(size(b) != 0);
typedef typename mtl::Collection<Vector>::value_type Scalar;
@@ -50,7 +25,4 @@ namespace itl
return iter;
}
-} // namespace itl
-
-#endif // ITL_PREONLY_INCLUDE
-
+} // end namespace itl
diff --git a/src/amdis/linearalgebra/mtl/itl/umfpack2_solve.hpp b/src/amdis/linearalgebra/mtl/itl/umfpack2_solve.hpp
deleted file mode 100644
index 876bb4edc0a696a77800525a02d9cb66e5efbfc0..0000000000000000000000000000000000000000
--- a/src/amdis/linearalgebra/mtl/itl/umfpack2_solve.hpp
+++ /dev/null
@@ -1,637 +0,0 @@
-// 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 MTL_MATRIX_UMFPACK2_SOLVE_INCLUDE
-#define MTL_MATRIX_UMFPACK2_SOLVE_INCLUDE
-
-#ifdef MTL_HAS_UMFPACK
-
-#include <iostream>
-
-
-#include <cassert>
-#include <algorithm>
-#include <boost/mpl/bool.hpp>
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-#include <boost/numeric/mtl/matrix/parameter.hpp>
-#include <boost/numeric/mtl/concept/collection.hpp>
-#include <boost/numeric/mtl/utility/exception.hpp>
-#include <boost/numeric/mtl/utility/make_copy_or_reference.hpp>
-#include <boost/numeric/mtl/operation/merge_complex_vector.hpp>
-#include <boost/numeric/mtl/operation/split_complex_vector.hpp>
-#include <boost/numeric/mtl/interface/vpt.hpp>
-
-extern "C" {
-# include <umfpack.h>
-}
-
-namespace mtl
-{
- namespace mat
- {
-
- /// Namespace for Umfpack solver
- namespace umfpack
- {
-
- // conversion for value_type needed if not double or complex<double> (where possible)
- template <typename Value> struct value {};
- template<> struct value<long double>
- {
- typedef double type;
- };
- template<> struct value<double>
- {
- typedef double type;
- };
- template<> struct value<float>
- {
- typedef double type;
- };
- template<> struct value<std::complex<long double>>
- {
- typedef std::complex<double> type;
- };
- template<> struct value<std::complex<double>>
- {
- typedef std::complex<double> type;
- };
- template<> struct value<std::complex<float>>
- {
- typedef std::complex<double> type;
- };
-
- template <typename Value> struct use_long
- {
- static const bool value= sizeof(Value) > sizeof(int);
- };
-
- template <bool Larger> struct index_aux
- {
- typedef int type;
- };
-#if defined(UF_long)
- template<> struct index_aux<true>
- {
- typedef UF_long type;
- };
-#elif defined(SuiteSparse_long)
- template<> struct index_aux<true>
- {
- typedef SuiteSparse_long type;
- };
-#else
- template<> struct index_aux<true>
- {
- typedef long type;
- };
-#endif
-
- template <typename Value> struct index
- : index_aux<use_long<Value>::value> {};
-
- template <typename Matrix, typename Value, typename Orientation>
- struct matrix_copy {};
-
- // If arbitrary compressed matrix -> copy
- template <typename Value, typename Parameters, typename Orientation>
- struct matrix_copy<compressed2D<Value, Parameters>, Value, Orientation>
- {
- typedef typename value<Value>::type value_type;
- typedef compressed2D<value_type, parameters<col_major>> matrix_type;
- typedef compressed2D<Value, Parameters> in_matrix_type;
-
- matrix_copy(const in_matrix_type& A) : matrix(A) {}
- matrix_type matrix;
- };
-
- struct error : public domain_error
- {
- error(const char* s, int code) : domain_error(s), code(code) {}
- int code;
- };
-
- inline void check(int res, const char* s)
- {
- MTL_THROW_IF(res != UMFPACK_OK, error(s, res));
- }
-
- /// Class for repeated Umfpack solutions
- /** Keeps symbolic and numeric preprocessing. Numeric part can be updated.
- Only defined for compressed2D<double> and compressed2D<complex<double> >. **/
- template <typename T>
- class solver
- {
- public:
- /// Constructor referring to matrix \p A (not changed) and optionally Umfpack's strategy and alloc_init (look for the specializations)
- // \ref solver<compressed2D<double, Parameters> > and \ref solver<compressed2D<std::complex<double>, Parameters> >)
- explicit solver(const T& A) {}
-
- /// Update numeric part, for matrices that kept the sparsity and changed the values
- void update_numeric() {}
-
- /// Update symbolic and numeric part
- void update() {}
-
- /// Solve system A*x == b with matrix passed in constructor
- /** Please note that the order of b and x is different than in solve() !!! **/
- template <typename VectorX, typename VectorB>
- int operator()(VectorX& x, const VectorB& b) const
- {
- return 0;
- }
-
- /// Solve system A*x == b with matrix passed in constructor
- /** Please note that the order of b and x is different than in operator() !!! **/
- template <typename VectorB, typename VectorX>
- int operator()(const VectorB& b, VectorX& x) const
- {
- return 0;
- }
- };
-
- /// Speciatization of solver for \ref matrix::compressed2D with double values
- template <typename Parameters>
- class solver<compressed2D<double, Parameters>>
- {
- typedef double value_type;
- typedef compressed2D<value_type, Parameters> matrix_type;
- typedef typename matrix_type::size_type size_type;
- typedef typename index<size_type>::type index_type;
-
- static const bool copy_indices= sizeof(index_type) != sizeof(size_type),
- long_indices= use_long<size_type>::value;
- typedef boost::mpl::bool_<long_indices> blong;
- typedef boost::mpl::true_ true_;
- typedef boost::mpl::false_ false_;
-
- // typedef parameters<col_major> Parameters;
-
- void assign_pointers()
- {
- if (copy_indices)
- {
- if (Apc == 0) Apc= new index_type[n + 1];
- if (my_nnz != A.nnz() && Aic)
- {
- delete[] Aic;
- Aic= 0;
- }
- if (Aic == 0) Aic= new index_type[A.nnz()];
- std::copy(A.address_major(), A.address_major() + n + 1, Apc);
- std::copy(A.address_minor(), A.address_minor() + A.nnz(), Aic);
- Ap= Apc;
- Ai= Aic;
- }
- else
- {
- Ap= reinterpret_cast<const index_type*>(A.address_major());
- Ai= reinterpret_cast<const index_type*>(A.address_minor());
- }
- Ax= A.address_data();
- }
-
- void init_aux(true_)
- {
- check(umfpack_dl_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info), "Error in dl_symbolic");
- check(umfpack_dl_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info), "Error in dl_numeric");
- }
-
- void init_aux(false_)
- {
- check(umfpack_di_symbolic(n, n, Ap, Ai, Ax, &Symbolic, Control, Info), "Error in di_symbolic");
-#if 0
- std::cout << "=== INFO of umfpack_*_symbolic ===\n";
- std::cout << " UMFPACK_STATUS: " << (Info[UMFPACK_STATUS] == UMFPACK_OK ? "OK" : "ERROR") << "\n";
- std::cout << " UMFPACK_NROW: " << Info[UMFPACK_NROW] << "\n";
- std::cout << " UMFPACK_NCOL: " << Info[UMFPACK_NCOL] << "\n";
- std::cout << " UMFPACK_NZ: " << Info[UMFPACK_NZ] << "\n";
- std::cout << " UMFPACK_SIZE_OF_UNIT: " << Info[UMFPACK_SIZE_OF_UNIT] << "\n";
- std::cout << " UMFPACK_NDENSE_ROW: " << Info[UMFPACK_NDENSE_ROW] << "\n";
- std::cout << " UMFPACK_NEMPTY_ROW: " << Info[UMFPACK_NEMPTY_ROW] << "\n";
- std::cout << " UMFPACK_NDENSE_COL: " << Info[UMFPACK_NDENSE_COL] << "\n";
- std::cout << " UMFPACK_NEMPTY_COL: " << Info[UMFPACK_NEMPTY_COL] << "\n";
- std::cout << " UMFPACK_SYMBOLIC_DEFRAG: " << Info[UMFPACK_SYMBOLIC_DEFRAG] << "\n";
- std::cout << " UMFPACK_SYMBOLIC_PEAK_MEMORY: " << Info[UMFPACK_SYMBOLIC_PEAK_MEMORY] << "\n";
- std::cout << " UMFPACK_SYMBOLIC_SIZE: " << Info[UMFPACK_SYMBOLIC_SIZE] << "\n";
- std::cout << " UMFPACK_VARIABLE_PEAK_ESTIMATE: " << Info[UMFPACK_VARIABLE_PEAK_ESTIMATE] << "\n";
- std::cout << " UMFPACK_NUMERIC_SIZE_ESTIMATE: " << Info[UMFPACK_NUMERIC_SIZE_ESTIMATE] << "\n";
- std::cout << " UMFPACK_PEAK_MEMORY_ESTIMATE: " << Info[UMFPACK_PEAK_MEMORY_ESTIMATE] << "\n";
- std::cout << " UMFPACK_FLOPS_ESTIMATE: " << Info[UMFPACK_FLOPS_ESTIMATE] << "\n";
- std::cout << " UMFPACK_LNZ_ESTIMATE: " << Info[UMFPACK_LNZ_ESTIMATE] << "\n";
- std::cout << " UMFPACK_UNZ_ESTIMATE: " << Info[UMFPACK_UNZ_ESTIMATE] << "\n";
- std::cout << " UMFPACK_MAX_FRONT_SIZE_ESTIMATE: " << Info[UMFPACK_MAX_FRONT_SIZE_ESTIMATE] << "\n";
- std::cout << " UMFPACK_SYMBOLIC_TIME: " << Info[UMFPACK_SYMBOLIC_TIME] << "\n";
- std::cout << " UMFPACK_SYMBOLIC_WALLTIME: " << Info[UMFPACK_SYMBOLIC_WALLTIME] << "\n";
-
- if (Info[UMFPACK_STRATEGY_USED] == UMFPACK_STRATEGY_SYMMETRIC)
- std::cout << " UMFPACK_STRATEGY_USED: SYMMETRIC\n";
- else
- {
- if(Info[UMFPACK_STRATEGY_USED] == UMFPACK_STRATEGY_UNSYMMETRIC)
- std::cout << " UMFPACK_STRATEGY_USED: UNSYMMETRIC\n";
- else
- {
- if (Info[UMFPACK_STRATEGY_USED] == UMFPACK_STRATEGY_2BY2)
- std::cout << " UMFPACK_STRATEGY_USED: 2BY2\n";
- else
- std::cout << " UMFPACK_STRATEGY_USED: UNKOWN STRATEGY " << Info[UMFPACK_STRATEGY_USED] << "\n";
- }
- }
-
- std::cout << " UMFPACK_ORDERING_USED: " << Info[UMFPACK_ORDERING_USED] << "\n";
- std::cout << " UMFPACK_QFIXED: " << Info[UMFPACK_QFIXED] << "\n";
- std::cout << " UMFPACK_DIAG_PREFERRED: " << Info[UMFPACK_DIAG_PREFERRED] << "\n";
- std::cout << " UMFPACK_ROW_SINGLETONS: " << Info[UMFPACK_ROW_SINGLETONS] << "\n";
- std::cout << " UMFPACK_COL_SINGLETONS: " << Info[UMFPACK_COL_SINGLETONS] << "\n";
- std::cout << " UMFPACK_PATTERN_SYMMETRY: " << Info[UMFPACK_PATTERN_SYMMETRY] << "\n";
- std::cout << " UMFPACK_NZ_A_PLUS_AT: " << Info[UMFPACK_NZ_A_PLUS_AT] << "\n";
- std::cout << " UMFPACK_NZDIAG: " << Info[UMFPACK_NZDIAG] << "\n";
- std::cout << " UMFPACK_N2: " << Info[UMFPACK_N2] << "\n";
- std::cout << " UMFPACK_S_SYMMETRIC: " << Info[UMFPACK_S_SYMMETRIC] << "\n";
- std::cout << " UMFPACK_MAX_FRONT_NROWS_ESTIMATE: " << Info[UMFPACK_MAX_FRONT_NROWS_ESTIMATE] << "\n";
- std::cout << " UMFPACK_MAX_FRONT_NCOLS_ESTIMATE: " << Info[UMFPACK_MAX_FRONT_NCOLS_ESTIMATE] << "\n";
- std::cout << " UMFPACK_SYMMETRIC_LUNZ: " << Info[UMFPACK_SYMMETRIC_LUNZ] << "\n";
- std::cout << " UMFPACK_SYMMETRIC_FLOPS: " << Info[UMFPACK_SYMMETRIC_FLOPS] << "\n";
- std::cout << " UMFPACK_SYMMETRIC_NDENSE: " << Info[UMFPACK_SYMMETRIC_NDENSE] << "\n";
- std::cout << " UMFPACK_SYMMETRIC_DMAX: " << Info[UMFPACK_SYMMETRIC_DMAX] << "\n";
-#endif
-
- check(umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info), "Error in di_numeric");
-
-#if 0
- std::cout << "=== INFO of umfpack_*_numeric ===\n";
- std::cout << " UMFPACK_STATUS: " << (Info[UMFPACK_STATUS] == UMFPACK_OK ? "OK" : "ERROR") << "\n";
- std::cout << " UMFPACK_VARIABLE_PEAK: " << Info[UMFPACK_VARIABLE_PEAK] << "\n";
- std::cout << " UMFPACK_PEAK_MEMORY: " << Info[UMFPACK_PEAK_MEMORY] << "\n";
- std::cout << " UMFPACK_FLOPS: " << Info[UMFPACK_FLOPS] << "\n";
- std::cout << " UMFPACK_LNZ: " << Info[UMFPACK_LNZ] << "\n";
- std::cout << " UMFPACK_UNZ: " << Info[UMFPACK_UNZ] << "\n";
- std::cout << " UMFPACK_NUMERIC_DEFRAG: " << Info[UMFPACK_NUMERIC_DEFRAG] << "\n";
- std::cout << " UMFPACK_NUMERIC_REALLOC: " << Info[UMFPACK_NUMERIC_REALLOC] << "\n";
- std::cout << " UMFPACK_NUMERIC_COSTLY_REALLOC: " << Info[UMFPACK_NUMERIC_COSTLY_REALLOC] << "\n";
- std::cout << " UMFPACK_COMPRESSED_PATTERN: " << Info[UMFPACK_COMPRESSED_PATTERN] << "\n";
- std::cout << " UMFPACK_LU_ENTRIES: " << Info[UMFPACK_LU_ENTRIES] << "\n";
- std::cout << " UMFPACK_NUMERIC_TIME: " << Info[UMFPACK_NUMERIC_TIME] << "\n";
- std::cout << " UMFPACK_RCOND: " << Info[UMFPACK_RCOND] << "\n";
- std::cout << " UMFPACK_UDIAG_NZ: " << Info[UMFPACK_UDIAG_NZ] << "\n";
- std::cout << " UMFPACK_UMIN: " << Info[UMFPACK_UMIN] << "\n";
- std::cout << " UMFPACK_UMAX: " << Info[UMFPACK_UMAX] << "\n";
- std::cout << " UMFPACK_MAX_FRONT_NROWS: " << Info[UMFPACK_MAX_FRONT_NROWS] << "\n";
- std::cout << " UMFPACK_MAX_FRONT_NCOLS: " << Info[UMFPACK_MAX_FRONT_NCOLS] << "\n";
- std::cout << " UMFPACK_ALL_LNZ: " << Info[UMFPACK_ALL_LNZ] << "\n";
- std::cout << " UMFPACK_ALL_UNZ: " << Info[UMFPACK_ALL_UNZ] << "\n";
-#endif
- }
-
- void init()
- {
- MTL_THROW_IF(num_rows(A) != num_cols(A), matrix_not_square());
- n= num_rows(A);
- assign_pointers();
- init_aux(blong());
- }
-
-
-
- public:
- /// Constructor referring to matrix \p A (not changed) and optionally Umfpack's strategy and alloc_init
- solver(const matrix_type& A, int strategy = UMFPACK_STRATEGY_AUTO, double alloc_init = 0.7)
- : A(A), Apc(0), Aic(0), my_nnz(0), Symbolic(0), Numeric(0)
- {
- vampir_trace<5060> trace;
- // Use default setings.
- if (long_indices)
- umfpack_dl_defaults(Control);
- else
- umfpack_di_defaults(Control);
-
- Control[UMFPACK_STRATEGY] = strategy;
- Control[UMFPACK_ALLOC_INIT] = alloc_init;
- init();
- }
-
- ~solver()
- {
- vampir_trace<5061> trace;
- if (long_indices)
- {
- umfpack_dl_free_numeric(&Numeric);
- umfpack_dl_free_symbolic(&Symbolic);
- }
- else
- {
- umfpack_di_free_numeric(&Numeric);
- umfpack_di_free_symbolic(&Symbolic);
- }
- if (Apc) delete[] Apc;
- if (Aic) delete[] Aic;
- }
-
- void update_numeric_aux(true_)
- {
- umfpack_dl_free_numeric(&Numeric);
- check(umfpack_dl_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info), "Error in dl_numeric");
- }
-
- void update_numeric_aux(false_)
- {
- umfpack_di_free_numeric(&Numeric);
- check(umfpack_di_numeric(Ap, Ai, Ax, Symbolic, &Numeric, Control, Info), "Error in di_numeric");
- }
-
- /// Update numeric part, for matrices that kept the sparsity and changed the values
- void update_numeric()
- {
- assign_pointers();
- update_numeric_aux(blong());
- }
-
- /// Update symbolic and numeric part
- void update()
- {
- if (long_indices)
- {
- umfpack_dl_free_numeric(&Numeric);
- umfpack_dl_free_symbolic(&Symbolic);
- }
- else
- {
- umfpack_di_free_numeric(&Numeric);
- umfpack_di_free_symbolic(&Symbolic);
- }
- init();
- }
-
- template <typename VectorX, typename VectorB>
- void solve_aux(int sys, VectorX& xx, const VectorB& bb, true_)
- {
- check(umfpack_dl_solve(sys, Ap, Ai, Ax, &xx.value[0], &bb.value[0], Numeric, Control, Info), "Error in dl_solve");
- }
-
- template <typename VectorX, typename VectorB>
- void solve_aux(int sys, VectorX& xx, const VectorB& bb, false_)
- {
- check(umfpack_di_solve(sys, Ap, Ai, Ax, &xx.value[0], &bb.value[0], Numeric, Control, Info), "Error in di_solve");
- }
-
- /// Solve double system
- template <typename VectorX, typename VectorB>
- int operator()(VectorX& x, const VectorB& b)
- {
- vampir_trace<5062> trace;
- MTL_THROW_IF(num_rows(A) != size(x) || num_rows(A) != size(b), incompatible_size());
- make_in_out_copy_or_reference<dense_vector<value_type>, VectorX> xx(x);
- make_in_copy_or_reference<dense_vector<value_type>, VectorB> bb(b);
- int sys= mtl::traits::is_row_major<Parameters>::value ? UMFPACK_At : UMFPACK_A;
- solve_aux(sys, xx, bb, blong());
- return UMFPACK_OK;
- }
-
- /// Solve double system
- template <typename VectorB, typename VectorX>
- int solve(const VectorB& b, VectorX& x) const
- {
- // return (*this)(x, b);
- return const_cast<solver&>(*this)(x, b); // evil hack because Umfpack has no const
- }
-
- private:
- const matrix_type& A;
- int n;
- const index_type* Ap, *Ai;
- index_type* Apc, *Aic;
- size_type my_nnz;
- const double* Ax;
- double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
- void* Symbolic, *Numeric;
- };
-
- /// Speciatization of solver for \ref matrix::compressed2D with double values
- template <typename Parameters>
- class solver<compressed2D<std::complex<double>, Parameters>>
- {
- typedef std::complex<double> value_type;
- typedef compressed2D<value_type, Parameters> matrix_type;
- typedef typename matrix_type::size_type size_type;
- typedef typename index<size_type>::type index_type;
-
- static const bool copy_indices= sizeof(index_type) != sizeof(size_type),
- long_indices= use_long<size_type>::value;
-
- typedef boost::mpl::bool_<long_indices> blong;
- typedef boost::mpl::true_ true_;
- typedef boost::mpl::false_ false_;
-
-
- void assign_pointers()
- {
- if (copy_indices)
- {
- if (Apc == 0) Apc= new index_type[n + 1];
- if (Aic == 0) Aic= new index_type[A.nnz()];
- std::copy(A.address_major(), A.address_major() + n + 1, Apc);
- std::copy(A.address_minor(), A.address_minor() + A.nnz(), Aic);
- Ap= Apc;
- Ai= Aic;
- }
- else
- {
- Ap= reinterpret_cast<const index_type*>(A.address_major());
- Ai= reinterpret_cast<const index_type*>(A.address_minor());
- }
- split_complex_vector(A.data, Ax, Az);
- }
-
- void init_aux(true_)
- {
- check(umfpack_zl_symbolic(n, n, Ap, Ai, &Ax[0], &Az[0], &Symbolic, Control, Info), "Error in zl_symbolic");
- check(umfpack_zl_numeric(Ap, Ai, &Ax[0], &Az[0], Symbolic, &Numeric, Control, Info), "Error in zl_numeric");
- }
-
- void init_aux(false_)
- {
- check(umfpack_zi_symbolic(n, n, Ap, Ai, &Ax[0], &Az[0], &Symbolic, Control, Info), "Error in zi_symbolic");
- check(umfpack_zi_numeric(Ap, Ai, &Ax[0], &Az[0], Symbolic, &Numeric, Control, Info), "Error in zi_numeric");
- }
-
- void initialize()
- {
- MTL_THROW_IF(num_rows(A) != num_cols(A), matrix_not_square());
- n= num_rows(A);
- assign_pointers();
- init_aux(blong());
- }
- public:
- /// Constructor referring to matrix \p A (not changed) and optionally Umfpack's strategy and alloc_init (look for the specializations)
- explicit solver(const compressed2D<value_type, Parameters>& A, int strategy = UMFPACK_STRATEGY_AUTO, double alloc_init = 0.7)
- : A(A), Apc(0), Aic(0)
- {
- vampir_trace<5060> trace;
- // Use default setings.
- if (long_indices)
- umfpack_zl_defaults(Control);
- else
- umfpack_zi_defaults(Control);
- // umfpack_zi_defaults(Control);
-
- Control[UMFPACK_STRATEGY] = strategy;
- Control[UMFPACK_ALLOC_INIT] = alloc_init;
- initialize();
- }
-
- ~solver()
- {
- vampir_trace<5061> trace;
- if (long_indices)
- {
- umfpack_zl_free_numeric(&Numeric);
- umfpack_zl_free_symbolic(&Symbolic);
- }
- else
- {
- umfpack_zi_free_numeric(&Numeric);
- umfpack_zi_free_symbolic(&Symbolic);
- }
- if (Apc) delete[] Apc;
- if (Aic) delete[] Aic;
- }
-
- void update_numeric_aux(true_)
- {
- umfpack_zl_free_numeric(&Numeric);
- check(umfpack_zl_numeric(Ap, Ai, &Ax[0], &Az[0], Symbolic, &Numeric, Control, Info), "Error in dl_numeric D");
- }
-
- void update_numeric_aux(false_)
- {
- umfpack_zi_free_numeric(&Numeric);
- check(umfpack_zi_numeric(Ap, Ai, &Ax[0], &Az[0], Symbolic, &Numeric, Control, Info), "Error in di_numeric");
- }
-
- /// Update numeric part, for matrices that kept the sparsity and changed the values
- void update_numeric()
- {
- assign_pointers();
- update_numeric_aux(blong());
- }
-
- /// Update symbolic and numeric part
- void update()
- {
- Ax.change_dim(0);
- Az.change_dim(0);
- if (long_indices)
- {
- umfpack_zl_free_numeric(&Numeric);
- umfpack_zl_free_symbolic(&Symbolic);
- }
- else
- {
- umfpack_zi_free_numeric(&Numeric);
- umfpack_zi_free_symbolic(&Symbolic);
- }
- initialize();
- }
-
- template <typename VectorX, typename VectorB>
- void solve_aux(int sys, VectorX& Xx, VectorX& Xz, const VectorB& Bx, const VectorB& Bz, true_)
- {
- check(umfpack_zl_solve(sys, Ap, Ai, &Ax[0], &Az[0], &Xx[0], &Xz[0], &Bx[0], &Bz[0], Numeric, Control, Info),
- "Error in zi_solve");
- }
-
- template <typename VectorX, typename VectorB>
- void solve_aux(int sys, VectorX& Xx, VectorX& Xz, const VectorB& Bx, const VectorB& Bz, false_)
- {
- check(umfpack_zi_solve(sys, Ap, Ai, &Ax[0], &Az[0], &Xx[0], &Xz[0], &Bx[0], &Bz[0], Numeric, Control, Info),
- "Error in zi_solve");
- }
-
- /// Solve complex system
- template <typename VectorX, typename VectorB>
- int operator()(VectorX& x, const VectorB& b)
- {
- vampir_trace<5062> trace;
- MTL_THROW_IF(num_rows(A) != size(x) || num_rows(A) != size(b), incompatible_size());
- dense_vector<double> Xx(size(x)), Xz(size(x)), Bx, Bz;
- split_complex_vector(b, Bx, Bz);
- int sys= mtl::traits::is_row_major<Parameters>::value ? UMFPACK_Aat : UMFPACK_A;
- solve_aux(sys, Xx, Xz, Bx, Bz, blong());
- merge_complex_vector(Xx, Xz, x);
- return UMFPACK_OK;
- }
-
- /// Solve complex system
- template <typename VectorB, typename VectorX>
- int solve(const VectorB& b, VectorX& x)
- {
- return (*this)(x, b);
- }
-
- private:
- const matrix_type& A;
- int n;
- const index_type* Ap, *Ai;
- index_type* Apc, *Aic;
- dense_vector<double> Ax, Az;
- double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
- void* Symbolic, *Numeric;
- };
-
- template <typename Value, typename Parameters>
- class solver<compressed2D<Value, Parameters>>
- : matrix_copy<compressed2D<Value, Parameters>, Value, typename Parameters::orientation>,
- public solver<typename matrix_copy<compressed2D<Value, Parameters>, Value, typename Parameters::orientation>::matrix_type>
- {
- typedef matrix_copy<compressed2D<Value, Parameters>, Value, typename Parameters::orientation> copy_type;
- typedef solver<typename matrix_copy<compressed2D<Value, Parameters>, Value, typename Parameters::orientation>::matrix_type> solver_type;
- public:
- explicit solver(const compressed2D<Value, Parameters>& A)
- : copy_type(A), solver_type(copy_type::matrix), A(A)
- {}
-
- void update()
- {
- copy_type::matrix= A;
- solver_type::update();
- }
-
- void update_numeric()
- {
- copy_type::matrix= A;
- solver_type::update_numeric();
- }
- private:
- const compressed2D<Value, Parameters>& A;
- };
- } // umfpack
-
- /// Solve A*x == b with umfpack
- /** Only available when compiled with enabled macro MTL_HAS_UMFPACK.
- Uses classes umfpack::solver internally.
- If you want more control on single operations or to keep umfpack's
- internal factorization, use this class.
- **/
- template <typename Value, typename Parameters, typename VectorX, typename VectorB>
- int umfpack2_solve(const compressed2D<Value, Parameters>& A, VectorX& x, const VectorB& b)
- {
- umfpack::solver<compressed2D<Value, Parameters>> solver(A);
- return solver(x, b);
- }
-
- }
-} // namespace mtl::mat
-
-#endif
-
-#endif // MTL_MATRIX_UMFPACK_SOLVE_INCLUDE