update htool.hpp and examples

htool-ddm · Aug 6, 2024 · ff33b04 · ff33b04
1 parent 037871f
commit ff33b04
Show file tree

Hide file tree

Showing 5 changed files with 66 additions and 55 deletions.
diff --git a/examples/compression_comparison.cpp b/examples/compression_comparison.cpp
@@ -10,19 +10,17 @@ using namespace htool;
 class MyMatrix : public VirtualGenerator<double> {
     const vector<int> &m_target_permutation;
     const vector<int> &m_source_permutation;
-    const vector<double> &m_p1;
-    const vector<double> &m_p2;
+    const vector<double> &m_target_coordinates;
+    const vector<double> &m_source_coordinates;
     int m_space_dim;
     int m_nr;
     int m_nc;
 
   public:
     // Constructor
-    MyMatrix(int space_dim0, const vector<int> &target_permutation, const vector<int> &source_permutation, const vector<double> &p10, const vector<double> &p20) : m_target_permutation(target_permutation), m_source_permutation(source_permutation), m_p1(p10), m_p2(p20), m_space_dim(space_dim0), m_nr(target_permutation.size()), m_nc(source_permutation.size()) {}
-    double get_coef(const int &k, const int &j) const {
-        return 1. / (1 + std::sqrt(std::inner_product(m_p1.begin() + m_space_dim * k, m_p1.begin() + m_space_dim * k + m_space_dim, m_p2.begin() + m_space_dim * j, double(0), std::plus<double>(), [](double u, double v) { return (u - v) * (u - v); })));
-    }
+    MyMatrix(int spatial_dimension, const vector<int> &target_permutation, const vector<int> &source_permutation, const vector<double> &target_coordinates, const vector<double> &source_coordinates) : m_target_permutation(target_permutation), m_source_permutation(source_permutation), m_target_coordinates(target_coordinates), m_source_coordinates(source_coordinates), m_space_dim(spatial_dimension), m_nr(target_permutation.size()), m_nc(source_permutation.size()) {}
 
+    // Virtual function to overload, necessary
     void copy_submatrix(int M, int N, int row_offset, int col_offset, double *ptr) const override {
         for (int j = 0; j < M; j++) {
             for (int k = 0; k < N; k++) {
@@ -31,6 +29,11 @@ class MyMatrix : public VirtualGenerator<double> {
         }
     }
 
+    double get_coef(const int &k, const int &j) const {
+        return 1. / (1 + std::sqrt(std::inner_product(m_target_coordinates.begin() + m_space_dim * k, m_target_coordinates.begin() + m_space_dim * k + m_space_dim, m_source_coordinates.begin() + m_space_dim * j, double(0), std::plus<double>(), [](double u, double v) { return (u - v) * (u - v); })));
+    }
+
+    // Matrix vector product, not necessary
     std::vector<double> operator*(std::vector<double> a) {
         std::vector<double> result(m_nr, 0);
         for (int i = 0; i < m_nr; i++) {
@@ -41,6 +44,7 @@ class MyMatrix : public VirtualGenerator<double> {
         return result;
     }
 
+    // Frobenius norm, not necessary
     double normFrob() {
         double norm = 0;
         for (int i = 0; i < m_nr; i++) {
@@ -78,17 +82,18 @@ int main(int argc, char *argv[]) {
     int nc          = 100;
 
     // Geometry
-    vector<double> p1(3 * nr);
-    vector<double> p2(3 * nc);
-    create_sphere(nr, p1.data());
-    create_sphere(nc, p2.data(), {distance, 0, 0});
+    int spatial_dimension = 3;
+    vector<double> target_coordinates(spatial_dimension * nr);
+    vector<double> source_coordinates(spatial_dimension * nc);
+    create_sphere(nr, target_coordinates.data());
+    create_sphere(nc, source_coordinates.data(), {distance, 0, 0});
 
     // Clustering
     ClusterTreeBuilder<double> recursive_build_strategy;
-    Cluster<double> target_cluster = recursive_build_strategy.create_cluster_tree(nr, 3, p1.data(), 2, 1);
-    Cluster<double> source_cluster = recursive_build_strategy.create_cluster_tree(nc, 3, p2.data(), 2, 1);
+    Cluster<double> target_cluster = recursive_build_strategy.create_cluster_tree(nr, spatial_dimension, target_coordinates.data(), 2, 1);
+    Cluster<double> source_cluster = recursive_build_strategy.create_cluster_tree(nc, spatial_dimension, source_coordinates.data(), 2, 1);
 
-    MyMatrix A(3, target_cluster.get_permutation(), source_cluster.get_permutation(), p1, p2);
+    MyMatrix A(spatial_dimension, target_cluster.get_permutation(), source_cluster.get_permutation(), target_coordinates, source_coordinates);
     double norm_A = A.normFrob();
 
     // SVD with fixed rank

diff --git a/examples/smallest_example.cpp b/examples/smallest_example.cpp
@@ -1,37 +1,34 @@
-#include <htool/clustering/tree_builder/recursive_build.hpp>
-#include <htool/distributed_operator/utility.hpp>
-#include <htool/hmatrix/hmatrix_distributed_output.hpp>
+#include <htool/htool.hpp>
 #include <htool/testing/geometry.hpp>
 
 using namespace std;
 using namespace htool;
 
-class MyMatrix : public VirtualGenerator<double> {
-    const vector<int> &m_target_permutation;
-    const vector<int> &m_source_permutation;
-    const vector<double> &m_p1;
-    const vector<double> &m_p2;
-    int m_space_dim;
+class MyMatrix : public VirtualGeneratorInUserNumbering<double> {
+    const vector<double> &m_target_coordinates;
+    const vector<double> &m_source_coordinates;
+    int m_spatial_dimension;
     int m_nr;
     int m_nc;
 
   public:
     // Constructor
-    MyMatrix(int space_dim0, const vector<int> &target_permutation, const vector<int> &source_permutation, const vector<double> &p10, const vector<double> &p20) : m_target_permutation(target_permutation), m_source_permutation(source_permutation), m_p1(p10), m_p2(p20), m_space_dim(space_dim0), m_nr(target_permutation.size()), m_nc(source_permutation.size()) {}
-    double get_coef(const int &k, const int &j) const {
-        return 1. / (1 + std::sqrt(std::inner_product(m_p1.begin() + m_space_dim * k, m_p1.begin() + m_space_dim * k + m_space_dim, m_p2.begin() + m_space_dim * j, double(0), std::plus<double>(), [](double u, double v) { return (u - v) * (u - v); })));
-    }
+    MyMatrix(int space_dim, const vector<double> &target_coordinates, const vector<double> &source_coordinates) : m_target_coordinates(target_coordinates), m_source_coordinates(source_coordinates), m_spatial_dimension(space_dim), m_nr(target_coordinates.size() / m_spatial_dimension), m_nc(source_coordinates.size() / m_spatial_dimension) {}
 
-    // Virtual function to overload
-    void copy_submatrix(int M, int N, int row_offset, int col_offset, double *ptr) const override {
+    // Virtual function to overload, necessary
+    void copy_submatrix(int M, int N, const int *rows, const int *cols, double *ptr) const override {
         for (int j = 0; j < M; j++) {
             for (int k = 0; k < N; k++) {
-                ptr[j + M * k] = this->get_coef(m_target_permutation[j + row_offset], m_source_permutation[k + col_offset]);
+                ptr[j + M * k] = this->get_coef(rows[j], cols[k]);
             }
         }
     }
 
-    // Matrix vector product
+    double get_coef(const int &k, const int &j) const {
+        return 1. / (1 + std::sqrt(std::inner_product(m_target_coordinates.begin() + m_spatial_dimension * k, m_target_coordinates.begin() + m_spatial_dimension * k + m_spatial_dimension, m_source_coordinates.begin() + m_spatial_dimension * j, double(0), std::plus<double>(), [](double u, double v) { return (u - v) * (u - v); })));
+    }
+
+    // Matrix vector product, not necessary
     std::vector<double> operator*(std::vector<double> a) {
         std::vector<double> result(m_nr, 0);
         for (int j = 0; j < m_nr; j++) {
@@ -42,7 +39,7 @@ class MyMatrix : public VirtualGenerator<double> {
         return result;
     }
 
-    // Frobenius norm
+    // Frobenius norm, not necessary
     double norm() {
         double norm = 0;
         for (int j = 0; j < m_nr; j++) {
@@ -75,23 +72,26 @@ int main(int argc, char *argv[]) {
     std::string outputpath = argv[1];
 
     // Htool parameters
-    double epsilon = 0.001;
-    double eta     = 100;
+    const double epsilon = 0.001;
+    const double eta     = 100;
 
     // Geometry
-    int size = 5000;
-    vector<double> p(3 * size);
+    const int size              = 5000;
+    const int spatial_dimension = 3;
+    vector<double> p(spatial_dimension * size);
     create_sphere(size, p.data());
 
     // Clustering
     ClusterTreeBuilder<double> recursive_build_strategy;
-    Cluster<double> cluster = recursive_build_strategy.create_cluster_tree(size, 3, p.data(), 2, sizeWorld);
+    Cluster<double> cluster = recursive_build_strategy.create_cluster_tree(size, spatial_dimension, p.data(), 2, sizeWorld);
 
     // Generator
-    MyMatrix A(3, cluster.get_permutation(), cluster.get_permutation(), p, p);
+    MyMatrix A(spatial_dimension, p, p);
 
     // Distributed operator
-    DefaultApproximationBuilder<double, double> default_approximation_builder(A, cluster, cluster, epsilon, eta, 'S', 'U', MPI_COMM_WORLD);
+    char symmetry = 'S';
+    char UPLO     = 'U';
+    DefaultApproximationBuilder<double, double> default_approximation_builder(A, cluster, cluster, epsilon, eta, symmetry, UPLO, MPI_COMM_WORLD);
     DistributedOperator<double> &distributed_operator = default_approximation_builder.distributed_operator;
 
     // Matrix vector product

diff --git a/examples/visucluster.cpp b/examples/visucluster.cpp
@@ -1,5 +1,4 @@
-#include <htool/clustering/cluster_output.hpp>
-#include <htool/clustering/tree_builder/recursive_build.hpp>
+#include <htool/htool.hpp>
 #include <htool/testing/geometry.hpp>
 
 using namespace std;
@@ -19,13 +18,14 @@ int main(int argc, char *argv[]) {
     std::string outputname = argv[1];
 
     // Geometry
-    int size = 1000;
-    vector<double> p(3 * size);
+    const int size              = 1000;
+    const int spatial_dimension = 3;
+    vector<double> p(spatial_dimension * size);
     create_sphere(size, p.data());
 
     // Clustering
     ClusterTreeBuilder<double> recursive_build_strategy;
-    Cluster<double> cluster = recursive_build_strategy.create_cluster_tree(size, 3, p.data(), 2, 2);
+    Cluster<double> cluster = recursive_build_strategy.create_cluster_tree(size, spatial_dimension, p.data(), 2, 2);
 
     // Output
     save_clustered_geometry(cluster, 3, p.data(), outputname + "/clustering_output", {1, 2, 3});

diff --git a/include/htool/hmatrix/tree_builder/tree_builder.hpp b/include/htool/hmatrix/tree_builder/tree_builder.hpp
@@ -306,14 +306,14 @@ void HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision>::reset_root_o
 template <typename CoefficientPrecision, typename CoordinatePrecision>
 void HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision>::compute_blocks(const VirtualGenerator<CoefficientPrecision> &generator) const {
 
-#if defined(_OPENMP) && !defined(PYTHON_INTERFACE)
+#if defined(_OPENMP) && !defined(HTOOL_WITH_PYTHON_INTERFACE)
 #    pragma omp parallel
 #endif
     {
         // std::vector<HMatrixType *> local_dense_leaves{};
         // std::vector<HMatrixType *> local_low_rank_leaves{};
         int local_false_positive = 0;
-#if defined(_OPENMP) && !defined(PYTHON_INTERFACE)
+#if defined(_OPENMP) && !defined(HTOOL_WITH_PYTHON_INTERFACE)
 #    pragma omp for schedule(guided) nowait
 #endif
         for (int p = 0; p < m_admissible_tasks.size(); p++) {
@@ -330,15 +330,15 @@ void HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision>::compute_bloc
             }
         }
         if (m_dense_blocks_generator.get() == nullptr) {
-#if defined(_OPENMP) && !defined(PYTHON_INTERFACE)
+#if defined(_OPENMP) && !defined(HTOOL_WITH_PYTHON_INTERFACE)
 #    pragma omp for schedule(guided) nowait
 #endif
             for (int p = 0; p < m_dense_tasks.size(); p++) {
                 m_dense_tasks[p]->compute_dense_data(generator);
             }
             // local_dense_leaves.emplace_back(m_dense_tasks[p]);
         }
-#if defined(_OPENMP) && !defined(PYTHON_INTERFACE)
+#if defined(_OPENMP) && !defined(HTOOL_WITH_PYTHON_INTERFACE)
 #    pragma omp critical
 #endif
         {

diff --git a/include/htool/htool.hpp b/include/htool/htool.hpp
@@ -4,24 +4,30 @@
 #include "htool_version.hpp"
 #include "misc/define.hpp"
 
+#include "basic_types/vector.hpp"
 #include "clustering/clustering.hpp"
-
+#include "distributed_operator/distributed_operator.hpp"
+#include "distributed_operator/utility.hpp"
+#include "hmatrix/hmatrix.hpp"
+#include "hmatrix/hmatrix_distributed_output.hpp"
+#include "hmatrix/interfaces/virtual_generator.hpp"
+#include "hmatrix/linalg/interface.hpp"
 #include "hmatrix/lrmat/SVD.hpp"
 #include "hmatrix/lrmat/fullACA.hpp"
+#include "hmatrix/lrmat/linalg/interface.hpp"
 #include "hmatrix/lrmat/lrmat.hpp"
 #include "hmatrix/lrmat/partialACA.hpp"
 #include "hmatrix/lrmat/sympartialACA.hpp"
-
+#include "matrix/linalg/interface.hpp"
+#include "matrix/matrix.hpp"
 #include "misc/misc.hpp"
 #include "misc/user.hpp"
 
-#include "basic_types/vector.hpp"
-#include "hmatrix/hmatrix.hpp"
-#include "hmatrix/interfaces/virtual_generator.hpp"
-#include "matrix/matrix.hpp"
-
-#ifdef WITH_HPDDM
+#ifdef HTOOL_WITH_HPDDM
 #    include "solvers/ddm.hpp"
+#    include "solvers/geneo/coarse_operator_builder.hpp"
+#    include "solvers/geneo/coarse_space_builder.hpp"
+#    include "solvers/utility.hpp"
 #    include "wrappers/wrapper_hpddm.hpp"
 #endif
 #include "wrappers/wrapper_blas.hpp"