diff --git a/example/create_geometry.py b/example/create_geometry.py index 698096f..33e27a2 100644 --- a/example/create_geometry.py +++ b/example/create_geometry.py @@ -1,6 +1,3 @@ -import math - -import mpi4py import numpy as np diff --git a/example/define_custom_dense_blocks_generator.py b/example/define_custom_dense_blocks_generator.py index 5dea850..1778d0c 100644 --- a/example/define_custom_dense_blocks_generator.py +++ b/example/define_custom_dense_blocks_generator.py @@ -1,5 +1,3 @@ -import mpi4py - import Htool diff --git a/example/define_custom_generators.py b/example/define_custom_generators.py index 7b65e9e..334e7a6 100644 --- a/example/define_custom_generators.py +++ b/example/define_custom_generators.py @@ -1,7 +1,5 @@ -import mpi4py -import numpy as np - import Htool +import numpy as np class CustomGenerator(Htool.VirtualGenerator): diff --git a/example/define_custom_local_operator.py b/example/define_custom_local_operator.py index eedebce..f1c441b 100644 --- a/example/define_custom_local_operator.py +++ b/example/define_custom_local_operator.py @@ -1,7 +1,5 @@ -import mpi4py -import numpy as np - import Htool +import numpy as np class CustomLocalOperator(Htool.LocalOperator): diff --git a/example/define_custom_low_rank_generator.py b/example/define_custom_low_rank_generator.py index 32a44b2..6b02f5b 100644 --- a/example/define_custom_low_rank_generator.py +++ b/example/define_custom_low_rank_generator.py @@ -1,8 +1,7 @@ import math -import numpy as np - import Htool +import numpy as np class CustomSVD(Htool.VirtualLowRankGenerator): diff --git a/example/use_cluster.py b/example/use_cluster.py index fdc56a1..e975e89 100644 --- a/example/use_cluster.py +++ b/example/use_cluster.py @@ -1,10 +1,8 @@ +import Htool import matplotlib.pyplot as plt import mpi4py -import numpy as np from create_geometry import create_random_geometries -import Htool - # Random geometry nb_rows = 500 nb_cols = 500 diff --git a/example/use_cluster_with_given_partition.py b/example/use_cluster_with_given_partition.py index d66c43d..279e44a 100644 --- a/example/use_cluster_with_given_partition.py +++ b/example/use_cluster_with_given_partition.py @@ -1,6 +1,5 @@ import matplotlib.pyplot as plt import mpi4py -import numpy as np from create_geometry import create_partitionned_geometries import Htool diff --git a/example/use_custom_dense_block_generator.py b/example/use_custom_dense_block_generator.py index 2ee94bc..83fa788 100644 --- a/example/use_custom_dense_block_generator.py +++ b/example/use_custom_dense_block_generator.py @@ -1,4 +1,3 @@ -import matplotlib.pyplot as plt import mpi4py import numpy as np from create_geometry import create_partitionned_geometries diff --git a/example/use_custom_low_rank_approximation.py b/example/use_custom_low_rank_approximation.py index 15cf7b9..1ed76b5 100644 --- a/example/use_custom_low_rank_approximation.py +++ b/example/use_custom_low_rank_approximation.py @@ -1,11 +1,10 @@ +import Htool import mpi4py import numpy as np from create_geometry import create_partitionned_geometries from define_custom_generators import CustomGenerator from define_custom_low_rank_generator import CustomSVD -import Htool - # Random geometry nb_rows = 500 nb_cols = 500 @@ -52,8 +51,11 @@ -1, mpi4py.MPI.COMM_WORLD.rank, ) -# hmatrix_builder.set_low_rank_generator(low_rank_generator) +hmatrix_builder.set_low_rank_generator(low_rank_generator) hmatrix: Htool.HMatrix = hmatrix_builder.build(generator) +del hmatrix_builder +low_rank_generator.test() +del low_rank_generator # Build local operator diff --git a/example/use_default_build.py b/example/use_default_build.py index 86570db..1697cc0 100644 --- a/example/use_default_build.py +++ b/example/use_default_build.py @@ -1,10 +1,9 @@ -import logging import matplotlib.pyplot as plt import mpi4py import numpy as np from create_geometry import create_partitionned_geometries -from define_custom_generators import CustomGenerator, CustomGeneratorWithPermutation +from define_custom_generators import CustomGenerator import Htool diff --git a/example/use_local_build.py b/example/use_local_build.py index a7250a9..9b365c9 100644 --- a/example/use_local_build.py +++ b/example/use_local_build.py @@ -1,9 +1,8 @@ -import logging import matplotlib.pyplot as plt import mpi4py import numpy as np -from create_geometry import create_partitionned_geometries, create_random_geometries +from create_geometry import create_random_geometries from define_custom_generators import CustomGenerator from define_custom_local_operator import CustomLocalOperator diff --git a/lib/htool b/lib/htool index ca6582c..e2105b0 160000 --- a/lib/htool +++ b/lib/htool @@ -1 +1 @@ -Subproject commit ca6582cb4787e3c8c72fb94b72ebd299011293a1 +Subproject commit e2105b070a92c281b9b5f756c040f1764e5ec4ec diff --git a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp index ee11fee..18fa23b 100644 --- a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp +++ b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp @@ -13,7 +13,6 @@ using namespace htool; template class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator { - int m_rank; mutable std::vector> m_mats_U; // owned by Python mutable std::vector> m_mats_V; // owned by Python @@ -35,7 +34,9 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator &A, int target_size, int source_size, int target_offset, int source_offset, underlying_type epsilon) const = 0; // LCOV_EXCL_LINE - void set_U(py::array_t U0) { m_mats_U.push_back(U0); } + void set_U(py::array_t U0) { + m_mats_U.push_back(U0); // no copy here + } void set_V(py::array_t V0) { m_mats_V.push_back(V0); } void set_rank(int rank) { m_rank = rank; } }; diff --git a/src/htool/main.cpp b/src/htool/main.cpp index 37c6cbb..84df0f9 100644 --- a/src/htool/main.cpp +++ b/src/htool/main.cpp @@ -25,12 +25,17 @@ #include "distributed_operator/implementation/partition_from_cluster.hpp" #include "distributed_operator/utility.hpp" +#include "solver/geneo/coarse_operator_builder.hpp" +#include "solver/geneo/coarse_space_builder.hpp" +#include "solver/interfaces/virtual_coarse_operator_builder.hpp" +#include "solver/interfaces/virtual_coarse_space_builder.hpp" #include "solver/solver.hpp" #include "solver/utility.hpp" #include "matplotlib/cluster.hpp" #include "matplotlib/hmatrix.hpp" #include "misc/logger.hpp" +#include "misc/wrapper_mpi.hpp" // #include "ddm_solver.hpp" // #include "dense_blocks_generator.hpp" // #include "hmatrix.hpp" @@ -39,7 +44,7 @@ PYBIND11_MODULE(Htool, m) { // Delegate logging to python logging module - Logger::get_instance().set_current_writer(std::make_shared()); + htool::Logger::get_instance().set_current_writer(std::make_shared()); // import the mpi4py API if (import_mpi4py() < 0) { @@ -72,6 +77,11 @@ PYBIND11_MODULE(Htool, m) { declare_distributed_operator_utility(m); declare_DDM(m, "Solver"); + declare_virtual_coarse_space_builder(m, "VirtualCoarseSpaceBuilder", "ICoarseSpaceBuilder"); + declare_virtual_coarse_operator_builder(m, "", "ICoarseOperatorBuilder"); + declare_geneo_coarse_operator_builder(m, "GeneoCoarseOperatorBuilder"); + declare_geneo_coarse_space_dense_builder(m, "GeneoCoarseSpaceDenseBuilder"); + declare_virtual_geneo_coarse_space_dense_builder(m, "VirtualGeneoCoarseSpaceDenseBuilder"); declare_solver_utility(m); declare_solver_utility(m); @@ -86,5 +96,11 @@ PYBIND11_MODULE(Htool, m) { declare_distributed_operator_utility, double>(m, "Complex"); declare_DDM>(m, "ComplexSolver"); + declare_virtual_coarse_space_builder>(m, "ComplexVirtualCoarseSpaceBuilder", "IComplexCoarseSpaceBuilder"); + declare_virtual_coarse_operator_builder>(m, "", "IComplexCoarseOperatorBuilder"); + declare_geneo_coarse_operator_builder>(m, "ComplexGeneoCoarseOperatorBuilder"); + declare_geneo_coarse_space_dense_builder>(m, "ComplexGeneoCoarseSpaceDenseBuilder"); + declare_virtual_geneo_coarse_space_dense_builder>(m, "VirtualComplexGeneoCoarseSpaceDenseBuilder"); + declare_solver_utility, double>(m, "Complex"); } diff --git a/src/htool/misc/logger.hpp b/src/htool/misc/logger.hpp index ae95a79..3f5d72d 100644 --- a/src/htool/misc/logger.hpp +++ b/src/htool/misc/logger.hpp @@ -3,30 +3,30 @@ #include #include - +namespace py = pybind11; // Delegates formatting and filtering to python logging module // mapping python logging level to Htool logging level -class PythonLoggerWriter : public IObjectWriter { +class PythonLoggerWriter : public htool::IObjectWriter { public: - void write(LogLevel logging_level, const std::string &message) override { + void write(htool::LogLevel logging_level, const std::string &message) override { py::object logging_module = py::module::import("logging"); py::object logger = logging_module.attr("getLogger")("Htool"); switch (logging_level) { - case LogLevel::CRITICAL: + case htool::LogLevel::CRITICAL: logger.attr("critical")(message); break; - case LogLevel::ERROR: + case htool::LogLevel::ERROR: logger.attr("error")(message); break; - case LogLevel::WARNING: + case htool::LogLevel::WARNING: logger.attr("warning")(message); break; - case LogLevel::DEBUG: + case htool::LogLevel::DEBUG: logger.attr("debug")(message); break; - case LogLevel::INFO: + case htool::LogLevel::INFO: logger.attr("info")(message); break; default: @@ -34,7 +34,7 @@ class PythonLoggerWriter : public IObjectWriter { } } - void set_log_level(LogLevel log_level) override { + void set_log_level(htool::LogLevel log_level) override { } }; diff --git a/src/htool/solver/geneo/coarse_operator_builder.hpp b/src/htool/solver/geneo/coarse_operator_builder.hpp new file mode 100644 index 0000000..c16e301 --- /dev/null +++ b/src/htool/solver/geneo/coarse_operator_builder.hpp @@ -0,0 +1,16 @@ +#ifndef HTOOL_PYTHON_GENEO_COARSE_OPERATOR_BUILDER_HPP +#define HTOOL_PYTHON_GENEO_COARSE_OPERATOR_BUILDER_HPP + +#include +#include +namespace py = pybind11; + +template +void declare_geneo_coarse_operator_builder(py::module &m, const std::string &className) { + + using Class = GeneoCoarseOperatorBuilder; + py::class_> py_class(m, className.c_str()); + py_class.def(py::init &>()); +} + +#endif diff --git a/src/htool/solver/geneo/coarse_space_builder.hpp b/src/htool/solver/geneo/coarse_space_builder.hpp new file mode 100644 index 0000000..9d9b7e6 --- /dev/null +++ b/src/htool/solver/geneo/coarse_space_builder.hpp @@ -0,0 +1,131 @@ +#ifndef HTOOL_PYTHON_GENEO_COARSE_SPACE_BUILDER_HPP +#define HTOOL_PYTHON_GENEO_COARSE_SPACE_BUILDER_HPP + +#include +#include +namespace py = pybind11; + +template +class GeneoCoarseSpaceDenseBuilderPython : public GeneoCoarseSpaceDenseBuilder { + py::array_t m_coarse_space; + + using GeneoCoarseSpaceDenseBuilder::GeneoCoarseSpaceDenseBuilder; + + public: + char get_symmetry() { return this->m_symmetry; } + int get_geneo_nu() { return this->m_geneo_nu; } + htool::underlying_type get_geneo_threshold() { return this->m_geneo_threshold; } + + GeneoCoarseSpaceDenseBuilderPython(int size_wo_overlap, Matrix Ai, Matrix Bi, char symmetry, char uplo, int geneo_nu, htool::underlying_type geneo_threshold) : GeneoCoarseSpaceDenseBuilder(size_wo_overlap, Ai, Bi, symmetry, uplo, geneo_nu, geneo_threshold) {} + + Matrix build_coarse_space() override { + py::array_t Ai(std::array{this->m_DAiD.nb_rows(), this->m_DAiD.nb_cols()}, this->m_DAiD.data(), py::capsule(this->m_DAiD.data())); + py::array_t Bi(std::array{this->m_Bi.nb_rows(), this->m_Bi.nb_cols()}, this->m_Bi.data(), py::capsule(this->m_Bi.data())); + compute_coarse_space(Ai, Bi); + Matrix coarse_space_mat(m_coarse_space.shape()[0], m_coarse_space.shape()[1]); + std::copy_n(m_coarse_space.data(), m_coarse_space.shape()[0] * m_coarse_space.shape()[1], coarse_space_mat.data()); // HPDDM deletes the coarse space, so we have to copy. + // coarse_space_mat.assign(m_coarse_space.shape()[0], m_coarse_space.shape()[1], m_coarse_space.mutable_data(), false); + // coarse_space_mat.print(std::cout, ","); + return coarse_space_mat; + } + + virtual void compute_coarse_space(py::array_t Ai, py::array_t Bi) = 0; + + void set_coarse_space(py::array_t coarse_space) { + m_coarse_space = coarse_space; + } +}; + +template +class PyGeneoCoarseSpaceDenseBuilder : public GeneoCoarseSpaceDenseBuilderPython { + + public: + /* Inherit the constructors */ + using GeneoCoarseSpaceDenseBuilderPython::GeneoCoarseSpaceDenseBuilderPython; + + /* Trampoline (need one for each virtual function) */ + void compute_coarse_space(py::array_t Ai, py::array_t Bi) override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + GeneoCoarseSpaceDenseBuilderPython, /* Parent class */ + compute_coarse_space, /* Name of function in C++ (must match Python name) */ + Ai, + Bi /* Argument(s) */ + ); + } +}; + +template +void declare_geneo_coarse_space_dense_builder(py::module &m, const std::string &className) { + + using Class = GeneoCoarseSpaceDenseBuilder; + py::class_> py_class(m, className.c_str()); + py_class.def(py::init([](int size_wo_overlap, py::array_t Ai, py::array_t Bi, char symmetry, char uplo, int geneo_nu) { + Matrix Ai_mat; + Ai_mat.assign(Ai.shape()[0], Ai.shape()[1], Ai.mutable_data(), false); + Matrix Bi_mat; + Bi_mat.assign(Bi.shape()[0], Bi.shape()[1], Bi.mutable_data(), false); + return Class::GeneoWithNu(size_wo_overlap, Ai_mat, Bi_mat, symmetry, uplo, geneo_nu); + }), + py::arg("size_wo_overlap"), + py::arg("Ai"), + py::arg("Bi"), + py::arg("symmetry"), + py::arg("uplo"), + py::kw_only(), + py::arg("geneo_nu")); + py_class.def(py::init([](int size_wo_overlap, py::array_t Ai, py::array_t Bi, char symmetry, char uplo, double geneo_threshold) { + Matrix Ai_mat; + Ai_mat.assign(Ai.shape()[0], Ai.shape()[1], Ai.mutable_data(), false); + Matrix Bi_mat; + Bi_mat.assign(Bi.shape()[0], Bi.shape()[1], Bi.mutable_data(), false); + return Class::GeneoWithThreshold(size_wo_overlap, Ai_mat, Bi_mat, symmetry, uplo, geneo_threshold); + }), + py::arg("size_wo_overlap"), + py::arg("Ai"), + py::arg("Bi"), + py::arg("symmetry"), + py::arg("uplo"), + py::kw_only(), + py::arg("geneo_threshold")); +} + +template +void declare_virtual_geneo_coarse_space_dense_builder(py::module &m, const std::string &className) { + + using Class = GeneoCoarseSpaceDenseBuilderPython; + py::class_, VirtualCoarseSpaceBuilder> py_class(m, className.c_str()); + py_class.def(py::init([](int size_wo_overlap, py::array_t Ai, py::array_t Bi, char symmetry, char uplo, int geneo_nu) { + Matrix Ai_mat; + Ai_mat.assign(Ai.shape()[0], Ai.shape()[1], Ai.mutable_data(), false); + Matrix Bi_mat; + Bi_mat.assign(Bi.shape()[0], Bi.shape()[1], Bi.mutable_data(), false); + return PyGeneoCoarseSpaceDenseBuilder(size_wo_overlap, Ai_mat, Bi_mat, symmetry, uplo, geneo_nu, -1); + }), + py::arg("size_wo_overlap"), + py::arg("Ai"), + py::arg("Bi"), + py::arg("symmetry"), + py::arg("uplo"), + py::kw_only(), + py::arg("geneo_nu")); + py_class.def(py::init([](int size_wo_overlap, py::array_t Ai, py::array_t Bi, char symmetry, char uplo, double geneo_threshold) { + Matrix Ai_mat; + Ai_mat.assign(Ai.shape()[0], Ai.shape()[1], Ai.mutable_data(), false); + Matrix Bi_mat; + Bi_mat.assign(Bi.shape()[0], Bi.shape()[1], Bi.mutable_data(), false); + return PyGeneoCoarseSpaceDenseBuilder(size_wo_overlap, Ai_mat, Bi_mat, symmetry, uplo, 0, geneo_threshold); + }), + py::arg("size_wo_overlap"), + py::arg("Ai"), + py::arg("Bi"), + py::arg("symmetry"), + py::arg("uplo"), + py::kw_only(), + py::arg("geneo_threshold")); + py_class.def("set_coarse_space", &Class::set_coarse_space); + py_class.def_property_readonly("symmetry", &Class::get_symmetry); + py_class.def_property_readonly("geneo_nu", &Class::get_geneo_nu); + py_class.def_property_readonly("geneo_threshold", &Class::get_geneo_threshold); +} +#endif diff --git a/src/htool/solver/interfaces/virtual_coarse_operator_builder.hpp b/src/htool/solver/interfaces/virtual_coarse_operator_builder.hpp new file mode 100644 index 0000000..1a11382 --- /dev/null +++ b/src/htool/solver/interfaces/virtual_coarse_operator_builder.hpp @@ -0,0 +1,14 @@ +#ifndef HTOOL_PYTHON_VIRTUAL_COARSE_OPERATOR_BUILDER_HPP +#define HTOOL_PYTHON_VIRTUAL_COARSE_OPERATOR_BUILDER_HPP + +#include +#include +namespace py = pybind11; + +template +void declare_virtual_coarse_operator_builder(py::module &m, const std::string &className, const std::string &base_class_name) { + using BaseClass = VirtualCoarseOperatorBuilder; + py::class_(m, base_class_name.c_str()); +} + +#endif diff --git a/src/htool/solver/interfaces/virtual_coarse_space_builder.hpp b/src/htool/solver/interfaces/virtual_coarse_space_builder.hpp new file mode 100644 index 0000000..93d6804 --- /dev/null +++ b/src/htool/solver/interfaces/virtual_coarse_space_builder.hpp @@ -0,0 +1,61 @@ +#ifndef HTOOL_PYTHON_VIRTUAL_COARSE_SPACE_BUILDER_HPP +#define HTOOL_PYTHON_VIRTUAL_COARSE_SPACE_BUILDER_HPP + +#include +#include +#include + +namespace py = pybind11; + +template +class VirtualCoarseSpaceBuilderPython : public VirtualCoarseSpaceBuilder { + py::array_t m_coarse_space; + + public: + Matrix build_coarse_space() { + compute_coarse_space(); + if (m_coarse_space.ndim() != 2) { + htool::Logger::get_instance() + .log(LogLevel::ERROR, "Wrong dimension for coarse space matrix when building coarse space."); // LCOV_EXCL_LINE + } + + Matrix coarse_space(m_coarse_space.shape()[0], m_coarse_space.shape()[1]); + std::copy_n(m_coarse_space.data(), m_coarse_space.shape()[0] * m_coarse_space.shape()[1], coarse_space.data()); + return coarse_space; + } + + // lcov does not see it because of trampoline I assume + virtual void compute_coarse_space() const = 0; // LCOV_EXCL_LINE + + void set_coarse_space(py::array_t coarse_space) { + m_coarse_space = coarse_space; + } +}; + +template +class PyVirtualCoarseSpaceBuilder : public VirtualCoarseSpaceBuilderPython { + public: + using VirtualCoarseSpaceBuilderPython::VirtualCoarseSpaceBuilderPython; + + /* Trampoline (need one for each virtual function) */ + virtual void compute_coarse_space() const override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + VirtualCoarseSpaceBuilderPython, /* Parent class */ + compute_coarse_space /* Name of function in C++ (must match Python name) */ + ); + } +}; + +template +void declare_virtual_coarse_space_builder(py::module &m, const std::string &className, const std::string &base_class_name) { + using BaseClass = VirtualCoarseSpaceBuilder; + py::class_(m, base_class_name.c_str()); + + using Class = VirtualCoarseSpaceBuilderPython; + py::class_> py_class(m, className.c_str()); + py_class.def(py::init<>()); + py_class.def("compute_coarse_space", &Class::compute_coarse_space); + py_class.def("set_coarse_space", &Class::set_coarse_space); +} +#endif diff --git a/src/htool/solver/solver.hpp b/src/htool/solver/solver.hpp index 1f6fe48..826cb0c 100644 --- a/src/htool/solver/solver.hpp +++ b/src/htool/solver/solver.hpp @@ -20,19 +20,19 @@ void declare_DDM(py::module &m, const std::string &className) { const std::vector &, const std::vector> &>()); py_class.def("facto_one_level", &Class::facto_one_level); - py_class.def("build_coarse_space", [](Class &self, py::array_t Ki) { - if (Ki.ndim() != 2) { - throw std::invalid_argument("Wrong dimension for local matrix when building coarse space\n"); // LCOV_EXCL_LINE - } - if (Ki.shape()[0] != self.get_local_size() && Ki.shape()[1] != self.get_local_size()) { - throw std::invalid_argument("Wrong size for local matrix when building coarse space: (" + std::to_string(Ki.shape()[0]) + "," + std::to_string(Ki.shape()[1]) + ") vs (" + std::to_string(self.get_local_size()) + "," + std::to_string(self.get_local_size()) + ")\n"); // LCOV_EXCL_LINE - } + // py_class.def("build_coarse_space", [](Class &self, py::array_t Ki) { + // if (Ki.ndim() != 2) { + // throw std::invalid_argument("Wrong dimension for local matrix when building coarse space\n"); // LCOV_EXCL_LINE + // } + // if (Ki.shape()[0] != self.get_local_size() && Ki.shape()[1] != self.get_local_size()) { + // throw std::invalid_argument("Wrong size for local matrix when building coarse space: (" + std::to_string(Ki.shape()[0]) + "," + std::to_string(Ki.shape()[1]) + ") vs (" + std::to_string(self.get_local_size()) + "," + std::to_string(self.get_local_size()) + ")\n"); // LCOV_EXCL_LINE + // } - Matrix Ki_mat(Ki.shape()[0], Ki.shape()[1]); - std::copy_n(Ki.data(), Ki.shape()[0] * Ki.shape()[1], Ki_mat.data()); - self.build_coarse_space(Ki_mat); - }); - // py_class.def("get_local_to_global_numbering", &Class::get_local_to_global_numbering); + // Matrix Ki_mat(Ki.shape()[0], Ki.shape()[1]); + // std::copy_n(Ki.data(), Ki.shape()[0] * Ki.shape()[1], Ki_mat.data()); + // self.build_coarse_space(Ki_mat); + // }); + py_class.def("build_coarse_space", py::overload_cast &, VirtualCoarseOperatorBuilder &>(&Class::build_coarse_space)); py_class.def( "solve", [](Class &self, py::array_t x, const py::array_t b, std::string hpddm_args) { // HPDDM arguments diff --git a/src/htool/solver/utility.hpp b/src/htool/solver/utility.hpp index 0af8e18..f32ad8d 100644 --- a/src/htool/solver/utility.hpp +++ b/src/htool/solver/utility.hpp @@ -34,10 +34,22 @@ void declare_solver_utility(py::module &m, std::string prefix = "") { "solver", [](const DefaultDDMSolverBuilderAddingOverlap &self) { return &self.solver; }, py::return_value_policy::reference_internal); default_ddm_solver_adding_overlap_class.def_property_readonly( "local_to_global_numbering", [](const DefaultDDMSolverBuilderAddingOverlap &self) { return &self.local_to_global_numbering; }, py::return_value_policy::reference_internal); + default_ddm_solver_adding_overlap_class.def_property_readonly( + "block_diagonal_dense_matrix", [](const DefaultDDMSolverBuilderAddingOverlap &self) { + py::array_t mat({self.block_diagonal_dense_matrix.nb_cols(), self.block_diagonal_dense_matrix.nb_rows()}, self.block_diagonal_dense_matrix.data(), py::capsule(self.block_diagonal_dense_matrix.data())); + return mat; + }, + py::return_value_policy::reference_internal); py::class_ default_ddm_solver_class(m, default_ddm_solver_name.c_str()); default_ddm_solver_class.def(py::init &, const HMatrix &, const std::vector &, const std::vector> &>(), py::keep_alive<1, 2>(), py::keep_alive<1, 4>(), py::keep_alive<1, 5>()); default_ddm_solver_class.def_property_readonly( "solver", [](const DefaultDDMSolverBuilder &self) { return &self.solver; }, py::return_value_policy::reference_internal); + default_ddm_solver_class.def_property_readonly( + "block_diagonal_dense_matrix", [](const DefaultDDMSolverBuilder &self) { + py::array_t mat({self.block_diagonal_dense_matrix.nb_cols(), self.block_diagonal_dense_matrix.nb_rows()}, self.block_diagonal_dense_matrix.data(), py::capsule(self.block_diagonal_dense_matrix.data())); + return mat; + }, + py::return_value_policy::reference_internal); } #endif diff --git a/tests/conftest.py b/tests/conftest.py index ea68427..8f0dafe 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,11 +2,11 @@ import pathlib import struct +import Htool import mpi4py import numpy as np import pytest -import Htool from example.define_custom_dense_blocks_generator import CustomDenseBlocksGenerator from example.define_custom_generators import CustomGenerator from example.define_custom_local_operator import CustomLocalOperator @@ -375,7 +375,7 @@ def load_data_solver(symmetry, mu): (m, n) = struct.unpack("@II", data[:8]) # print(m,n) local_neumann_matrix = np.frombuffer(data[8:], dtype=np.dtype("complex128")) - local_neumann_matrix = np.transpose(local_neumann_matrix.reshape((m, n))) + local_neumann_matrix = np.transpose(local_neumann_matrix.reshape((m, n),order='C')).copy("F") return [ A, x_ref, diff --git a/tests/test_ddm_solver.py b/tests/test_ddm_solver.py index 855d830..119773b 100644 --- a/tests/test_ddm_solver.py +++ b/tests/test_ddm_solver.py @@ -1,60 +1,162 @@ -import math -import os -import pathlib -import struct +import time -import matplotlib.pyplot as plt +import Htool import mpi4py import numpy as np import pytest from conftest import GeneratorFromMatrix, LocalGeneratorFromMatrix -from mpi4py import MPI +from scipy.linalg import eig, eigh -import Htool + +class CustomGeneoCoarseSpaceDenseBuilder( + Htool.VirtualComplexGeneoCoarseSpaceDenseBuilder +): + def compute_coarse_space(self, Ai, Bi): + coarse_space = None + + if self.symmetry == "S" or self.symmetry == "H": + if self.geneo_threshold > 0: + [_, coarse_space] = eigh( + Ai, Bi, subset_by_value=[self.geneo_threshold, np.inf] + ) + else: + n = Ai.shape[0] + [_, coarse_space] = eigh( + Ai, Bi, subset_by_index=[n - self.geneo_nu, n - 1] + ) + else: + [w, v] = eig(Ai, Bi) + if self.geneo_threshold > 0: + nb_eig = (w > self.geneo_threshold).sum() + coarse_space = v[:, 0:nb_eig] + else: + coarse_space = v[:, 0 : self.geneo_nu] + + self.set_coarse_space(coarse_space) @pytest.mark.parametrize("epsilon", [1e-6]) @pytest.mark.parametrize("eta", [10]) @pytest.mark.parametrize("tol", [1e-6]) @pytest.mark.parametrize( - "mu,symmetry,ddm_builder,hpddm_schwarz_method,hpddm_schwarz_coarse_correction", + "mu,symmetry,ddm_builder,hpddm_schwarz_method,hpddm_schwarz_coarse_correction,geneo_type", [ - (1, "N", "SolverBuilder", "none", "none"), - (1, "N", "SolverBuilder", "asm", "none"), - (1, "N", "SolverBuilder", "ras", "none"), - (1, "N", "DDMSolverBuilderAddingOverlap", "asm", "none"), - (1, "N", "DDMSolverBuilderAddingOverlap", "ras", "none"), - (1, "N", "DDMSolverBuilder", "asm", "none"), - (1, "N", "DDMSolverBuilder", "ras", "none"), - (10, "N", "SolverBuilder", "none", "none"), - (10, "N", "SolverBuilder", "asm", "none"), - (10, "N", "SolverBuilder", "ras", "none"), - (10, "N", "DDMSolverBuilderAddingOverlap", "asm", "none"), - (10, "N", "DDMSolverBuilderAddingOverlap", "ras", "none"), - (10, "N", "DDMSolverBuilder", "asm", "none"), - (10, "N", "DDMSolverBuilder", "ras", "none"), - (1, "S", "SolverBuilder", "none", "none"), - (1, "S", "SolverBuilder", "asm", "none"), - (1, "S", "SolverBuilder", "ras", "none"), - (1, "S", "DDMSolverBuilderAddingOverlap", "asm", "none"), - (1, "S", "DDMSolverBuilderAddingOverlap", "ras", "none"), - (1, "S", "DDMSolverBuilder", "asm", "none"), - (1, "S", "DDMSolverBuilder", "ras", "none"), - (10, "S", "SolverBuilder", "none", "none"), - (10, "S", "SolverBuilder", "asm", "none"), - (10, "S", "SolverBuilder", "ras", "none"), - (10, "S", "DDMSolverBuilderAddingOverlap", "asm", "none"), - (10, "S", "DDMSolverBuilderAddingOverlap", "ras", "none"), - (1, "S", "DDMSolverBuilderAddingOverlap", "asm", "additive"), - (1, "S", "DDMSolverBuilderAddingOverlap", "ras", "additive"), - (10, "S", "DDMSolverBuilderAddingOverlap", "asm", "additive"), - (10, "S", "DDMSolverBuilderAddingOverlap", "ras", "additive"), - (10, "S", "DDMSolverBuilder", "asm", "none"), - (10, "S", "DDMSolverBuilder", "ras", "none"), - (1, "S", "DDMSolverBuilder", "asm", "additive"), - (1, "S", "DDMSolverBuilder", "ras", "additive"), - (10, "S", "DDMSolverBuilder", "asm", "additive"), - (10, "S", "DDMSolverBuilder", "ras", "additive"), + (1, "N", "SolverBuilder", "none", "none", "none"), + (1, "N", "SolverBuilder", "asm", "none", "none"), + (1, "N", "SolverBuilder", "ras", "none", "none"), + (1, "N", "DDMSolverBuilderAddingOverlap", "asm", "none", "none"), + (1, "N", "DDMSolverBuilderAddingOverlap", "ras", "none", "none"), + (1, "N", "DDMSolverBuilder", "asm", "none", "none"), + (1, "N", "DDMSolverBuilder", "ras", "none", "none"), + (10, "N", "SolverBuilder", "none", "none", "none"), + (10, "N", "SolverBuilder", "asm", "none", "none"), + (10, "N", "SolverBuilder", "ras", "none", "none"), + (10, "N", "DDMSolverBuilderAddingOverlap", "asm", "none", "none"), + (10, "N", "DDMSolverBuilderAddingOverlap", "ras", "none", "none"), + (10, "N", "DDMSolverBuilder", "asm", "none", "none"), + (10, "N", "DDMSolverBuilder", "ras", "none", "none"), + (1, "S", "SolverBuilder", "none", "none", "none"), + (1, "S", "SolverBuilder", "asm", "none", "none"), + (1, "S", "SolverBuilder", "ras", "none", "none"), + (1, "S", "DDMSolverBuilderAddingOverlap", "asm", "none", "none"), + (1, "S", "DDMSolverBuilderAddingOverlap", "ras", "none", "none"), + (1, "S", "DDMSolverBuilder", "asm", "none", "none"), + (1, "S", "DDMSolverBuilder", "ras", "none", "none"), + (10, "S", "SolverBuilder", "none", "none", "none"), + (10, "S", "SolverBuilder", "asm", "none", "none"), + (10, "S", "SolverBuilder", "ras", "none", "none"), + (10, "S", "DDMSolverBuilderAddingOverlap", "asm", "none", "none"), + (10, "S", "DDMSolverBuilderAddingOverlap", "ras", "none", "none"), + (1, "S", "DDMSolverBuilderAddingOverlap", "asm", "additive", "geneo_nu"), + (1, "S", "DDMSolverBuilderAddingOverlap", "ras", "additive", "geneo_nu"), + (10, "S", "DDMSolverBuilderAddingOverlap", "asm", "additive", "geneo_nu"), + (10, "S", "DDMSolverBuilderAddingOverlap", "ras", "additive", "geneo_nu"), + (1, "S", "DDMSolverBuilderAddingOverlap", "asm", "additive", "geneo_threshold"), + (1, "S", "DDMSolverBuilderAddingOverlap", "ras", "additive", "geneo_threshold"), + ( + 10, + "S", + "DDMSolverBuilderAddingOverlap", + "asm", + "additive", + "geneo_threshold", + ), + ( + 10, + "S", + "DDMSolverBuilderAddingOverlap", + "ras", + "additive", + "geneo_threshold", + ), + (1, "S", "DDMSolverBuilderAddingOverlap", "asm", "additive", "custom_geneo_nu"), + (1, "S", "DDMSolverBuilderAddingOverlap", "ras", "additive", "custom_geneo_nu"), + ( + 10, + "S", + "DDMSolverBuilderAddingOverlap", + "asm", + "additive", + "custom_geneo_nu", + ), + ( + 10, + "S", + "DDMSolverBuilderAddingOverlap", + "ras", + "additive", + "custom_geneo_nu", + ), + ( + 1, + "S", + "DDMSolverBuilderAddingOverlap", + "asm", + "additive", + "custom_geneo_threshold", + ), + ( + 1, + "S", + "DDMSolverBuilderAddingOverlap", + "ras", + "additive", + "custom_geneo_threshold", + ), + ( + 10, + "S", + "DDMSolverBuilderAddingOverlap", + "asm", + "additive", + "custom_geneo_threshold", + ), + ( + 10, + "S", + "DDMSolverBuilderAddingOverlap", + "ras", + "additive", + "custom_geneo_threshold", + ), + (10, "S", "DDMSolverBuilder", "asm", "none", "none"), + (10, "S", "DDMSolverBuilder", "ras", "none", "none"), + (1, "S", "DDMSolverBuilder", "asm", "additive", "geneo_nu"), + (1, "S", "DDMSolverBuilder", "ras", "additive", "geneo_nu"), + (10, "S", "DDMSolverBuilder", "asm", "additive", "geneo_nu"), + (10, "S", "DDMSolverBuilder", "ras", "additive", "geneo_nu"), + (1, "S", "DDMSolverBuilder", "asm", "additive", "geneo_threshold"), + (1, "S", "DDMSolverBuilder", "ras", "additive", "geneo_threshold"), + (10, "S", "DDMSolverBuilder", "asm", "additive", "geneo_threshold"), + (10, "S", "DDMSolverBuilder", "ras", "additive", "geneo_threshold"), + (1, "S", "DDMSolverBuilder", "asm", "additive", "custom_geneo_nu"), + (1, "S", "DDMSolverBuilder", "ras", "additive", "custom_geneo_nu"), + (10, "S", "DDMSolverBuilder", "asm", "additive", "custom_geneo_nu"), + (10, "S", "DDMSolverBuilder", "ras", "additive", "custom_geneo_nu"), + (1, "S", "DDMSolverBuilder", "asm", "additive", "custom_geneo_threshold"), + (1, "S", "DDMSolverBuilder", "ras", "additive", "custom_geneo_threshold"), + (10, "S", "DDMSolverBuilder", "asm", "additive", "custom_geneo_threshold"), + (10, "S", "DDMSolverBuilder", "ras", "additive", "custom_geneo_threshold"), ], # indirect=["setup_solver_dependencies"], ) @@ -68,6 +170,7 @@ def test_ddm_solver( tol, hpddm_schwarz_method, hpddm_schwarz_coarse_correction, + geneo_type, ): # ( # solver, @@ -104,18 +207,15 @@ def test_ddm_solver( UPLO, mpi4py.MPI.COMM_WORLD, ) - # print("Geometry", geometry) - # fig = plt.figure() - # ax = fig.add_subplot(1, 1, 1, projection="3d") - # ax.scatter(geometry[0, :], geometry[1, :], geometry[2, :], marker="o") - # plt.show() + solver = None + default_solver_builder = None if ddm_builder == "SolverBuilder": default_solver_builder = Htool.ComplexDefaultSolverBuilder( default_approximation.distributed_operator, default_approximation.block_diagonal_hmatrix, ) - solver = default_solver_builder.solver + elif ddm_builder == "DDMSolverBuilderAddingOverlap": default_solver_builder = Htool.ComplexDefaultDDMSolverBuilderAddingOverlap( default_approximation.distributed_operator, @@ -126,7 +226,7 @@ def test_ddm_solver( neighbors, intersections, ) - solver = default_solver_builder.solver + elif ddm_builder == "DDMSolverBuilder": local_numbering_builder = Htool.LocalNumberingBuilder( ovr_subdomain_to_global, @@ -163,8 +263,8 @@ def test_ddm_solver( neighbors, intersections, ) - solver = default_solver_builder.solver + solver = default_solver_builder.solver distributed_operator = default_approximation.distributed_operator # Solver @@ -192,7 +292,62 @@ def test_ddm_solver( solver.set_hpddm_args( "-hpddm_schwarz_coarse_correction " + hpddm_schwarz_coarse_correction ) - solver.build_coarse_space(local_neumann_matrix) + # solver.build_coarse_space(local_neumann_matrix) + # print(default_solver_builder.block_diagonal_dense_matrix.flags) + # print(local_neumann_matrix.flags) + # geneo_space_operator_builder = ( + # Htool.ComplexGeneoCoarseSpaceDenseBuilder.GeneoWithNu( + # cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank).get_size(), + # default_solver_builder.block_diagonal_dense_matrix, + # local_neumann_matrix, + # symmetry, + # UPLO, + # 2, + # ) + # ) + geneo_space_operator_builder = None + if geneo_type == "geneo_nu": + geneo_space_operator_builder = Htool.ComplexGeneoCoarseSpaceDenseBuilder( + cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank).get_size(), + default_solver_builder.block_diagonal_dense_matrix, + local_neumann_matrix, + symmetry, + UPLO, + geneo_nu=2, + ) + elif geneo_type == "geneo_threshold": + geneo_space_operator_builder = Htool.ComplexGeneoCoarseSpaceDenseBuilder( + cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank).get_size(), + default_solver_builder.block_diagonal_dense_matrix, + local_neumann_matrix, + symmetry, + UPLO, + geneo_threshold=100, + ) + elif geneo_type == "custom_geneo_nu": + geneo_space_operator_builder = CustomGeneoCoarseSpaceDenseBuilder( + cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank).get_size(), + default_solver_builder.block_diagonal_dense_matrix, + local_neumann_matrix, + symmetry, + UPLO, + geneo_nu=2, + ) + elif geneo_type == "custom_geneo_threshold": + geneo_space_operator_builder = Htool.ComplexGeneoCoarseSpaceDenseBuilder( + cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank).get_size(), + default_solver_builder.block_diagonal_dense_matrix, + local_neumann_matrix, + symmetry, + UPLO, + geneo_threshold=100, + ) + geneo_coarse_operator_builder = Htool.ComplexGeneoCoarseOperatorBuilder( + distributed_operator + ) + solver.build_coarse_space( + geneo_space_operator_builder, geneo_coarse_operator_builder + ) if hpddm_schwarz_method == "asm" or hpddm_schwarz_method == "ras": solver.facto_one_level() @@ -213,180 +368,8 @@ def test_ddm_solver( distributed_operator @ x - f ) / np.linalg.norm(f) solution_error = np.linalg.norm(x[:, 1] - x_ref) / np.linalg.norm(x_ref) - # error = np.linalg.norm(distributed_operator * x - f) - # if mpi4py.MPI.COMM_WORLD.rank == 0: - # print( - # iterative_solver, - # hpddm_schwarz_method, - # hpddm_schwarz_coarse_correction, - # epsilon, - # solver.get_information("Nb_it"), - # # error, - # # np.linalg.norm(f), - # # error / np.linalg.norm(f), - # # hpddm_args, - # ) - # print( - # np.linalg.norm(distributed_operator * x - f), - # np.linalg.norm(f), - # tol, - # mpi4py.MPI.COMM_WORLD.rank, - # ) + if mpi4py.MPI.COMM_WORLD.rank == 0: print(solver.get_information()) assert convergence_error < tol assert solution_error < epsilon * 10 - - # # DDM one level ASM wo overlap - # if rank == 0: - # print("ASM one level without overlap:") - # comm.Barrier() - # block_jacobi.set_hpddm_args("-hpddm_schwarz_method asm") - # block_jacobi.facto_one_level() - # block_jacobi.solve(x, f) - # block_jacobi.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # # DDM one level ASM wo overlap - # if rank == 0: - # print("RAS one level without overlap:") - # comm.Barrier() - # block_jacobi.set_hpddm_args("-hpddm_schwarz_method ras") - # block_jacobi.solve(x, f) - # block_jacobi.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # # Check infos - # if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: - # assert mpi4py.MPI.COMM_WORLD.Get_size() == int( - # block_jacobi.get_infos("Nb_subdomains") - # ) - - # # No precond with overlap - # if rank == 0: - # print("No precond with overlap:") - # DDM_solver.solve( - # x, - # f, - # hpddm_args="-hpddm_schwarz_method none -hpddm_max_it 200 -hpddm_tol " - # + str(tol), - # ) - # DDM_solver.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # # DDM one level ASM with overlap - # if rank == 0: - # print("ASM one level with overlap:") - # comm.Barrier() - # DDM_solver.set_hpddm_args("-hpddm_schwarz_method asm") - # if Symmetry == "S" and size > 1: - # # Matrix - # with open( - # os.path.join( - # os.path.dirname(__file__) - # + "/../lib/htool/data/data_test/" - # + folder - # + "/Ki_" - # + str(size) - # + "_" - # + str(rank) - # + ".bin" - # ), - # "rb", - # ) as input: - # data = input.read() - # (m, n) = struct.unpack("@II", data[:8]) - # # print(m,n) - # Ki = np.frombuffer(data[8:], dtype=np.dtype("complex128")) - # Ki = np.transpose(Ki.reshape((m, n))) - # DDM_solver.build_coarse_space(Ki) - # DDM_solver.facto_one_level() - # DDM_solver.solve(x, f) - # DDM_solver.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # # DDM one level RAS with overlap - # if rank == 0: - # print("RAS one level with overlap:") - # comm.Barrier() - # DDM_solver.set_hpddm_args("-hpddm_schwarz_method ras") - # DDM_solver.solve(x, f) - # DDM_solver.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # # DDM two level ASM with overlap - # if Symmetry == "S" and size > 1: - # if rank == 0: - # print("ASM two level with overlap:") - # comm.Barrier() - # DDM_solver.set_hpddm_args( - # "-hpddm_schwarz_method asm -hpddm_schwarz_coarse_correction additive" - # ) - # DDM_solver.solve(x, f) - # DDM_solver.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # if rank == 0: - # print("RAS two level with overlap:") - # comm.Barrier() - # DDM_solver.set_hpddm_args( - # "-hpddm_schwarz_method asm -hpddm_schwarz_coarse_correction additive" - # ) - # DDM_solver.solve(x, f) - # DDM_solver.print_infos() - # if mu == 1: - # error = np.linalg.norm(hmat * x - f) / np.linalg.norm(f) - # elif mu > 1: - # error = np.linalg.norm(hmat @ x - f) / np.linalg.norm(f) - # if rank == 0: - # print(error) - # assert error < tol - # x.fill(0) - - # # Check infos - # if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: - # assert mpi4py.MPI.COMM_WORLD.Get_size() == int( - # DDM_solver.get_infos("Nb_subdomains") - # )