Skip to content

Commit

Permalink
add conversion in user numbering
Browse files Browse the repository at this point in the history
  • Loading branch information
@PierreMarchand20
PierreMarchand20 committed Jul 16, 2024
1 parent ee9660c commit 8de449c
Show file tree
Hide file tree
Showing 5 changed files with 211 additions and 71 deletions.
148 changes: 118 additions & 30 deletions include/htool/hmatrix/hmatrix.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ class HMatrix : public TreeNode<HMatrix<CoefficientPrecision, CoordinatePrecisio
mutable std::vector<const HMatrix *> m_leaves{};
mutable std::vector<const HMatrix *> m_leaves_for_symmetry{};
mutable char m_symmetry_type_for_leaves{'N'};
mutable char m_UPLO_for_leaves{'N'};
// std::vector<HMatrix *> m_leaves_in_diagonal_block{};

StorageType m_storage_type{StorageType::Hierarchical};
Expand All @@ -71,6 +72,7 @@ class HMatrix : public TreeNode<HMatrix<CoefficientPrecision, CoordinatePrecisio

if (m_symmetry_type_for_leaves == 'N' && current_hmatrix->get_symmetry() != 'N') {
m_symmetry_type_for_leaves = current_hmatrix->get_symmetry();
m_UPLO_for_leaves = current_hmatrix->get_UPLO();
}

for (const auto &child : current_hmatrix->get_children()) {
Expand Down Expand Up @@ -251,6 +253,7 @@ class HMatrix : public TreeNode<HMatrix<CoefficientPrecision, CoordinatePrecisio

// HMatrix Tree setters
char get_symmetry_for_leaves() const { return m_symmetry_type_for_leaves; }
char get_UPLO_for_leaves() const { return m_UPLO_for_leaves; }

// Infos
// const std::map<std::string, std::string> &get_infos() const { return infos; }
Expand Down Expand Up @@ -1992,37 +1995,85 @@ void copy_to_dense(const HMatrix<CoefficientPrecision, CoordinatePrecision> &hma
}
}
}
}

template <typename CoefficientPrecision, typename CoordinatePrecision>
void copy_to_dense_in_user_numbering(const HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix, CoefficientPrecision *ptr) {

const auto &target_cluster = hmatrix.get_target_cluster();
const auto &source_cluster = hmatrix.get_source_cluster();
if (!target_cluster.is_root() && !is_cluster_on_partition(target_cluster)) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Target cluster is neither root nor local, permutation is not stable and copy_to_dense_in_user_numbering cannot be used."); // LCOV_EXCL_LINE
}

if (!source_cluster.is_root() && !is_cluster_on_partition(source_cluster)) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Source cluster is neither root nor local, permutation is not stable and copy_to_dense_in_user_numbering cannot be used."); // LCOV_EXCL_LINE
}

if (is_cluster_on_partition(target_cluster) && !(target_cluster.is_permutation_local())) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Target cluster is local, but permutation is not local, copy_to_dense_in_user_numbering cannot be used"); // LCOV_EXCL_LINE}
}

if (is_cluster_on_partition(source_cluster) && !(source_cluster.is_permutation_local())) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Source cluster is local, but permutation is not local, copy_to_dense_in_user_numbering cannot be used"); // LCOV_EXCL_LINE}
}
int target_offset = target_cluster.get_offset();
int source_offset = source_cluster.get_offset();
int target_size = target_cluster.get_size();
const auto &target_permutation = target_cluster.get_permutation();
const auto &source_permutation = source_cluster.get_permutation();

for (auto leaf : hmatrix.get_leaves()) {
int local_nr = leaf->get_target_cluster().get_size();
int local_nc = leaf->get_source_cluster().get_size();
int row_offset = leaf->get_target_cluster().get_offset();
int col_offset = leaf->get_source_cluster().get_offset();

if (leaf->is_dense()) {
for (int k = 0; k < local_nc; k++) {
for (int j = 0; j < local_nr; j++) {
ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] = (*leaf->get_dense_data())(j, k);
}
}
} else {

// char symmetry_type = hmatrix.get_symmetry_for_leaves();
// for (auto leaf : hmatrix.get_leaves_for_symmetry()) {
// int local_nr = leaf->get_target_cluster().get_size();
// int local_nc = leaf->get_source_cluster().get_size();
// int col_offset = leaf->get_target_cluster().get_offset() - source_offset;
// int row_offset = leaf->get_source_cluster().get_offset() - target_offset;
// if (leaf->is_dense()) {
// for (int j = 0; j < local_nr; j++) {
// for (int k = 0; k < local_nc; k++) {
// ptr[k + row_offset + (j + col_offset) * target_size] = (*leaf->get_dense_data())(j, k);
// }
// }
// } else {

// Matrix<CoefficientPrecision> low_rank_to_dense(local_nr, local_nc);
// leaf->get_low_rank_data()->copy_to_dense(low_rank_to_dense.data());
// for (int j = 0; j < local_nr; j++) {
// for (int k = 0; k < local_nc; k++) {
// ptr[k + row_offset + (j + col_offset) * target_size] = low_rank_to_dense(j, k);
// }
// }
// }
// if (symmetry_type == 'H') {
// for (int k = 0; k < local_nc; k++) {
// for (int j = 0; j < local_nr; j++) {
// ptr[k + row_offset + (j + col_offset) * target_size] = conj_if_complex(ptr[k + row_offset + (j + col_offset) * target_size]);
// }
// }
// }
// }
Matrix<CoefficientPrecision> low_rank_to_dense(local_nr, local_nc);
leaf->get_low_rank_data()->copy_to_dense(low_rank_to_dense.data());
for (int k = 0; k < local_nc; k++) {
for (int j = 0; j < local_nr; j++) {
ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] = low_rank_to_dense(j, k);
}
}
}
}
if (hmatrix.get_symmetry_for_leaves() != 'N') {
for (auto leaf : hmatrix.get_leaves_for_symmetry()) {
int local_nr = leaf->get_target_cluster().get_size();
int local_nc = leaf->get_source_cluster().get_size();
int row_offset = leaf->get_target_cluster().get_offset();
int col_offset = leaf->get_source_cluster().get_offset();

if (leaf->get_target_cluster().get_offset() == leaf->get_source_cluster().get_offset() and hmatrix.get_UPLO_for_leaves() == 'L') {
for (int k = 0; k < local_nc; k++) {
for (int j = k + 1; j < local_nr; j++) {
ptr[source_permutation[k + col_offset] - target_offset + (target_permutation[j + row_offset] - source_offset) * target_size] = hmatrix.get_symmetry_for_leaves() == 'S' ? ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] : conj_if_complex(ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size]);
}
}
} else if (leaf->get_target_cluster().get_offset() == leaf->get_source_cluster().get_offset() and hmatrix.get_UPLO_for_leaves() == 'U') {
for (int k = 0; k < local_nc; k++) {
for (int j = 0; j < k; j++) {
ptr[source_permutation[k + col_offset] - target_offset + (target_permutation[j + row_offset] - source_offset) * target_size] = hmatrix.get_symmetry_for_leaves() == 'S' ? ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] : conj_if_complex(ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size]);
}
}
} else {
for (int k = 0; k < local_nc; k++) {
for (int j = 0; j < local_nr; j++) {
ptr[source_permutation[k + col_offset] - target_offset + (target_permutation[j + row_offset] - source_offset) * target_size] = hmatrix.get_symmetry_for_leaves() == 'S' ? ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] : conj_if_complex(ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size]);
}
}
}
}
}
}

template <typename CoefficientPrecision, typename CoordinatePrecision>
Expand Down Expand Up @@ -2054,6 +2105,43 @@ void copy_diagonal(const HMatrix<CoefficientPrecision, CoordinatePrecision> &hma
}
}

template <typename CoefficientPrecision, typename CoordinatePrecision>
void copy_diagonal_in_user_numbering(const HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix, CoefficientPrecision *ptr) {
if (hmatrix.get_target_cluster().get_offset() != hmatrix.get_source_cluster().get_offset() || hmatrix.get_target_cluster().get_size() != hmatrix.get_source_cluster().get_size()) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Matrix is not square a priori, get_local_diagonal cannot be used"); // LCOV_EXCL_LINE
}

if (!(hmatrix.get_target_cluster().is_root() && hmatrix.get_source_cluster().is_root()) && !(is_cluster_on_partition(hmatrix.get_target_cluster()) && is_cluster_on_partition(hmatrix.get_source_cluster()))) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Clusters are neither root nor local, permutations are not stable and copy_diagonal_in_user_numbering cannot be used."); // LCOV_EXCL_LINE
}

if ((is_cluster_on_partition(hmatrix.get_target_cluster()) && is_cluster_on_partition(hmatrix.get_source_cluster())) && !(hmatrix.get_target_cluster().is_permutation_local() && hmatrix.get_source_cluster().is_permutation_local())) {
htool::Logger::get_instance().log(LogLevel::ERROR, "Clusters are local, but permutations are not local, copy_diagonal_in_user_numbering cannot be used"); // LCOV_EXCL_LINE}
}

const auto &permutation = hmatrix.get_target_cluster().get_permutation();
int target_offset = hmatrix.get_target_cluster().get_offset();
// int source_offset = hmatrix.get_source_cluster().get_offset();

for (auto leaf : hmatrix.get_leaves()) {
int local_nr = leaf->get_target_cluster().get_size();
int local_nc = leaf->get_source_cluster().get_size();
int offset_i = leaf->get_target_cluster().get_offset();
int offset_j = leaf->get_source_cluster().get_offset();
if (leaf->is_dense() && offset_i == offset_j) {
for (int k = 0; k < std::min(local_nr, local_nc); k++) {
ptr[permutation[k + offset_i] - target_offset] = (*leaf->get_dense_data())(k, k);
}
} else if (leaf->is_low_rank() && offset_i == offset_j) { // pretty rare...
Matrix<CoefficientPrecision> low_rank_to_dense(local_nc, local_nr);
leaf->get_low_rank_data()->copy_to_dense(low_rank_to_dense.data());
for (int k = 0; k < std::min(local_nr, local_nc); k++) {
ptr[permutation[k + offset_i] - target_offset] = (low_rank_to_dense)(k, k);
}
}
}
}

// template <typename CoefficientPrecision>
// Matrix<CoefficientPrecision> HMatrix<CoefficientPrecision,CoordinatePrecision>::get_local_interaction(bool permutation) const {
// int local_size_source = cluster_tree_s->get_masteroffset(rankWorld).second;
Expand Down
7 changes: 4 additions & 3 deletions include/htool/misc/petsc.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,14 @@
#include "define.hpp"
#include "misc.hpp"

#include "../clustering/bounding_box_1.hpp"
#include "../clustering/pca.hpp"
#include "../clustering/clustering.hpp"
#include "../distributed_operator/distributed_operator.hpp"
#include "../distributed_operator/implementations/partition_from_cluster.hpp"
#include "../hmatrix/hmatrix.hpp"
#include "../hmatrix/lrmat/SVD.hpp"
#include "../hmatrix/lrmat/fullACA.hpp"
#include "../hmatrix/lrmat/lrmat.hpp"
#include "../hmatrix/lrmat/sympartialACA.hpp"
#include "../types/hmatrix.hpp"

#if defined(__clang__)
# pragma clang diagnostic pop
Expand Down
40 changes: 19 additions & 21 deletions tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_complex_double.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,27 +6,25 @@ using namespace htool;
int main(int argc, char *argv[]) {
MPI_Init(&argc, &argv);

bool is_error = false;
const int number_of_rows = 200;
const int number_of_rows_increased = 400;
const int number_of_columns = 200;
const int number_of_columns_increased = 400;

for (auto use_local_cluster : {true, false}) {
for (auto epsilon : {1e-14, 1e-6}) {
for (auto use_dense_block_generator : {true, false}) {
// Square matrix
is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);

is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'S', 'U', epsilon, use_dense_block_generator);
is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'S', 'L', epsilon, use_dense_block_generator);

is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexHermitian>(number_of_rows, number_of_columns, use_local_cluster, 'H', 'U', epsilon, use_dense_block_generator);
is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexHermitian>(number_of_rows, number_of_columns, use_local_cluster, 'H', 'L', epsilon, use_dense_block_generator);

// Rectangle matrix
is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplex>(number_of_rows_increased, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);
is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplex>(number_of_rows, number_of_columns_increased, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);
bool is_error = false;

for (auto nr : {200, 400}) {
for (auto nc : {200, 400}) {
for (auto use_local_cluster : {true, false}) {
for (auto epsilon : {1e-14, 1e-6}) {
for (auto use_dense_block_generator : {true, false}) {
std::cout << nr << " " << nc << " " << use_local_cluster << " " << epsilon << " " << use_dense_block_generator << "\n";

is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexSymmetric>(nr, nc, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);
if (nr == nc) {
for (auto UPLO : {'U', 'L'}) {
is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexSymmetric>(nr, nr, use_local_cluster, 'S', UPLO, epsilon, use_dense_block_generator);

is_error = is_error || test_hmatrix_build<std::complex<double>, GeneratorTestComplexHermitian>(nr, nr, use_local_cluster, 'H', UPLO, epsilon, use_dense_block_generator);
}
}
}
}
}
}
}
Expand Down
36 changes: 19 additions & 17 deletions tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_double.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,23 +6,25 @@ using namespace htool;
int main(int argc, char *argv[]) {
MPI_Init(&argc, &argv);

bool is_error = false;
const int number_of_rows = 200;
const int number_of_rows_increased = 400;
const int number_of_columns = 200;
const int number_of_columns_increased = 400;

for (auto use_local_cluster : {true, false}) {
for (auto epsilon : {1e-14, 1e-6}) {
for (auto use_dense_block_generator : {true, false}) {
// Square matrix
is_error = is_error || test_hmatrix_build<double, GeneratorTestDoubleSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);
is_error = is_error || test_hmatrix_build<double, GeneratorTestDoubleSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'S', 'U', epsilon, use_dense_block_generator);
is_error = is_error || test_hmatrix_build<double, GeneratorTestDoubleSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'S', 'L', epsilon, use_dense_block_generator);

// Rectangle matrix
is_error = is_error || test_hmatrix_build<double, GeneratorTestDouble>(number_of_rows_increased, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);
is_error = is_error || test_hmatrix_build<double, GeneratorTestDouble>(number_of_rows, number_of_columns_increased, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);
bool is_error = false;

for (auto nr : {200, 400}) {
for (auto nc : {200, 400}) {
for (auto use_local_cluster : {true, false}) {
for (auto epsilon : {1e-14, 1e-6}) {
for (auto use_dense_block_generator : {true, false}) {
std::cout << nr << " " << nc << " " << use_local_cluster << " " << epsilon << " " << use_dense_block_generator << "\n";

is_error = is_error || test_hmatrix_build<double, GeneratorTestDoubleSymmetric>(nr, nc, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator);

if (nr == nc) {
for (auto UPLO : {'U', 'L'}) {
std::cout << UPLO << "\n";
is_error = is_error || test_hmatrix_build<double, GeneratorTestDoubleSymmetric>(nr, nc, use_local_cluster, 'S', UPLO, epsilon, use_dense_block_generator);
}
}
}
}
}
}
}
Expand Down
Loading

0 comments on commit 8de449c

Please sign in to comment.