diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7d84eeb..56c3bce 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -29,27 +29,27 @@ jobs: - name: "Ubuntu OpenMPI g++" CC_COMPILER: gcc CXX_COMPILER: g++ - container: "ubuntu_openmpi" + container: "ubuntu_gcc_openmpi" MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe" USE_SANITIZER: "" CMAKE_BUILD_TYPE: Release - DOCKER_OPTIONS: " " + DOCKER_OPTIONS: "--cap-add SYS_PTRACE" CODE_COVERAGE: "OFF" - name: "Ubuntu OpenMPI clang++" CC_COMPILER: clang CXX_COMPILER: clang++ - container: "ubuntu_openmpi" + container: "ubuntu_clang_openmpi" MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe" USE_SANITIZER: "" CMAKE_BUILD_TYPE: Release - DOCKER_OPTIONS: " " + DOCKER_OPTIONS: "--cap-add SYS_PTRACE" CODE_COVERAGE: "OFF" - name: "Ubuntu MPICH g++" CC_COMPILER: gcc CXX_COMPILER: g++ - container: "ubuntu_mpich" + container: "ubuntu_gcc_mpich" MPIEXEC_PREFLAGS: "" CMAKE_BUILD_TYPE: Debug DOCKER_OPTIONS: " " @@ -58,35 +58,37 @@ jobs: - name: "Ubuntu MPICH clang++" CC_COMPILER: clang CXX_COMPILER: clang++ - container: "ubuntu_mpich" + container: "ubuntu_clang_mpich" MPIEXEC_PREFLAGS: "" CMAKE_BUILD_TYPE: Release DOCKER_OPTIONS: " " CODE_COVERAGE: "OFF" # Hangs on github - # - name: "Debian OpenMPI g++" - # CC_COMPILER: gcc - # CXX_COMPILER: g++ - # container: "debian_openmpi" - # MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe --mca btl_vader_single_copy_mechanism none" - # USE_SANITIZER: "" - # CMAKE_BUILD_TYPE: Debug - # DOCKER_OPTIONS: "--cap-add SYS_PTRACE" - - # - name: "Debian OpenMPI clang++" - # CC_COMPILER: clang - # CXX_COMPILER: clang++ - # container: "debian_openmpi" - # MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe --mca btl_vader_single_copy_mechanism none" - # USE_SANITIZER: "" - # CMAKE_BUILD_TYPE: Debug - # DOCKER_OPTIONS: "--cap-add SYS_PTRACE" + - name: "Debian OpenMPI g++" + CC_COMPILER: gcc + CXX_COMPILER: g++ + container: "debian_gcc_openmpi" + MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe" + USE_SANITIZER: "" + CMAKE_BUILD_TYPE: Release + DOCKER_OPTIONS: "--cap-add SYS_PTRACE" + CODE_COVERAGE: "OFF" + + - name: "Debian OpenMPI clang++" + CC_COMPILER: clang + CXX_COMPILER: clang++ + container: "debian_clang_openmpi" + MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe" + USE_SANITIZER: "" + CMAKE_BUILD_TYPE: Release + DOCKER_OPTIONS: "--cap-add SYS_PTRACE" + CODE_COVERAGE: "OFF" - name: "Debian MPICH g++" CC_COMPILER: gcc CXX_COMPILER: g++ - container: "debian_mpich" + container: "debian_gcc_mpich" MPIEXEC_PREFLAGS: "" CMAKE_BUILD_TYPE: Release DOCKER_OPTIONS: " " @@ -95,7 +97,7 @@ jobs: - name: "Debian MPICH clang++" CC_COMPILER: clang CXX_COMPILER: clang++ - container: "debian_mpich" + container: "debian_clang_mpich" MPIEXEC_PREFLAGS: "" CMAKE_BUILD_TYPE: Release DOCKER_OPTIONS: " " @@ -115,46 +117,75 @@ jobs: if: "!contains(github.event.head_commit.message, '[ci skip]')" steps: - name: Checkout htool-python - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: submodules: "true" + - name: Activate virtualenv + run: | + . /usr/local/share/venv/bin/activate + echo PATH=$PATH >> $GITHUB_ENV + - name: Build package in release if: matrix.CODE_COVERAGE == 'OFF' run: | - pip3 install . --user + python3 -m pip install .[dev] - name: Build package in debug for coverage if: matrix.CODE_COVERAGE == 'ON' run: | mkdir build && cd build - pip3 install numpy scipy matplotlib + python3 -m pip install numpy matplotlib ruff pytest mpi4py scipy CC=${{ matrix.CC_COMPILER }} CXX=${{ matrix.CXX_COMPILER }} cmake -DUSE_SANITIZER=${{ matrix.USE_SANITIZER }} -DCMAKE_BUILD_TYPE=${{ matrix.CMAKE_BUILD_TYPE }} -DCODE_COVERAGE=${{ matrix.CODE_COVERAGE }} ../ make - name: Run tests run: | export PYTHONPATH="$PWD/build" - mpirun -np 2 ${{ matrix.MPIEXEC_PREFLAGS }} pytest tests - mpirun -np 3 ${{ matrix.MPIEXEC_PREFLAGS }} pytest tests - mpirun -np 4 ${{ matrix.MPIEXEC_PREFLAGS }} pytest tests + mpirun -np 1 ${{ matrix.MPIEXEC_PREFLAGS }} python3 -m pytest tests + mpirun -np 2 ${{ matrix.MPIEXEC_PREFLAGS }} python3 -m pytest tests + mpirun -np 3 ${{ matrix.MPIEXEC_PREFLAGS }} python3 -m pytest tests + mpirun -np 4 ${{ matrix.MPIEXEC_PREFLAGS }} python3 -m pytest tests - name: Check format if: matrix.CODE_COVERAGE == 'ON' + # https://github.com/actions/checkout/pull/762 run: | + git config --global --add safe.directory `pwd` cd build make format make cmake-format git diff --exit-code + - name: Check with ruff + run: | + ruff check example/ tests/ + - name: Generate coverage reports if: matrix.CODE_COVERAGE == 'ON' run: | - lcov --capture --base-directory ./ --directory build/ --output-file coverage.info - lcov --remove coverage.info '/usr/*' '*/hpddm/*' '*/pybind11/*' '*/lib/htool/*' --output-file coverage.info + lcov --ignore-errors mismatch --capture --base-directory ./ --directory build/ --output-file coverage.info + lcov --remove coverage.info '/usr/*' '*/hpddm/*' '*/pybind11/*' '*/lib/htool/*' --output-file ../coverage.info - - name: Upload coverage to Codecov + - uses: actions/upload-artifact@v4 + with: + path: coverage.info + if-no-files-found: error if: matrix.CODE_COVERAGE == 'ON' - uses: codecov/codecov-action@v2 + + coverage: + runs-on: ubuntu-latest + if: ${{ success() }} + needs: [ linux ] + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + - uses: actions/download-artifact@v4 + - name: Upload coverage report + uses: codecov/codecov-action@v4 with: - file: coverage.info + fail_ci_if_error: true + file: ./coverage.info + token: ${{ secrets.CODECOV_TOKEN }} + verbose: true diff --git a/.gitignore b/.gitignore index 2cee99d..d963fdc 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,37 @@ build* -__pycache__ \ No newline at end of file + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ diff --git a/.gitmodules b/.gitmodules index dc67137..d57ce97 100644 --- a/.gitmodules +++ b/.gitmodules @@ -7,3 +7,6 @@ [submodule "lib/pybind11"] path = lib/pybind11 url = https://github.com/pybind/pybind11.git +[submodule "lib/htool_generate_data_test"] + path = lib/htool_generate_data_test + url = https://github.com/PierreMarchand20/htool_generate_data_test.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 41d53ce..e8419d0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,13 +13,13 @@ else() LANGUAGES CXX) endif() -# To force c++11 +# To force c++14 if(${CMAKE_VERSION} VERSION_LESS 3.1) - add_compile_options(-std=c++11) + add_compile_options(-std=c++14) elseif(${CMAKE_VERSION} VERSION_LESS 3.6.3 AND ${CMAKE_CXX_COMPILER_ID} STREQUAL "Intel") - add_compile_options(-std=c++11) + add_compile_options(-std=c++14) else() - set(CMAKE_CXX_STANDARD 11) + set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) endif() @@ -78,11 +78,20 @@ target_include_directories( ${MPI4PY_INCLUDE_DIR}) target_link_libraries(Htool PRIVATE ${MPI_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${ARPACK_LIBRARIES} ${OpenMP_CXX_LIBRARIES}) -target_compile_definitions(Htool PRIVATE "-DPYTHON_INTERFACE" "-DWITH_HPDDM") +if("${BLA_VENDOR}" STREQUAL "Intel10_32" + OR "${BLA_VENDOR}" STREQUAL "Intel10_64lp" + OR "${BLA_VENDOR}" STREQUAL "Intel10_64lp_seq" + OR "${BLA_VENDOR}" STREQUAL "Intel10_64ilp" + OR "${BLA_VENDOR}" STREQUAL "Intel10_64ilp_seq" + OR "${BLA_VENDOR}" STREQUAL "Intel10_64_dyn") + target_compile_definitions(htool PRIVATE "-DHPDDM_MKL -DHTOOL_MKL") +endif() + +target_compile_definitions(Htool PRIVATE "-DHTOOL_WITH_PYTHON_INTERFACE" "-DHTOOL_WITH_HPDDM") if(CODE_COVERAGE AND (CMAKE_C_COMPILER_ID MATCHES "GNU" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU")) target_compile_options(Htool PRIVATE -fprofile-arcs -ftest-coverage) target_link_libraries(Htool PRIVATE gcov) endif() -# target_compile_features(Htool INTERFACE cxx_std_11) +# target_compile_features(Htool INTERFACE cxx_std_14) diff --git a/cmake_modules/FindMPI4PY.cmake b/cmake_modules/FindMPI4PY.cmake index 64f3703..0a88cbc 100644 --- a/cmake_modules/FindMPI4PY.cmake +++ b/cmake_modules/FindMPI4PY.cmake @@ -6,29 +6,32 @@ # https://compacc.fnal.gov/projects/repositories/entry/synergia2/CMake/FindMPI4PY.cmake?rev=c147eafb60728606af4fe7b1b161a660df142e9a -if(NOT MPI4PY_INCLUDE_DIR) - execute_process(COMMAND - "python3" "-c" "import mpi4py; print(mpi4py.get_include())" - OUTPUT_VARIABLE MPI4PY_INCLUDE_DIR - RESULT_VARIABLE MPI4PY_COMMAND_RESULT - OUTPUT_STRIP_TRAILING_WHITESPACE) - if(MPI4PY_COMMAND_RESULT) - message("jfa: mpi4py not found") +if(NOT PYTHON_EXECUTABLE) + set(PYTHON_EXECUTABLE "python3") +endif(NOT PYTHON_EXECUTABLE) +execute_process(COMMAND "which" "${PYTHON_EXECUTABLE}") +execute_process( + COMMAND "${PYTHON_EXECUTABLE}" "-c" "import mpi4py; print(mpi4py.get_include())" + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} + OUTPUT_VARIABLE MPI4PY_INCLUDE_DIR + RESULT_VARIABLE MPI4PY_COMMAND_RESULT + OUTPUT_STRIP_TRAILING_WHITESPACE) +if(MPI4PY_COMMAND_RESULT) + message("jfa: mpi4py not found") + set(MPI4PY_FOUND FALSE) +else(MPI4PY_COMMAND_RESULT) + if(MPI4PY_INCLUDE_DIR MATCHES "Traceback") + message("jfa: mpi4py matches traceback") + ## Did not successfully include MPI4PY set(MPI4PY_FOUND FALSE) - else(MPI4PY_COMMAND_RESULT) - if (MPI4PY_INCLUDE_DIR MATCHES "Traceback") - message("jfa: mpi4py matches traceback") - ## Did not successfully include MPI4PY - set(MPI4PY_FOUND FALSE) - else (MPI4PY_INCLUDE_DIR MATCHES "Traceback") - ## successful - set(MPI4PY_FOUND TRUE) - set(MPI4PY_INCLUDE_DIR ${MPI4PY_INCLUDE_DIR} CACHE STRING "mpi4py include path") - endif (MPI4PY_INCLUDE_DIR MATCHES "Traceback") - endif(MPI4PY_COMMAND_RESULT) -else(NOT MPI4PY_INCLUDE_DIR) - set(MPI4PY_FOUND TRUE) -endif(NOT MPI4PY_INCLUDE_DIR) + else(MPI4PY_INCLUDE_DIR MATCHES "Traceback") + ## successful + set(MPI4PY_FOUND TRUE) + set(MPI4PY_INCLUDE_DIR + ${MPI4PY_INCLUDE_DIR} + CACHE STRING "mpi4py include path") + endif(MPI4PY_INCLUDE_DIR MATCHES "Traceback") +endif(MPI4PY_COMMAND_RESULT) include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(MPI4PY DEFAULT_MSG MPI4PY_INCLUDE_DIR) \ No newline at end of file +find_package_handle_standard_args(MPI4PY DEFAULT_MSG MPI4PY_INCLUDE_DIR) diff --git a/example/cluster.py b/example/cluster.py deleted file mode 100644 index 42a5b47..0000000 --- a/example/cluster.py +++ /dev/null @@ -1,61 +0,0 @@ -import Htool -import numpy as np -import mpi4py - -# Parameters -minclustersize = 5 -Number_points = 500 -np.random.seed(0) -points_3D = np.zeros((3, Number_points)) -points_2D = np.zeros((2, Number_points)) - -# Random 3D geometry -points_3D[0, :] = np.random.random(Number_points) -points_3D[1, :] = np.random.random(Number_points) -points_3D[2, :] = np.random.random(Number_points) - - -# Cluster 3D -cluster = Htool.BoundingBox1RegularClustering(3) -cluster.set_minclustersize(minclustersize) -cluster.build(Number_points, points_3D, 2) -cluster.display(points_3D, 1) -cluster.display(points_3D, 2) - -# Cluster 3D with given partition - -sizeworld = mpi4py.MPI.COMM_WORLD.Get_size() -local_size = int(Number_points/sizeworld) -MasterOffset = np.zeros((2, sizeworld)) -for i in range(0, sizeworld-1): - MasterOffset[0, i] = i*local_size - MasterOffset[1, i] = local_size - points_3D[0, i*local_size:(i+1)*local_size] = i -points_3D[0, (sizeworld-1)*local_size:] = sizeworld-1 -MasterOffset[0, sizeworld-1] = (sizeworld-1)*local_size -MasterOffset[1, sizeworld-1] = Number_points-(sizeworld-1)*local_size - - -cluster = Htool.BoundingBox1RegularClustering(3) -cluster.set_minclustersize(minclustersize) -cluster.build(Number_points, points_3D, MasterOffset, 2) -cluster.display(points_3D, 1) -cluster.display(points_3D, 2) - -# Random 2D geometry -points_2D = np.zeros((2, Number_points)) -# points_2D[0, :] = np.random.random(Number_points) -# points_2D[1, :] = np.random.random(Number_points) - -r = np.random.uniform(low=0, high=1, size=Number_points) # radius -theta = np.random.uniform(low=0, high=2*np.pi, size=Number_points) # angle - -points_2D[0, :] = np.sqrt(r) * np.cos(theta) -points_2D[1, :] = np.sqrt(r) * np.sin(theta) - -# Cluster 2D -cluster = Htool.PCAGeometricClustering(2) -cluster.set_minclustersize(minclustersize) -cluster.build(Number_points, points_2D, 2) -cluster.display(points_2D, 1) -cluster.display(points_2D, 2) diff --git a/example/create_geometry.py b/example/create_geometry.py new file mode 100644 index 0000000..33e27a2 --- /dev/null +++ b/example/create_geometry.py @@ -0,0 +1,91 @@ +import numpy as np + + +def create_random_points_in_disk(nb_points: int): + u = np.random.rand(nb_points) + theta = 2 * np.pi * np.random.rand(nb_points) + + x = np.sqrt(u) * np.cos(theta) + y = np.sqrt(u) * np.sin(theta) + return np.array([x, y]) + + +def create_random_points_in_sphere(nb_points: int): + u = np.random.rand(nb_points) + theta = 2 * np.pi * np.random.rand(nb_points) + phi = np.arccos(2 * np.random.rand(nb_points) - 1) + + r = np.cbrt(u) + x = r * np.sin(theta) * np.cos(phi) + y = r * np.sin(theta) * np.sin(phi) + z = r * np.cos(theta) + return np.array([x, y, z]) + + +def create_random_geometries(dimension: int, nb_rows: int, nb_cols: int): + np.random.seed(0) + + if dimension == 3: + target_points = create_random_points_in_sphere(nb_rows) + source_points = create_random_points_in_sphere(nb_cols) + + if dimension == 2: + target_points = create_random_points_in_disk(nb_rows) + source_points = create_random_points_in_disk(nb_cols) + + source_points[0, :] += 2 + return [target_points, source_points] + + +def create_partitionned_geometries( + dimension: int, nb_rows: int, nb_cols: int, sizeWorld: int +): + np.random.seed(0) + target_local_size = int(nb_rows / sizeWorld) + target_partition = np.zeros((2, sizeWorld), dtype=int) + + if dimension == 2: + target_points = np.array([[], []]) + for i in range(0, sizeWorld - 1): + target_partition[0, i] = i * target_local_size + target_partition[1, i] = target_local_size + new_target_points = create_random_points_in_disk(target_partition[1, i]) + new_target_points[0, :] += 3 * i + target_points = np.concatenate((target_points, new_target_points), axis=1) + + target_partition[0, sizeWorld - 1] = (sizeWorld - 1) * target_local_size + target_partition[1, sizeWorld - 1] = ( + nb_rows - (sizeWorld - 1) * target_local_size + ) + new_target_points = create_random_points_in_disk( + target_partition[1, sizeWorld - 1] + ) + new_target_points[0, :] += 3 * (sizeWorld - 1) + target_points = np.concatenate((target_points, new_target_points), axis=1) + source_points = create_random_points_in_disk(nb_cols) + source_points[0, :] += 3 * (sizeWorld - 1) / 2 + source_points[1, :] += 3 + + if dimension == 3: + target_points = np.array([[], [], []]) + for i in range(0, sizeWorld - 1): + target_partition[0, i] = i * target_local_size + target_partition[1, i] = target_local_size + new_target_points = create_random_points_in_sphere(target_partition[1, i]) + new_target_points[0, :] += 3 * i + target_points = np.concatenate((target_points, new_target_points), axis=1) + + target_partition[0, sizeWorld - 1] = (sizeWorld - 1) * target_local_size + target_partition[1, sizeWorld - 1] = ( + nb_rows - (sizeWorld - 1) * target_local_size + ) + new_target_points = create_random_points_in_sphere( + target_partition[1, sizeWorld - 1] + ) + new_target_points[0, :] += 3 * (sizeWorld - 1) + target_points = np.concatenate((target_points, new_target_points), axis=1) + source_points = create_random_points_in_sphere(nb_cols) + source_points[0, :] += 3 * (sizeWorld - 1) / 2 + source_points[1, :] += 3 + + return [target_points, source_points, target_partition] diff --git a/example/ddm_solver.py b/example/ddm_solver.py deleted file mode 100644 index 7052b6f..0000000 --- a/example/ddm_solver.py +++ /dev/null @@ -1,88 +0,0 @@ -import Htool -import numpy as np -import mpi4py - -# MPI -comm = mpi4py.MPI.COMM_WORLD -rank = comm.Get_rank() - -# Custom generator - - -class Generator(Htool.VirtualGenerator): - - def __init__(self, points_target, points_source): - super().__init__(points_target.shape[1], points_source.shape[1]) - self.points_target = points_target - self.points_source = points_source - - def get_coef(self, i, j): - return 1.0 / (1e-5 + np.linalg.norm(self.points_target[i, :] - self.points_source[j, :])) - - def build_submatrix(self, J, K, mat): - for j in range(0, len(J)): - for k in range(0, len(K)): - mat[j, k] = 1.0 / (1.e-5 + np.linalg.norm( - self.points_target[:, J[j]] - self.points_source[:, K[k]])) - - def mat_vec(self, x): - y = np.zeros(self.nb_rows()) - for i in range(0, self.nb_rows()): - for j in range(0, self.nb_cols()): - y[i] += self.get_coef(i, j)*x[j] - return y - - def mat_mat(self, X): - Y = np.zeros((self.nb_rows(), X.shape[1])) - - for i in range(0, self.nb_rows()): - for j in range(0, X.shape[1]): - for k in range(0, self.nb_cols()): - Y[i, j] += self.get_coef(i, k)*X[k, j] - return Y - - -# Random geometry -Size = 500 -np.random.seed(0) -points_target = np.zeros((2, Size)) -points_target[0, :] = np.random.random(Size) -points_target[1, :] = np.random.random(Size) - -# Htool parameters -eta = 5 -epsilon = 1e-3 -minclustersize = 10 - -# Build H matrix -generator = Generator(points_target, points_target) -cluster = Htool.PCARegularClustering(2) -cluster.build(Size, points_target, 2) -cluster.set_minclustersize(minclustersize) -hmat = Htool.HMatrix(cluster, cluster, epsilon, eta, 'S', 'L') -hmat.build(generator, points_target) - -# Solver with block Jacobi -x_ref = np.random.random(Size) -b = hmat*x_ref -x = np.zeros(Size) -ddm_solver = Htool.DDM(hmat) - -hpddm_args = "-hpddm_compute_residual l2 " -if (rank == 0): - hpddm_args += "-hpddm_verbosity 10" -ddm_solver.set_hpddm_args(hpddm_args) -ddm_solver.facto_one_level() -ddm_solver.solve(x, b) - -# Several ways to display information -print(hmat) -hmat.print_infos() -hmat.display() -hmat.display_cluster(points_target, 2, "target") -ddm_solver.print_infos() - -if rank == 0: - nb_it = ddm_solver.get_infos("Nb_it") - print("Nb_it", nb_it) - print(np.linalg.norm(x-x_ref)/np.linalg.norm(x_ref)) diff --git a/example/define_custom_dense_blocks_generator.py b/example/define_custom_dense_blocks_generator.py new file mode 100644 index 0000000..d72cbac --- /dev/null +++ b/example/define_custom_dense_blocks_generator.py @@ -0,0 +1,14 @@ +import Htool + + +class CustomDenseBlocksGenerator(Htool.VirtualDenseBlocksGenerator): + def __init__( + self, generator, target_cluster: Htool.Cluster, source_cluster: Htool.Cluster + ): + super().__init__(target_cluster, source_cluster) + self.generator = generator + + def build_dense_blocks(self, rows_offsets, cols_offsets, blocks): + nb_blocks = len(blocks) + for i in range(nb_blocks): + self.generator.build_submatrix(rows_offsets[i], cols_offsets[i], blocks[i]) diff --git a/example/define_custom_generators.py b/example/define_custom_generators.py new file mode 100644 index 0000000..0703093 --- /dev/null +++ b/example/define_custom_generators.py @@ -0,0 +1,42 @@ +import Htool +import numpy as np + + +class CustomGenerator(Htool.VirtualGenerator): + def __init__(self, target_points, source_points): + super().__init__() + self.target_points = target_points + self.source_points = source_points + self.nb_rows = target_points.shape[1] + self.nb_cols = source_points.shape[1] + + def get_coef(self, i, j): + return 1.0 / ( + 1e-1 + np.linalg.norm(self.target_points[:, i] - self.source_points[:, j]) + ) + + def build_submatrix(self, J, K, mat): + for j in range(0, len(J)): + for k in range(0, len(K)): + mat[j, k] = 1.0 / ( + 1e-1 + + np.linalg.norm( + self.target_points[:, J[j]] - self.source_points[:, K[k]] + ) + ) + + def mat_vec(self, x): + y = np.zeros(self.nb_rows) + for i in range(0, self.nb_rows): + for j in range(0, self.nb_cols): + y[i] += self.get_coef(i, j) * x[j] + return y + + def mat_mat(self, X): + Y = np.zeros((self.nb_rows, X.shape[1])) + + for i in range(0, self.nb_rows): + for j in range(0, X.shape[1]): + for k in range(0, self.nb_cols): + Y[i, j] += self.get_coef(i, k) * X[k, j] + return Y diff --git a/example/define_custom_local_operator.py b/example/define_custom_local_operator.py new file mode 100644 index 0000000..e2bdb28 --- /dev/null +++ b/example/define_custom_local_operator.py @@ -0,0 +1,59 @@ +import Htool +import numpy as np + + +class CustomLocalOperator(Htool.LocalOperator): + def __init__( + self, + generator: Htool.VirtualGenerator, + target_cluster: Htool.Cluster, + source_cluster: Htool.Cluster, + symmetry: str = "N", + UPLO: str = "N", + target_use_permutation_to_mvprod: bool = False, + source_use_permutation_to_mvprod: bool = False, + ) -> None: + super().__init__( + target_cluster, + source_cluster, + symmetry, + UPLO, + target_use_permutation_to_mvprod, + source_use_permutation_to_mvprod, + ) + self.data = np.zeros((target_cluster.get_size(), source_cluster.get_size())) + generator.build_submatrix( + target_cluster.get_permutation()[ + target_cluster.get_offset() : target_cluster.get_offset() + + target_cluster.get_size() + ], + source_cluster.get_permutation()[ + source_cluster.get_offset() : source_cluster.get_offset() + + source_cluster.get_size() + ], + self.data, + ) + + def add_vector_product( + self, trans, alpha, input: np.array, beta, output: np.array + ) -> None: + # Beware, inplace operation needed for output to keep the underlying data + output *= beta + if trans == "N": + output += alpha * self.data.dot(input) + elif trans == "T": + output += alpha * np.transpose(self.data).dot(input) + elif trans == "C": + output += alpha * np.vdot(np.transpose(self.data), input) + + def add_matrix_product_row_major( + self, trans, alpha, input: np.array, beta, output: np.array + ) -> None: + output *= beta + if trans == "N": + output += alpha * self.data @ input + elif trans == "T": + output += alpha * np.transpose(self.data) @ input + elif trans == "C": + output += alpha * np.matrix.H(self.data) @ input + output = np.asfortranarray(output) diff --git a/example/define_custom_low_rank_generator.py b/example/define_custom_low_rank_generator.py new file mode 100644 index 0000000..a8d263b --- /dev/null +++ b/example/define_custom_low_rank_generator.py @@ -0,0 +1,31 @@ +import math + +import Htool +import numpy as np + + +class CustomSVD(Htool.VirtualLowRankGenerator): + def __init__(self, generator: Htool.VirtualGenerator): + super().__init__() + self.generator = generator + + def build_low_rank_approximation(self, rows, cols, epsilon): + submat = np.zeros((len(rows), len(cols)), order="F") + self.generator.build_submatrix(rows, cols, submat) + u, s, vh = np.linalg.svd(submat, full_matrices=False) + + norm = np.linalg.norm(submat) + svd_norm = 0 + truncated_rank = len(s) - 1 + while truncated_rank > 0 and math.sqrt(svd_norm) / norm < epsilon: + svd_norm += s[truncated_rank] ** 2 + truncated_rank -= 1 + truncated_rank += 1 + + if truncated_rank * (len(rows) + len(cols)) > (len(rows) * len(cols)): + print("coucou") + return False # the low rank approximation is not worthwhile + + self.set_U(u[:, 0:truncated_rank] * s[0:truncated_rank]) + self.set_V(vh[0:truncated_rank, :]) + return True diff --git a/example/partition_example.py b/example/partition_example.py deleted file mode 100644 index e47777b..0000000 --- a/example/partition_example.py +++ /dev/null @@ -1,105 +0,0 @@ -import Htool -import numpy as np -import mpi4py - -# Custom generator - - -class Generator(Htool.VirtualGenerator): - - def __init__(self, points_target, points_source): - super().__init__(points_target.shape[1], points_source.shape[1]) - self.points_target = points_target - self.points_source = points_source - - def get_coef(self, i, j): - return 1.0 / (1e-5 + np.linalg.norm(self.points_target[:, i] - self.points_source[:, j])) - - def build_submatrix(self, J, K, mat): - for j in range(0, len(J)): - for k in range(0, len(K)): - mat[j, k] = 1.0 / (1.e-5 + np.linalg.norm( - self.points_target[:, J[j]] - self.points_source[:, K[k]])) - - def mat_vec(self, x): - y = np.zeros(self.nb_rows()) - for i in range(0, self.nb_rows()): - for j in range(0, self.nb_cols()): - y[i] += self.get_coef(i, j)*x[j] - return y - - def mat_mat(self, X): - Y = np.zeros((self.nb_rows(), X.shape[1])) - - for i in range(0, self.nb_rows()): - for j in range(0, X.shape[1]): - for k in range(0, self.nb_cols()): - Y[i, j] += self.get_coef(i, k)*X[k, j] - return Y - - -# Htool parameters -eta = 10 -epsilon = 1e-3 -minclustersize = 10 - -# Random geometry -NbRows = 500 -NbCols = 500 -np.random.seed(0) -points_target = np.zeros((2, NbRows)) -sizeworld = mpi4py.MPI.COMM_WORLD.Get_size() -local_size = int(NbRows/sizeworld) -MasterOffset = np.zeros((2, sizeworld)) -for i in range(0, sizeworld-1): - MasterOffset[0, i] = i*local_size - MasterOffset[1, i] = local_size - points_target[0, i*local_size:(i+1)*local_size] = i - -points_target[0, (sizeworld-1)*local_size:] = sizeworld-1 -MasterOffset[0, sizeworld-1] = (sizeworld-1)*local_size -MasterOffset[1, sizeworld-1] = NbRows-(sizeworld-1)*local_size - -points_target[1, :] = np.random.random(NbRows) - -# Cluster target -cluster_target = Htool.PCARegularClustering(2) -cluster_target.set_minclustersize(minclustersize) -cluster_target.build(NbRows, points_target, MasterOffset, 2) - -if NbRows == NbCols: - points_source = points_target - cluster_source = cluster_target -else: - points_source = np.zeros((2, NbCols)) - points_source[0, :] = np.random.random(NbCols) - points_source[1, :] = np.random.random(NbCols) - cluster_source = None - - -# Build H matrix -generator = Generator(points_target, points_source) -HMatrix_test = Htool.HMatrix(cluster_target, cluster_source, epsilon, eta) -HMatrix_test.build(generator, points_target, points_source) - -# Test matrix vector product -x = np.random.rand(NbCols) -y_1 = HMatrix_test*x -y_2 = generator.mat_vec(x) - -print(np.linalg.norm(y_1-y_2)/np.linalg.norm(y_2)) - -# Test matrix matrix product -X = np.random.rand(NbCols, 2) -Y_1 = HMatrix_test @ X -Y_2 = generator.mat_mat(X) - -print(np.linalg.norm(Y_1-Y_2)/np.linalg.norm(Y_2)) - -# Several ways to display information -HMatrix_test.print_infos() -print(HMatrix_test) -print("Number_of_MPI_tasks:", HMatrix_test.get_infos("Number_of_MPI_tasks")) -HMatrix_test.display() -HMatrix_test.display_cluster(points_target, 2, "target") -HMatrix_test.display_cluster(points_source, 2, "source") diff --git a/example/smallest_example.py b/example/smallest_example.py deleted file mode 100644 index 34ca05b..0000000 --- a/example/smallest_example.py +++ /dev/null @@ -1,95 +0,0 @@ -import Htool -import numpy as np -import mpi4py - -# Custom generator - - -class Generator(Htool.VirtualGenerator): - - def __init__(self, points_target, points_source): - super().__init__(points_target.shape[1], points_source.shape[1]) - self.points_target = points_target - self.points_source = points_source - - def get_coef(self, i, j): - return 1.0 / (1e-5 + np.linalg.norm(self.points_target[:, i] - self.points_source[:, j])) - - def build_submatrix(self, J, K, mat): - for j in range(0, len(J)): - for k in range(0, len(K)): - mat[j, k] = 1.0 / (1.e-5 + np.linalg.norm( - self.points_target[:, J[j]] - self.points_source[:, K[k]])) - - def mat_vec(self, x): - y = np.zeros(self.nb_rows()) - for i in range(0, self.nb_rows()): - for j in range(0, self.nb_cols()): - y[i] += self.get_coef(i, j)*x[j] - return y - - def mat_mat(self, X): - Y = np.zeros((self.nb_rows(), X.shape[1])) - - for i in range(0, self.nb_rows()): - for j in range(0, X.shape[1]): - for k in range(0, self.nb_cols()): - Y[i, j] += self.get_coef(i, k)*X[k, j] - return Y - - -# Random geometry -NbRows = 500 -NbCols = 500 -np.random.seed(0) -points_target = np.zeros((2, NbRows)) -points_target[0, :] = np.random.random(NbRows) -points_target[1, :] = np.random.random(NbRows) - -if NbRows == NbCols: - points_source = points_target -else: - points_source = np.zeros((2, NbCols)) - points_source[0, :] = np.random.random(NbCols) - points_source[1, :] = np.random.random(NbCols) - -# Htool parameters -eta = -1000 -epsilon = 1e-1 -minclustersize = 10 - -# Build clusters -dimension = 2 -cluster_target = Htool.PCARegularClustering(dimension) -cluster_source = Htool.PCARegularClustering(dimension) -cluster_target.build(NbRows, points_target, 2) -cluster_source.build(NbCols, points_source, 2) -cluster_target.set_minclustersize(minclustersize) -cluster_source.set_minclustersize(minclustersize) - -# Build H matrix -generator = Generator(points_target, points_source) -HMatrix_test = Htool.HMatrix(cluster_target, cluster_source, epsilon, eta) -HMatrix_test.build(generator, points_target, points_source) - -# Test matrix vector product -x = np.random.rand(NbCols) -y_1 = HMatrix_test*x -y_2 = generator.mat_vec(x) - -print(np.linalg.norm(y_1-y_2)/np.linalg.norm(y_2)) - -# Test matrix matrix product -X = np.random.rand(NbCols, 2) -Y_1 = HMatrix_test @ X -Y_2 = generator.mat_mat(X) - -print(np.linalg.norm(Y_1-Y_2)/np.linalg.norm(Y_2)) - -# Several ways to display information -HMatrix_test.print_infos() -print(HMatrix_test) -print("Number_of_MPI_tasks:", HMatrix_test.get_infos("Number_of_MPI_tasks")) -HMatrix_test.display() -HMatrix_test.display_cluster(points_target, 2, "target") -HMatrix_test.display_cluster(points_source, 2, "source") diff --git a/example/solver_scipy.py b/example/solver_scipy.py deleted file mode 100644 index ca3d0e0..0000000 --- a/example/solver_scipy.py +++ /dev/null @@ -1,81 +0,0 @@ -import Htool -from mpi4py import MPI -import numpy as np -from numpy.linalg import norm -from scipy.sparse.linalg import gmres - -# Custom generator - - -class Generator(Htool.VirtualGenerator): - - def __init__(self, points): - super().__init__(points.shape[1], points.shape[1]) - self.points = points - - def get_coef(self, i, j): - return 1.0 / (1e-5 + np.linalg.norm(self.points[:, i] - self.points[:, j])) - - def build_submatrix(self, J, K, mat): - for j in range(0, len(J)): - for k in range(0, len(K)): - mat[j, k] = 1.0 / \ - (1.e-5 + - np.linalg.norm(self.points[:, J[j]] - self.points[:, K[k]])) - - def mat_vec(self, x): - y = np.zeros(self.nb_rows()) - for i in range(0, self.nb_rows()): - for j in range(0, self.nb_cols()): - y[i] += self.get_coef(i, j)*x[j] - return y - - -# SETUP -# n² points on a regular grid in a square -n = int(np.sqrt(4761)) -points = np.zeros((2, n*n)) -for j in range(0, n): - for k in range(0, n): - points[:, j+k*n] = (j, k) - -# Htool parameters -eta = 10 -epsilon = 1e-3 -minclustersize = 25 - -# Build H matrix -generator = Generator(points) -symmetric = 'S' -UPLO = 'L' -cluster = Htool.PCARegularClustering(2) -cluster.build(n*n, points, 2) -cluster.set_minclustersize(minclustersize) -HMatrix = Htool.HMatrix(cluster, cluster, epsilon, eta, symmetric, UPLO) -HMatrix.build(generator, points) - -# Dense matrix -Full_H = 1.0 / (1e-5 + norm(points.T.reshape(1, n*n, 2) - - points.T.reshape(n*n, 1, 2), axis=2)) - -# GMRES -y = np.ones((n*n,)) -x, _ = gmres(HMatrix, y) -x_full, _ = gmres(Full_H, y) - -err_gmres_hmat = norm(HMatrix @ x - y) -err_gmres_dense = norm(Full_H @ x_full - y) -err_comp = norm(x - x_full)/norm(x) - -comm = MPI.COMM_WORLD -rank = comm.Get_rank() -if rank == 0: - print("Error from gmres (Hmatrix):", err_gmres_hmat) - print("Error from gmres (full matrix):", err_gmres_dense) - print("Error between the two solutions:", err_comp) - -# Several ways to display information -print(HMatrix) -HMatrix.print_infos() -HMatrix.display() -HMatrix.display_cluster(points, 2, "target") diff --git a/example/use_block_jacobi_solver.py b/example/use_block_jacobi_solver.py new file mode 100644 index 0000000..a338d5b --- /dev/null +++ b/example/use_block_jacobi_solver.py @@ -0,0 +1,110 @@ +import copy +import logging + +import Htool +import matplotlib.pyplot as plt +import mpi4py +import numpy as np +from create_geometry import create_random_geometries +from define_custom_generators import CustomGenerator + +logging.basicConfig(level=logging.INFO) + +# Random geometry +size = 500 +dimension = 3 + +[points, _] = create_random_geometries(dimension, size, size) + + +# Htool parameters +eta = 10 +epsilon = 1e-3 +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + points, number_of_children, mpi4py.MPI.COMM_WORLD.size +) + +# Build generator +generator = CustomGenerator(points, points) + +# Build distributed operator +default_approximation = Htool.DefaultApproximationBuilder( + generator, + cluster, + cluster, + epsilon, + eta, + "S", + "L", + mpi4py.MPI.COMM_WORLD, +) +hmatrix = default_approximation.hmatrix +Htool.recompression(hmatrix) + +# Solver with block Jacobi preconditionner +block_diagonal_hmatrix = copy.deepcopy( + default_approximation.block_diagonal_hmatrix +) # inplace factorization requires deepcopy + +default_solver_builder = Htool.DDMSolverBuilder( + default_approximation.distributed_operator, block_diagonal_hmatrix +) +solver = default_solver_builder.solver + + +# Solver with block Jacobi +x_ref = np.random.random(size) +b = default_approximation.distributed_operator * x_ref +x = np.zeros(size) + +hpddm_args = "-hpddm_compute_residual l2 " +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + hpddm_args += "-hpddm_verbosity 10" +solver.set_hpddm_args(hpddm_args) +solver.facto_one_level() +solver.solve(x, b) + + +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +solver_information = solver.get_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref)) + hmatrix = default_approximation.hmatrix + local_block_hmatrix = default_approximation.block_diagonal_hmatrix + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + print(solver_information) + + fig = plt.figure() + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("cluster at depth 1") + ax2.set_title("cluster at depth 2") + ax3.set_title("Hmatrix on rank 0") + ax4.set_title("Block diagonal Hmatrix on rank 0") + Htool.plot(ax1, cluster, points, 1) + Htool.plot(ax2, cluster, points, 2) + Htool.plot(ax3, hmatrix) + Htool.plot(ax4, local_block_hmatrix) + plt.show() diff --git a/example/use_cluster.py b/example/use_cluster.py new file mode 100644 index 0000000..1f3589a --- /dev/null +++ b/example/use_cluster.py @@ -0,0 +1,56 @@ +import Htool +import matplotlib.pyplot as plt +import mpi4py +from create_geometry import create_random_geometries + +# Random geometry +nb_rows = 500 +nb_cols = 500 +dimension = 3 +[target_points, _] = create_random_geometries(dimension, nb_rows, nb_cols) +# [target_points, _, _] = create_partitionned_geometries_test(3, nb_rows, nb_cols, 1) + + +# Parameters +minclustersize = 10 +number_of_children = 2 + +# Cluster builder +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) + +# Strategies +# direction_computation_strategy = Htool.ComputeBoundingBox() +direction_computation_strategy = Htool.ComputeLargestExtent() # default +splitting_strategy = Htool.RegularSplitting() # default +# splitting_strategy = Htool.GeometricSplitting() +cluster_builder.set_direction_computation_strategy(direction_computation_strategy) +cluster_builder.set_splitting_strategy(splitting_strategy) + + +# Build cluster +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + target_points, number_of_children, mpi4py.MPI.COMM_WORLD.size +) + + +# Local cluster +local_target_cluster: Htool.Cluster = target_cluster.get_cluster_on_partition( + mpi4py.MPI.COMM_WORLD.rank +) + +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(1, 2, 1) + ax2 = fig.add_subplot(1, 2, 2) + elif dimension == 3: + ax1 = fig.add_subplot(1, 2, 1, projection="3d") + ax2 = fig.add_subplot(1, 2, 2, projection="3d") + + ax1.set_title("target cluster\ndepth 2") + ax2.set_title("local cluster\ntarget partition number 0\ndepth 2") + Htool.plot(ax1, target_cluster, target_points, 2) + Htool.plot(ax2, local_target_cluster, target_points, 2) + plt.show() diff --git a/example/use_cluster_with_given_partition.py b/example/use_cluster_with_given_partition.py new file mode 100644 index 0000000..8499605 --- /dev/null +++ b/example/use_cluster_with_given_partition.py @@ -0,0 +1,46 @@ +import Htool +import matplotlib.pyplot as plt +import mpi4py +from create_geometry import create_partitionned_geometries + +# Random geometry +nb_rows = 500 +nb_cols = 500 +dimension = 3 +[target_points, _, target_partition] = create_partitionned_geometries( + dimension, nb_rows, nb_cols, mpi4py.MPI.COMM_WORLD.size +) + + +# Parameters +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) + +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + target_points, number_of_children, mpi4py.MPI.COMM_WORLD.size, target_partition +) + +# Local cluster +local_target_cluster: Htool.Cluster = target_cluster.get_cluster_on_partition( + mpi4py.MPI.COMM_WORLD.rank +) + +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(1, 2, 1) + ax2 = fig.add_subplot(1, 2, 2) + elif dimension == 3: + ax1 = fig.add_subplot(1, 2, 1, projection="3d") + ax2 = fig.add_subplot(1, 2, 2, projection="3d") + + ax1.set_title("target cluster\ndepth 2") + ax2.set_title("local cluster\ntarget partition number 0\ndepth 1") + Htool.plot(ax1, target_cluster, target_points, 1) + Htool.plot(ax2, local_target_cluster, target_points, 1) + plt.show() diff --git a/example/use_custom_dense_block_generator.py b/example/use_custom_dense_block_generator.py new file mode 100644 index 0000000..64a6b78 --- /dev/null +++ b/example/use_custom_dense_block_generator.py @@ -0,0 +1,100 @@ +import logging + +import Htool +import mpi4py +import numpy as np +from create_geometry import create_random_geometries +from define_custom_dense_blocks_generator import CustomDenseBlocksGenerator +from define_custom_generators import CustomGenerator +from matplotlib import pyplot as plt + +logging.basicConfig(level=logging.INFO) + +# Random geometry +size = 500 +dimension = 3 +[points, _] = create_random_geometries(dimension, size, size) + + +# Htool parameters +eta = 10 +epsilon = 1e-3 +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + points, number_of_children, mpi4py.MPI.COMM_WORLD.size +) + + +# Build generator +generator = CustomGenerator(points, points) + +# +custom_dense_blocks_generator = CustomDenseBlocksGenerator( + generator, target_cluster, target_cluster +) + +# Build HMatrix +# low_rank_generator = CustomSVD(generator) +hmatrix_builder = Htool.HMatrixBuilder(epsilon, eta, "N", "N") + +hmatrix_builder.set_dense_blocks_generator(custom_dense_blocks_generator) + +# Build distributed operator +distributed_operator_from_hmatrix = Htool.DistributedOperatorFromHMatrix( + generator, target_cluster, target_cluster, hmatrix_builder, mpi4py.MPI.COMM_WORLD +) + +distributed_operator = distributed_operator_from_hmatrix.distributed_operator +hmatrix = distributed_operator_from_hmatrix.hmatrix +Htool.openmp_recompression(hmatrix) + +# Test matrix vector product +np.random.seed(0) +x = np.random.rand(size) +y_1 = distributed_operator * x +y_2 = generator.mat_vec(x) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(y_1 - y_2) / np.linalg.norm(y_2)) + + +# Test matrix matrix product +X = np.asfortranarray(np.random.rand(size, 2)) +Y_1 = distributed_operator @ X +Y_2 = generator.mat_mat(X) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + # ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + # ax3 = fig.add_subplot(2, 2, 3, projection="3d") + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("target cluster at depth 1") + ax2.set_title("target cluster at depth 2") + # ax3.set_title("source cluster at depth 1") + ax4.set_title("Hmatrix on rank 0") + Htool.plot(ax1, target_cluster, points, 1) + Htool.plot(ax2, target_cluster, points, 2) + Htool.plot(ax4, hmatrix) + plt.show() diff --git a/example/use_custom_local_operator.py b/example/use_custom_local_operator.py new file mode 100644 index 0000000..b97d031 --- /dev/null +++ b/example/use_custom_local_operator.py @@ -0,0 +1,64 @@ +import Htool +import mpi4py +import numpy as np +from create_geometry import create_partitionned_geometries +from define_custom_generators import CustomGenerator +from define_custom_local_operator import CustomLocalOperator + +# Random geometry +nb_rows = 500 +nb_cols = 500 +dimension = 3 +[target_points, source_points, target_partition] = create_partitionned_geometries( + dimension, nb_rows, nb_cols, mpi4py.MPI.COMM_WORLD.size +) + + +# Htool parameters +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + target_points, number_of_children, mpi4py.MPI.COMM_WORLD.size, target_partition +) +source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + source_points, number_of_children, mpi4py.MPI.COMM_WORLD.size +) + + +# Build generator +generator = CustomGenerator(target_points, source_points) + +# Build local operator +local_operator = CustomLocalOperator( + generator, + target_cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank), + source_cluster, + "N", + "N", + False, + False, +) + +# Build distributed operator +custom_local_approximation = Htool.CustomApproximationBuilder( + target_cluster, source_cluster, "N", "N", mpi4py.MPI.COMM_WORLD, local_operator +) +distributed_operator = custom_local_approximation.distributed_operator + +# Test matrix vector product +np.random.seed(0) +x = np.random.rand(nb_cols) +y_1 = distributed_operator * x +y_2 = generator.mat_vec(x) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(y_1 - y_2) / np.linalg.norm(y_2)) + + +# Test matrix matrix product +X = np.asfortranarray(np.random.rand(nb_cols, 2)) +Y_1 = distributed_operator @ X +Y_2 = generator.mat_mat(X) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) diff --git a/example/use_custom_low_rank_approximation.py b/example/use_custom_low_rank_approximation.py new file mode 100644 index 0000000..6565822 --- /dev/null +++ b/example/use_custom_low_rank_approximation.py @@ -0,0 +1,105 @@ +import logging + +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 +from matplotlib import pyplot as plt + +logging.basicConfig(level=logging.INFO) + +# Random geometry +nb_rows = 500 +nb_cols = 500 +dimension = 3 +[target_points, source_points, target_partition] = create_partitionned_geometries( + dimension, nb_rows, nb_cols, mpi4py.MPI.COMM_WORLD.size +) + + +# Htool parameters +eta = 100 +epsilon = 1e-3 +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + target_points, number_of_children, mpi4py.MPI.COMM_WORLD.size, target_partition +) +source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + source_points, number_of_children, mpi4py.MPI.COMM_WORLD.size +) + + +# Build generator +generator = CustomGenerator(target_points, source_points) + +# Low rank generator + +low_rank_generator = CustomSVD(generator) + +# Build HMatrix +hmatrix_builder = Htool.HMatrixBuilder(epsilon, eta, "N", "N") +# or hmatrix_builder.set_low_rank_generator(low_rank_generator) + +# Build distributed operator +distributed_operator_from_hmatrix = Htool.DistributedOperatorFromHMatrix( + generator, target_cluster, source_cluster, hmatrix_builder, mpi4py.MPI.COMM_WORLD +) + +distributed_operator = distributed_operator_from_hmatrix.distributed_operator +hmatrix = distributed_operator_from_hmatrix.hmatrix +Htool.openmp_recompression(hmatrix) + +# Test matrix vector product +np.random.seed(0) +x = np.random.rand(nb_cols) +y_1 = distributed_operator * x +y_2 = generator.mat_vec(x) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(y_1 - y_2) / np.linalg.norm(y_2)) + + +# Test matrix matrix product +X = np.asfortranarray(np.random.rand(nb_cols, 2)) +Y_1 = distributed_operator @ X +Y_2 = generator.mat_mat(X) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) + +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + ax3 = fig.add_subplot(2, 2, 3, projection="3d") + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("target cluster at depth 1") + ax2.set_title("target cluster at depth 2") + ax3.set_title("source cluster at depth 1") + ax4.set_title("Hmatrix on rank 0") + Htool.plot(ax1, target_cluster, target_points, 1) + Htool.plot(ax2, target_cluster, target_points, 2) + Htool.plot(ax3, source_cluster, source_points, 1) + Htool.plot(ax4, hmatrix) + plt.show() diff --git a/example/use_hmatrix_compression.py b/example/use_hmatrix_compression.py new file mode 100644 index 0000000..60b6fc5 --- /dev/null +++ b/example/use_hmatrix_compression.py @@ -0,0 +1,104 @@ +import logging + +import Htool +import matplotlib.pyplot as plt +import mpi4py +import numpy as np +from create_geometry import create_partitionned_geometries +from define_custom_generators import CustomGenerator + +logging.basicConfig(level=logging.INFO) + +# Random geometry +nb_rows = 500 +nb_cols = 500 +dimension = 3 +[target_points, source_points, target_partition] = create_partitionned_geometries( + dimension, nb_rows, nb_cols, mpi4py.MPI.COMM_WORLD.size +) + + +# Htool parameters +eta = 10 +epsilon = 1e-3 +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + target_points, number_of_children, mpi4py.MPI.COMM_WORLD.size, target_partition +) +source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + source_points, number_of_children, mpi4py.MPI.COMM_WORLD.size +) + + +# Build generator +generator = CustomGenerator(target_points, source_points) + +# Build distributed operator +default_approximation = Htool.DefaultApproximationBuilder( + generator, + target_cluster, + source_cluster, + epsilon, + eta, + "N", + "N", + mpi4py.MPI.COMM_WORLD, +) + +distributed_operator = default_approximation.distributed_operator +hmatrix = default_approximation.hmatrix +Htool.openmp_recompression(hmatrix) + + +# Test matrix vector product +np.random.seed(0) +x = np.random.rand(nb_cols) +y_1 = distributed_operator * x +y_2 = generator.mat_vec(x) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(y_1 - y_2) / np.linalg.norm(y_2)) + + +# Test matrix matrix product +X = np.asfortranarray(np.random.rand(nb_cols, 2)) +Y_1 = distributed_operator @ X +Y_2 = generator.mat_mat(X) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) + +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + ax3 = fig.add_subplot(2, 2, 3, projection="3d") + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("target cluster at depth 1") + ax2.set_title("target cluster at depth 2") + ax3.set_title("source cluster at depth 1") + ax4.set_title("Hmatrix on rank 0") + Htool.plot(ax1, target_cluster, target_points, 1) + Htool.plot(ax2, target_cluster, target_points, 2) + Htool.plot(ax3, source_cluster, source_points, 1) + Htool.plot(ax4, hmatrix) + plt.show() diff --git a/example/use_local_hmatrix_compression.py b/example/use_local_hmatrix_compression.py new file mode 100644 index 0000000..10e7387 --- /dev/null +++ b/example/use_local_hmatrix_compression.py @@ -0,0 +1,180 @@ +import logging + +import Htool +import matplotlib.pyplot as plt +import mpi4py +import numpy as np +from create_geometry import create_random_geometries +from define_custom_generators import CustomGenerator +from define_custom_local_operator import CustomLocalOperator + +logging.basicConfig(level=logging.INFO) + +# Random geometry +target_size = 500 +source_size = 500 +dimension = 3 +# [points, _, partition] = create_partitionned_geometries( +# dimension, size, size, mpi4py.MPI.COMM_WORLD.size +# ) + +[target_points, source_points] = create_random_geometries( + dimension, target_size, source_size +) + +# Htool parameters +eta = 10 +epsilon = 1e-3 +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + target_points, + number_of_children, + mpi4py.MPI.COMM_WORLD.size, +) + +source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + source_points, + number_of_children, + mpi4py.MPI.COMM_WORLD.size, +) + +local_target_cluster: Htool.Cluster = target_cluster.get_cluster_on_partition( + mpi4py.MPI.COMM_WORLD.rank +) +local_source_cluster: Htool.Cluster = source_cluster.get_cluster_on_partition( + mpi4py.MPI.COMM_WORLD.rank +) + +# +source_permutation = source_cluster.get_permutation() +permuted_source_points = np.zeros((dimension, source_size)) +for i in range(0, source_size): + permuted_source_points[:, i] = source_points[:, source_permutation[i]] + +# Build generator +generator = CustomGenerator(target_points, source_points) + + +# Build distributed operator +default_local_approximation = Htool.DefaultLocalApproximationBuilder( + generator, + target_cluster, + source_cluster, + epsilon, + eta, + "N", + "N", + mpi4py.MPI.COMM_WORLD, +) +distributed_operator = default_local_approximation.distributed_operator +hmatrix = default_local_approximation.hmatrix +Htool.recompression(hmatrix) + + +# Build off diagonal operators +off_diagonal_nc_1 = local_source_cluster.get_offset() +off_diagonal_nc_2 = ( + source_cluster.get_size() + - local_source_cluster.get_size() + - local_source_cluster.get_offset() +) +local_nc = local_source_cluster.get_size() + +off_diagonal_partition = np.zeros((2, 2), dtype=int) +off_diagonal_partition[0, 0] = 0 +off_diagonal_partition[1, 0] = off_diagonal_nc_1 +off_diagonal_partition[0, 1] = off_diagonal_nc_1 + local_nc +off_diagonal_partition[1, 1] = off_diagonal_nc_2 +off_diagonal_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + permuted_source_points, number_of_children, 2, off_diagonal_partition +) + +off_diagonal_generator = CustomGenerator(target_points, permuted_source_points) + +local_operator_1 = None +if off_diagonal_nc_1 > 0: + local_operator_1 = CustomLocalOperator( + off_diagonal_generator, + local_target_cluster, + off_diagonal_cluster.get_cluster_on_partition(0), + "N", + "N", + False, + True, + ) + +local_operator_2 = None +if off_diagonal_nc_2 > 0: + local_operator_2 = CustomLocalOperator( + off_diagonal_generator, + local_target_cluster, + off_diagonal_cluster.get_cluster_on_partition(1), + "N", + "N", + False, + True, + ) + +if local_operator_1: + distributed_operator.add_local_operator(local_operator_1) +if local_operator_2: + distributed_operator.add_local_operator(local_operator_2) + +# Test matrix vector product +np.random.seed(0) +x = np.random.rand(source_size) +x = np.ones(source_size) +y_1 = distributed_operator * x +y_2 = generator.mat_vec(x) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(y_1 - y_2) / np.linalg.norm(y_2)) + + +# Test matrix matrix product +X = np.asfortranarray(np.random.rand(source_size, 2)) +Y_1 = distributed_operator @ X +Y_2 = generator.mat_mat(X) +print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) + +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + + fig = plt.figure() + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + ax3 = fig.add_subplot(2, 2, 3, projection="3d") + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("source cluster at depth 1") + ax2.set_title("source cluster at depth 2") + ax3.set_title("off diagonal cluster on rank 0 at depth 2") + ax4.set_title("Hmatrix on rank 0") + Htool.plot(ax1, source_cluster, source_points, 1) + Htool.plot(ax2, source_cluster, source_points, 2) + if mpi4py.MPI.COMM_WORLD.Get_size() > 1: + Htool.plot( + ax3, + off_diagonal_cluster.get_cluster_on_partition(1), + permuted_source_points, + 2, + ) + Htool.plot(ax4, hmatrix) + plt.show() diff --git a/lib/hpddm b/lib/hpddm index 9c5092c..5890d5a 160000 --- a/lib/hpddm +++ b/lib/hpddm @@ -1 +1 @@ -Subproject commit 9c5092ca711678ebc257ffadc7d29057d4566a87 +Subproject commit 5890d5addf3962d539dc25c441ec3ff4af93b3ab diff --git a/lib/htool b/lib/htool index b6e9169..c71b1e0 160000 --- a/lib/htool +++ b/lib/htool @@ -1 +1 @@ -Subproject commit b6e91690f8d7c20d67dfb3b8db2e7ea674405a37 +Subproject commit c71b1e0a8982a0699843dad36c1c7bddb92774cb diff --git a/lib/htool_generate_data_test b/lib/htool_generate_data_test new file mode 160000 index 0000000..b3b8909 --- /dev/null +++ b/lib/htool_generate_data_test @@ -0,0 +1 @@ +Subproject commit b3b8909a5a831609216298a0c3219a3be275a8ae diff --git a/lib/pybind11 b/lib/pybind11 index af8849f..0e43fcc 160000 --- a/lib/pybind11 +++ b/lib/pybind11 @@ -1 +1 @@ -Subproject commit af8849f47eea35fe16df0c5df64221b33b655968 +Subproject commit 0e43fcc75e6b7429e3511dfb44343ec05a0ab843 diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..5537558 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,54 @@ +[project] +name = "htool" +version = "0.1" +authors = [{ name = "Pierre Marchand", email = "test@test.com" }] +description = "CLI utility to automate asciinema" +readme = "README.md" +requires-python = ">=3.7" +license = { text = "MIT License" } +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: MacOS", + "Operating System :: Unix", +] +keywords = [ + "kernel approximation", + "hierarchical matrix", + "domain decomposition method", +] +dependencies = ["mpi4py", "numpy"] + + +[project.optional-dependencies] +dev = ["ruff", "matplotlib>=3.0.0", "pytest", "scipy"] + + +[project.urls] +"Homepage" = "https://github.com/htool-ddm/htool_python" +"Bug Tracker" = "https://github.com/htool-ddm/htool_python/issues" + +[tool.ruff] +line-length = 88 +indent-width = 4 +exclude = ["lib", "cmake-format.py"] + +[tool.ruff.lint] +select = [ + "E", # pycodestyle + "F", # Pyflakes + "UP", # pyupgrade + "B", # flake8-bugbear + "SIM", # flake8-simplify + "I", # isort +] +ignore = [] + +# Allow fix for all enabled rules (when `--fix`) is provided. +fixable = ["ALL"] +unfixable = [] + + +[build-system] +requires = ["setuptools>=42", "wheel", "cmake", "mpi4py", "numpy"] +build-backend = "setuptools.build_meta" diff --git a/setup.py b/setup.py index b24df7e..5732f45 100644 --- a/setup.py +++ b/setup.py @@ -1,85 +1,173 @@ import os import re -import sys -import platform import subprocess +import sys +from pathlib import Path -from setuptools import setup, Extension +from setuptools import Extension, setup from setuptools.command.build_ext import build_ext -from distutils.version import LooseVersion + +# Convert distutils Windows platform specifiers to CMake -A arguments +PLAT_TO_CMAKE = { + "win32": "Win32", + "win-amd64": "x64", + "win-arm32": "ARM", + "win-arm64": "ARM64", +} class CMakeExtension(Extension): - def __init__(self, name, sourcedir=''): - Extension.__init__(self, name, sources=[]) - self.sourcedir = os.path.abspath(sourcedir) + def __init__(self, name: str, sourcedir: str = "") -> None: + super().__init__(name, sources=[]) + self.sourcedir = os.fspath(Path(sourcedir).resolve()) class CMakeBuild(build_ext): - def run(self): - try: - out = subprocess.check_output(['cmake', '--version']) - except OSError: - raise RuntimeError("CMake must be installed to build the following extensions: " + - ", ".join(e.name for e in self.extensions)) - - if platform.system() == "Windows": - cmake_version = LooseVersion( - re.search(r'version\s*([\d.]+)', out.decode()).group(1)) - if cmake_version < '3.1.0': - raise RuntimeError("CMake >= 3.1.0 is required on Windows") - - for ext in self.extensions: - self.build_extension(ext) - def build_extension(self, ext): - extdir = os.path.abspath(os.path.dirname( - self.get_ext_fullpath(ext.name))) - # required for auto-detection of auxiliary "native" libs - if not extdir.endswith(os.path.sep): - extdir += os.path.sep - - cmake_args = ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, - '-DPYTHON_EXECUTABLE=' + sys.executable] - - cfg = 'Debug' if self.debug else 'Release' - build_args = ['--config', cfg] - - if platform.system() == "Windows": - cmake_args += [ - '-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{}={}'.format(cfg.upper(), extdir)] - if sys.maxsize > 2**32: - cmake_args += ['-A', 'x64'] - build_args += ['--', '/m'] + # Must be in this form due to bug in .resolve() only fixed in Python 3.10+ + ext_fullpath = Path.cwd() / self.get_ext_fullpath(ext.name) + extdir = ext_fullpath.parent.resolve() + + # Using this requires trailing slash for auto-detection & inclusion of + # auxiliary "native" libs + + debug = int(os.environ.get("DEBUG", 0)) if self.debug is None else self.debug + cfg = "Debug" if debug else "Release" + + # CMake lets you override the generator - we need to check this. + # Can be set with Conda-Build, for example. + cmake_generator = os.environ.get("CMAKE_GENERATOR", "") + + # Set Python_EXECUTABLE instead if you use PYBIND11_FINDPYTHON + # EXAMPLE_VERSION_INFO shows you how to pass a value into the C++ code + # from Python. + cmake_args = [ + f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY={extdir}{os.sep}", + f"-DPYTHON_EXECUTABLE={sys.executable}", + f"-DCMAKE_BUILD_TYPE={cfg}", # not used on MSVC, but no harm + ] + build_args = [] + # Adding CMake arguments set as environment variable + # (needed e.g. to build for ARM OSx on conda-forge) + if "CMAKE_ARGS" in os.environ: + cmake_args += [item for item in os.environ["CMAKE_ARGS"].split(" ") if item] + + # # In this example, we pass in the version to C++. You might not need to. + # cmake_args += [f"-DEXAMPLE_VERSION_INFO={self.distribution.get_version()}"] + + if self.compiler.compiler_type != "msvc": + # Using Ninja-build since it a) is available as a wheel and b) + # multithreads automatically. MSVC would require all variables be + # exported for Ninja to pick it up, which is a little tricky to do. + # Users can override the generator with CMAKE_GENERATOR in CMake + # 3.15+. + if cmake_generator == "Ninja": + try: + import ninja + + ninja_executable_path = Path(ninja.BIN_DIR) / "ninja" + cmake_args += [ + "-GNinja", + f"-DCMAKE_MAKE_PROGRAM:FILEPATH={ninja_executable_path}", + ] + except ImportError: + pass + else: - cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] - build_args += ['--', '-j2'] + # Single config generators are handled "normally" + single_config = any(x in cmake_generator for x in {"NMake", "Ninja"}) + + # CMake allows an arch-in-generator style for backward compatibility + contains_arch = any(x in cmake_generator for x in {"ARM", "Win64"}) + + # Specify the arch if using MSVC generator, but only if it doesn't + # contain a backward-compatibility arch spec already in the + # generator name. + if not single_config and not contains_arch: + cmake_args += ["-A", PLAT_TO_CMAKE[self.plat_name]] + + # Multi-config generators have a different way to specify configs + if not single_config: + cmake_args += [ + f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{cfg.upper()}={extdir}" + ] + build_args += ["--config", cfg] + + if sys.platform.startswith("darwin"): + # Cross-compile support for macOS - respect ARCHFLAGS if set + archs = re.findall(r"-arch (\S+)", os.environ.get("ARCHFLAGS", "")) + if archs: + cmake_args += ["-DCMAKE_OSX_ARCHITECTURES={}".format(";".join(archs))] + + # Set CMAKE_BUILD_PARALLEL_LEVEL to control the parallel build level + # across all generators. + if ( + "CMAKE_BUILD_PARALLEL_LEVEL" not in os.environ + and hasattr(self, "parallel") + and self.parallel + ): + # self.parallel is a Python 3 only way to set parallel jobs by hand + # using -j in the build_ext call, not supported by pip or PyPA-build. + + # CMake 3.12+ only. + build_args += [f"-j{self.parallel}"] + + build_temp = Path(self.build_temp) / ext.name + if not build_temp.exists(): + build_temp.mkdir(parents=True) + + subprocess.run( + ["cmake", ext.sourcedir, *cmake_args], cwd=build_temp, check=True + ) + subprocess.run( + ["cmake", "--build", ".", *build_args], cwd=build_temp, check=True + ) + + # debug = int(os.environ.get("DEBUG", 0)) if self.debug is None else self.debug + # cfg = "Debug" if debug else "Release" + + # cmake_args = [ + # "-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=" + extdir, + # "-DPYTHON_EXECUTABLE=" + sys.executable, + # ] + # cfg = "Debug" if self.debug else "Release" + # build_args = ["--config", cfg] - env = os.environ.copy() - env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format( - env.get('CXXFLAGS', ''), self.distribution.get_version()) - if not os.path.exists(self.build_temp): - os.makedirs(self.build_temp) - subprocess.check_call(['cmake', ext.sourcedir] + - cmake_args, cwd=self.build_temp, env=env) - subprocess.check_call(['cmake', '--build', '.'] + - build_args, cwd=self.build_temp) + # if platform.system() == "Windows": + # cmake_args += [ + # "-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{}={}".format(cfg.upper(), extdir) + # ] + # if sys.maxsize > 2**32: + # cmake_args += ["-A", "x64"] + # build_args += ["--", "/m"] + # else: + # cmake_args += ["-DCMAKE_BUILD_TYPE=" + cfg] + # if self.debug: + # cmake_args += ["-DUSE_SANITIZER=Address"] + # build_args += ["--", "-j2"] + # print("BOUH", cmake_args) + # env = os.environ.copy() + # env["CXXFLAGS"] = '{} -DVERSION_INFO=\\"{}\\"'.format( + # env.get("CXXFLAGS", ""), self.distribution.get_version() + # ) + # if not os.path.exists(self.build_temp): + # os.makedirs(self.build_temp) + # subprocess.check_call( + # ["cmake", ext.sourcedir] + cmake_args, cwd=self.build_temp, env=env + # ) + # subprocess.check_call( + # ["cmake", "--build", "."] + build_args, cwd=self.build_temp + # ) setup( - name='Htool', - version='1.0.0', - author='Pierre Marchand', - author_email='', - description='A pybind11 interface to Htool, a header only c++ library that provides Hierarchical matrices.', - long_description='', - ext_modules=[CMakeExtension('Htool')], + name="Htool", + version="0.9.0", + author="Pierre Marchand", + author_email="", + description="""A pybind11 interface to Htool.""", + long_description="", + ext_modules=[CMakeExtension("Htool")], cmdclass=dict(build_ext=CMakeBuild), zip_safe=False, - install_requires=[ - 'numpy>=1.0.0', - 'scipy>=1.0.0', - 'mpi4py>=3.0.0', - 'matplotlib>=3.0.0' - ], ) diff --git a/src/htool/cluster.hpp b/src/htool/cluster.hpp deleted file mode 100644 index 6fbb7e5..0000000 --- a/src/htool/cluster.hpp +++ /dev/null @@ -1,146 +0,0 @@ -#ifndef HTOOL_CLUSTER_CPP -#define HTOOL_CLUSTER_CPP - -#include -#include -#include -#include - -#include "misc.hpp" -#include "wrapper_mpi.hpp" -#include - -namespace py = pybind11; -using namespace pybind11::literals; -using namespace htool; - -template -void declare_Cluster(py::module &m, const std::string &className) { - - py::class_, VirtualCluster> py_class(m, className.c_str()); - py_class.def(py::init()); - py_class.def( - "build", [](ClusterType &self, int nb_pt, py::array_t x, int nb_sons, MPI_Comm_wrapper comm) { - if (x.ndim() != 2 && x.shape()[0] != self.get_space_dim()) { - throw std::runtime_error("Wrong dimension for x"); // LCOV_EXCL_LINE - } - self.build(nb_pt, x.data(), nb_sons, comm); - }, - "nb_pt"_a, - "x"_a, - "nb_sons"_a = 2, - py::arg("comm") = MPI_Comm_wrapper(MPI_COMM_WORLD)); - py_class.def( - "build", [](ClusterType &self, int nb_pt, py::array_t x, py::array_t MasterOffset, int nb_sons, MPI_Comm_wrapper comm) { - if (x.ndim() != 2 && x.shape()[0] != self.get_space_dim()) { - throw std::runtime_error("Wrong dimension for x"); // LCOV_EXCL_LINE - } - if (MasterOffset.ndim() != 2 && MasterOffset.shape()[0] != 2) { - throw std::runtime_error("Wrong dimension for MasterOffset"); // LCOV_EXCL_LINE - } - self.build(nb_pt, x.data(), MasterOffset.data(), nb_sons, comm); - }, - "nb_pt"_a, - "x"_a, - "MasterOffset"_a, - "nb_sons"_a = 2, - py::arg("comm") = MPI_Comm_wrapper(MPI_COMM_WORLD)); - py_class.def("get_masteroffset", [](ClusterType &self) { - std::vector> masteroffset = self.get_masteroffset(); - std::vector masteroffset_raw(masteroffset.size() * 2); - std::array shape{2, masteroffset.size()}; - for (int i = 0; i < masteroffset.size(); i++) { - masteroffset_raw[2 * i] = masteroffset[i].first; - masteroffset_raw[2 * i + 1] = masteroffset[i].second; - } - return py::array_t(shape, masteroffset_raw.data()); - }); - py_class.def( - "display", [](ClusterType &self, py::array_t x, int depth, bool show, MPI_Comm_wrapper comm) { - int rankWorld; - MPI_Comm_rank(comm, &rankWorld); - - if (rankWorld == 0) { - - VirtualCluster const *root = self.get_root(); - - std::stack s; - s.push(root); - - int size = root->get_size(); - int space_dim = root->get_space_dim(); - std::vector output((space_dim + 1) * size); - - // Permuted geometric points - for (int i = 0; i < size; ++i) { - for (int p = 0; p < space_dim; p++) { - output[i + size * p] = x.at(p, root->get_perm(i)); - } - } - - int counter = 0; - while (!s.empty()) { - VirtualCluster const *curr = s.top(); - s.pop(); - - if (depth == curr->get_depth()) { - std::fill_n(&(output[space_dim * size + curr->get_offset()]), curr->get_size(), counter); - counter += 1; - } - - // Recursion - if (!curr->IsLeaf()) { - - for (int p = 0; p < curr->get_nb_sons(); p++) { - s.push(&(curr->get_son(p))); - } - } - } - - // Import - py::object plt = py::module::import("matplotlib.pyplot"); - py::object colors = py::module::import("matplotlib.colors"); - - // Create Color Map - py::object colormap = plt.attr("get_cmap")("Dark2"); - py::object norm = colors.attr("Normalize")("vmin"_a = (*std::min_element(output.begin() + space_dim * size, output.end())), "vmax"_a = (*std::max_element(output.begin() + space_dim * size, output.end()))); - - // Figure - py::object fig = plt.attr("figure")(); - - if (space_dim == 2) { - py::object ax = fig.attr("add_subplot")(111); - ax.attr("scatter")(std::vector(output.begin(), output.begin() + size), std::vector(output.begin() + size, output.begin() + 2 * size), "c"_a = colormap(norm(std::vector(output.begin() + 2 * size, output.end()))), "marker"_a = 'o'); - - } else if (space_dim == 3) { - py::object Axes3D = py::module::import("mpl_toolkits.mplot3d").attr("Axes3D"); - - py::object ax = fig.attr("add_subplot")(111, "projection"_a = "3d"); - ax.attr("scatter")(std::vector(output.begin(), output.begin() + size), std::vector(output.begin() + size, output.begin() + 2 * size), std::vector(output.begin() + 2 * size, output.begin() + 3 * size), "c"_a = colormap(norm(std::vector(output.begin() + 3 * size, output.end()))), "marker"_a = 'o'); - } - - plt.attr("draw")(); - if (show) { - plt.attr("show")(); // LCOV_EXCL_LINE - } - return 0; - } - - return 0; - }, - "x"_a, - "depth"_a, - py::arg("show") = true, - py::arg("comm") = MPI_Comm_wrapper(MPI_COMM_WORLD)); - py_class.def( - "read_cluster", [](ClusterType &self, std::string file_permutation, std::string file_tree, MPI_Comm_wrapper comm) { - self.read_cluster(file_permutation, file_tree, comm); - }, - py::arg("file_permutation"), - py::arg("file_tree"), - py::arg("comm") = MPI_Comm_wrapper(MPI_COMM_WORLD)); - - py_class.def("set_minclustersize", &ClusterType::set_minclustersize); -} - -#endif \ No newline at end of file diff --git a/src/htool/clustering/cluster_builder.hpp b/src/htool/clustering/cluster_builder.hpp new file mode 100644 index 0000000..76557b4 --- /dev/null +++ b/src/htool/clustering/cluster_builder.hpp @@ -0,0 +1,36 @@ +#ifndef HTOOL_CLUSTERING_CLUSTER_BUILDER_CPP +#define HTOOL_CLUSTERING_CLUSTER_BUILDER_CPP + +#include "cluster_node.hpp" +#include +#include +#include +#include +#include + +namespace py = pybind11; +using namespace pybind11::literals; +using namespace htool; + +template +void declare_cluster_builder(py::module &m, const std::string &className) { + + using Class = ClusterTreeBuilder; + py::class_ py_class(m, className.c_str()); + + py_class.def(py::init<>()); + py_class.def("create_cluster_tree", [](Class &self, py::array_t coordinates, int number_of_children, int size_of_partition) { + return self.create_cluster_tree(coordinates.shape()[1], coordinates.shape()[0], coordinates.data(), number_of_children, size_of_partition); + }); + py_class.def("create_cluster_tree", [](Class &self, py::array_t coordinates, int number_of_children, int size_of_partition, py::array_t partition) { + if (partition.ndim() != 2 && partition.shape()[0] != 2) { + throw std::runtime_error("Wrong format for partition"); // LCOV_EXCL_LINE + } + return self.create_cluster_tree(coordinates.shape()[1], coordinates.shape()[0], coordinates.data(), number_of_children, size_of_partition, partition.data()); + }); + py_class.def("set_minclustersize", &Class::set_minclustersize); + py_class.def("set_direction_computation_strategy", &Class::set_direction_computation_strategy); + py_class.def("set_splitting_strategy", &Class::set_splitting_strategy); +} + +#endif diff --git a/src/htool/clustering/cluster_node.hpp b/src/htool/clustering/cluster_node.hpp new file mode 100644 index 0000000..c1600a4 --- /dev/null +++ b/src/htool/clustering/cluster_node.hpp @@ -0,0 +1,30 @@ +#ifndef HTOOL_CLUSTER_CPP +#define HTOOL_CLUSTER_CPP +#define PYBIND11_DETAILED_ERROR_MESSAGES +#include +#include +#include +#include +#include +namespace py = pybind11; +using namespace pybind11::literals; +using namespace htool; + +template +void declare_cluster_node(py::module &m, const std::string &className) { + + using Class = Cluster; + py::class_ py_class(m, className.c_str()); + + py_class.def("get_size", &Class::get_size); + py_class.def("get_offset", &Class::get_offset); + py_class.def("get_minclustersize", &Class::get_minclustersize); + py_class.def("get_permutation", [](const Class &self) { + auto &permutation = self.get_permutation(); + return py::array_t(std::array{permutation.size()}, permutation.data(), py::capsule(permutation.data())); + ; + }); + py_class.def("get_cluster_on_partition", &Class::get_cluster_on_partition, py::return_value_policy::reference_internal); +} + +#endif diff --git a/src/htool/clustering/implementation/direction_computation.hpp b/src/htool/clustering/implementation/direction_computation.hpp new file mode 100644 index 0000000..005df13 --- /dev/null +++ b/src/htool/clustering/implementation/direction_computation.hpp @@ -0,0 +1,19 @@ +#ifndef HTOOL_PYBIND11_CLUSTERING_IMPLEMENTATION_DIRECTION_COMPUTATION_HPP +#define HTOOL_PYBIND11_CLUSTERING_IMPLEMENTATION_DIRECTION_COMPUTATION_HPP + +#include + +template +void declare_compute_largest_extent(py::module &m) { + using Class = htool::ComputeLargestExtent; + py::class_, htool::VirtualDirectionComputationStrategy> py_class(m, "ComputeLargestExtent"); + py_class.def(py::init<>()); +} + +template +void declare_compute_bounding_box(py::module &m) { + using Class = htool::ComputeBoundingBox; + py::class_, htool::VirtualDirectionComputationStrategy> py_class(m, "ComputeBoundingBox"); + py_class.def(py::init<>()); +} +#endif diff --git a/src/htool/clustering/implementation/splitting.hpp b/src/htool/clustering/implementation/splitting.hpp new file mode 100644 index 0000000..107afe0 --- /dev/null +++ b/src/htool/clustering/implementation/splitting.hpp @@ -0,0 +1,19 @@ +#ifndef HTOOL_PYBIND11_CLUSTERING_IMPLEMENTATION_SPLITTING_HPP +#define HTOOL_PYBIND11_CLUSTERING_IMPLEMENTATION_SPLITTING_HPP + +#include + +template +void declare_regular_splitting(py::module &m) { + using Class = htool::RegularSplitting; + py::class_, htool::VirtualSplittingStrategy> py_class(m, "RegularSplitting"); + py_class.def(py::init<>()); +} + +template +void declare_geometric_splitting(py::module &m) { + using Class = htool::GeometricSplitting; + py::class_, htool::VirtualSplittingStrategy> py_class(m, "GeometricSplitting"); + py_class.def(py::init<>()); +} +#endif diff --git a/src/htool/clustering/interface/direction_computation.hpp b/src/htool/clustering/interface/direction_computation.hpp new file mode 100644 index 0000000..df60a04 --- /dev/null +++ b/src/htool/clustering/interface/direction_computation.hpp @@ -0,0 +1,11 @@ +#ifndef HTOOL_PYBIND11_CLUSTERING_INTERFACE_DIRECTION_COMPUTATION_HPP +#define HTOOL_PYBIND11_CLUSTERING_INTERFACE_DIRECTION_COMPUTATION_HPP + +#include + +template +void declare_interface_direction_computation(py::module &m) { + using Class = htool::VirtualDirectionComputationStrategy; + py::class_> py_class(m, "IDirectionComputation"); +} +#endif diff --git a/src/htool/clustering/interface/splitting.hpp b/src/htool/clustering/interface/splitting.hpp new file mode 100644 index 0000000..ce97737 --- /dev/null +++ b/src/htool/clustering/interface/splitting.hpp @@ -0,0 +1,11 @@ +#ifndef HTOOL_PYBIND11_CLUSTERING_INTERFACE_SPLITTING_HPP +#define HTOOL_PYBIND11_CLUSTERING_INTERFACE_SPLITTING_HPP + +#include + +template +void declare_interface_splitting(py::module &m) { + using Class = htool::VirtualSplittingStrategy; + py::class_> py_class(m, "ISplitting"); +} +#endif diff --git a/src/htool/clustering/utility.hpp b/src/htool/clustering/utility.hpp new file mode 100644 index 0000000..ee7c00e --- /dev/null +++ b/src/htool/clustering/utility.hpp @@ -0,0 +1,13 @@ +#ifndef HTOOL_PYBIND11_CLUSTER_UTILITY_CPP +#define HTOOL_PYBIND11_CLUSTER_UTILITY_CPP + +#include +#include +#include + +template +void declare_cluster_utility(py::module &m) { + m.def("read_cluster_from", &htool::read_cluster_tree); +} + +#endif diff --git a/src/htool/ddm_solver.hpp b/src/htool/ddm_solver.hpp deleted file mode 100644 index 770fa67..0000000 --- a/src/htool/ddm_solver.hpp +++ /dev/null @@ -1,101 +0,0 @@ -#ifndef HTOOL_DDM_SOLVER_CPP -#define HTOOL_DDM_SOLVER_CPP - -#include -#include -#include -#include - -#include "htool/solvers/ddm.hpp" -#include "matrix.hpp" -#include "wrapper_mpi.hpp" -#include - -namespace py = pybind11; -using namespace pybind11::literals; -using namespace htool; - -template -void declare_DDM(py::module &m, const std::string &className) { - - using Class = DDM; - py::class_ py_class(m, className.c_str()); - py_class.def(py::init *const>()); - py_class.def(py::init< - const VirtualGeneratorCpp &, - const VirtualHMatrix *const, - const std::vector &, - const std::vector &, - 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\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); - py_class.def( - "solve", [](Class &self, py::array_t x, const py::array_t b, std::string hpddm_args) { - int rank; - MPI_Comm_rank(self.get_comm(), &rank); - - // HPDDM arguments - HPDDM::Option &opt = *HPDDM::Option::get(); - opt.parse(hpddm_args); - int mu; - - if (b.ndim() == 1 && x.ndim() == 1) { - mu = 1; - } else if ((b.ndim() == 2 && x.ndim() == 2) && b.shape()[1] == x.shape()[1]) { - mu = b.shape()[1]; - } else { - std::string rhs = "("; // LCOV_EXCL_START - std::string sol = "("; - for (int p = 0; p < b.ndim(); p++) { - rhs += htool::NbrToStr(b.shape()[p]); - if (p != b.ndim() - 1) { - rhs += ","; - } - } - rhs += ")"; - for (int p = 0; p < x.ndim(); p++) { - sol += htool::NbrToStr(x.shape()[p]); - if (p != x.ndim() - 1) { - sol += ","; - } - } - sol += ")"; - throw std::invalid_argument("Wrong dimension for right-hand side or solution\nright-hand side: " + rhs + "\n" + "solution: " + sol + "\n"); - // LCOV_EXCL_STOP - } - if (b.shape()[0] != self.get_nb_cols()) { - throw std::invalid_argument("Wrong size for right-hand side"); // LCOV_EXCL_LINE - } - if (x.shape()[0] != self.get_nb_rows()) { - throw std::invalid_argument("Wrong size for solution"); // LCOV_EXCL_LINE - } - - self.solve(b.data(), x.mutable_data(), mu); - }, - py::arg("x").noconvert(true), - py::arg("b"), - py::arg("hpddm_args") = ""); - py_class.def("set_hpddm_args", [](Class &self, std::string hpddm_args) { - int rank; - MPI_Comm_rank(self.get_comm(), &rank); - HPDDM::Option &opt = *HPDDM::Option::get(); - opt.parse(hpddm_args); - }); - py_class.def("print_infos", &Class::print_infos); - py_class.def("get_infos", &Class::get_infos); -} - -#endif \ No newline at end of file diff --git a/src/htool/dense_blocks_generator.hpp b/src/htool/dense_blocks_generator.hpp deleted file mode 100644 index 75e4ad6..0000000 --- a/src/htool/dense_blocks_generator.hpp +++ /dev/null @@ -1,67 +0,0 @@ -#ifndef HTOOL_DENSE_BLOCKS_GENERATOR_CPP -#define HTOOL_DENSE_BLOCKS_GENERATOR_CPP - -#include -#include -#include -#include -#include -#include - -namespace py = pybind11; -using namespace pybind11::literals; -using namespace htool; - -template -class VirtualDenseBlocksGeneratorCpp : public VirtualDenseBlocksGenerator { - - public: - using VirtualDenseBlocksGenerator::VirtualDenseBlocksGenerator; - - void copy_dense_blocks(const std::vector &M, const std::vector &N, const std::vector &rows, const std::vector &cols, std::vector &ptr) const override { - - int nb_blocks = M.size(); - std::vector> vec_ptr; - std::vector> rows_ptr; - std::vector> cols_ptr; - for (int i = 0; i < nb_blocks; i++) { - rows_ptr.emplace_back(std::array{M[i]}, rows[i], py::capsule(rows[i])); - cols_ptr.emplace_back(std::array{N[i]}, cols[i], py::capsule(cols[i])); - vec_ptr.emplace_back(std::array{M[i], N[i]}, ptr[i], py::capsule(ptr[i])); - } - - build_dense_blocks(rows_ptr, cols_ptr, vec_ptr); - } - - // lcov does not see it because of trampoline I assume - virtual void build_dense_blocks(const std::vector> &rows, const std::vector> &cols, std::vector> &blocks) const = 0; // LCOV_EXCL_LINE -}; - -template -class PyVirtualDenseBlocksGenerator : public VirtualDenseBlocksGeneratorCpp { - public: - using VirtualDenseBlocksGeneratorCpp::VirtualDenseBlocksGeneratorCpp; - // PyVirtualGenerator(int nr0, int nc0): IMatrix(nr0,nc0){} - - /* Trampoline (need one for each virtual function) */ - virtual void build_dense_blocks(const std::vector> &rows, const std::vector> &cols, std::vector> &blocks) const override { - PYBIND11_OVERRIDE_PURE( - void, /* Return type */ - VirtualDenseBlocksGeneratorCpp, /* Parent class */ - build_dense_blocks, /* Name of function in C++ (must match Python name) */ - rows, - cols, - blocks /* Argument(s) */ - ); - } -}; - -template -void declare_custom_VirtualDenseBlocksGenerator(py::module &m, const std::string &className) { - using Class = VirtualDenseBlocksGeneratorCpp; - py::class_, PyVirtualDenseBlocksGenerator> py_class(m, className.c_str()); - py_class.def(py::init<>()); - py_class.def("build_dense_blocks", &Class::build_dense_blocks); -} - -#endif diff --git a/src/htool/distributed_operator/distributed_operator.hpp b/src/htool/distributed_operator/distributed_operator.hpp new file mode 100644 index 0000000..08b6215 --- /dev/null +++ b/src/htool/distributed_operator/distributed_operator.hpp @@ -0,0 +1,55 @@ +#ifndef HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_HPP +#define HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_HPP +#include "../misc/utility.hpp" +#include "../misc/wrapper_mpi.hpp" +#include +#include +#include + +template +void declare_distributed_operator(py::module &m, const std::string &class_name) { + using Class = DistributedOperator; + + py::class_ py_class(m, class_name.c_str()); + py_class.def(py::init &, IPartition &, char, char, MPI_Comm_wrapper>(), py::keep_alive<1, 2>(), py::keep_alive<1, 3>()); + py_class.def("add_local_operator", &Class::add_local_operator, py::keep_alive<1, 2>()); + + // Linear algebra + py_class.def( + "__mul__", [](const Class &self, const py::array_t input) { + if (input.ndim() != 1) { + throw std::runtime_error("Wrong dimension for HMatrix-vector product"); // LCOV_EXCL_LINE + } + if (input.shape()[0] != self.get_source_partition().get_global_size()) { + throw std::runtime_error("Wrong size for HMatrix-vector product"); // LCOV_EXCL_LINE + } + py::array_t result(std::array{self.get_target_partition().get_global_size()}); + self.vector_product_global_to_global(input.data(), result.mutable_data()); + + return result; + }, + "in"_a); + + py_class.def( + "__matmul__", [](const Class &self, py::array_t input) { + int mu; + if (input.ndim() == 2) { + mu = input.shape()[1]; + } else { + throw std::runtime_error("Wrong dimension for HMatrix-matrix product"); // LCOV_EXCL_LINE + } + if (input.shape()[0] != self.get_source_partition().get_global_size()) { + throw std::runtime_error("Wrong size for HMatrix-matrix product"); // LCOV_EXCL_LINE + } + + std::array shape{self.get_target_partition().get_global_size(), mu}; + py::array_t result(shape); + + self.matrix_product_global_to_global(input.data(), result.mutable_data(), mu); + + return result; + }, + py::arg("input").noconvert(true)); +} + +#endif diff --git a/src/htool/distributed_operator/utility.hpp b/src/htool/distributed_operator/utility.hpp new file mode 100644 index 0000000..13f3034 --- /dev/null +++ b/src/htool/distributed_operator/utility.hpp @@ -0,0 +1,53 @@ +#ifndef HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_UTILITY_HPP +#define HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_UTILITY_HPP + +#include "../hmatrix/interfaces/virtual_generator.hpp" +#include "distributed_operator.hpp" +#include + +template +void declare_distributed_operator_utility(py::module &m, std::string prefix = "") { + + using CustomApproximation = CustomApproximationBuilder; + using DefaultApproximation = DefaultApproximationBuilder; + using LocalDefaultApproximation = DefaultLocalApproximationBuilder; + using DistributedOperatorFromHMatrix = DistributedOperatorFromHMatrix; + + std::string custom_approximation_name = prefix + "CustomApproximationBuilder"; + std::string default_approximation_name = prefix + "DefaultApproximationBuilder"; + std::string default_local_approximation_name = prefix + "DefaultLocalApproximationBuilder"; + std::string distributed_operator_from_hmatrix_name = prefix + "DistributedOperatorFromHMatrix"; + + py::class_ custom_approximation_class(m, custom_approximation_name.c_str()); + custom_approximation_class.def(py::init &, const Cluster &, char, char, MPI_Comm_wrapper, const VirtualLocalOperator &>()); + custom_approximation_class.def_property_readonly( + "distributed_operator", [](const CustomApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); + + py::class_ default_approximation_class(m, default_approximation_name.c_str()); + default_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, char, char, MPI_Comm_wrapper>()); + default_approximation_class.def_property_readonly( + "distributed_operator", [](const DefaultApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); + default_approximation_class.def_property( + "hmatrix", [](const DefaultApproximation &self) { return &self.hmatrix; }, [](DefaultApproximation &self, const HMatrix &hmatrix) { self.hmatrix = hmatrix; }); + default_approximation_class.def_property_readonly( + "block_diagonal_hmatrix", [](const DefaultApproximation &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); + + py::class_ default_local_approximation_class(m, default_local_approximation_name.c_str()); + default_local_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, char, char, MPI_Comm_wrapper>()); + default_local_approximation_class.def_property_readonly( + "distributed_operator", [](const LocalDefaultApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); + default_local_approximation_class.def_readwrite( + "hmatrix", &LocalDefaultApproximation::hmatrix); + default_local_approximation_class.def_property_readonly( + "block_diagonal_hmatrix", [](const LocalDefaultApproximation &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); + + py::class_ distributed_operator_from_hmatrix_class(m, distributed_operator_from_hmatrix_name.c_str()); + distributed_operator_from_hmatrix_class.def(py::init &, const Cluster &, const Cluster &, const HMatrixTreeBuilder &, MPI_Comm_wrapper>()); + distributed_operator_from_hmatrix_class.def_property_readonly( + "distributed_operator", [](const DistributedOperatorFromHMatrix &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); + distributed_operator_from_hmatrix_class.def_readwrite( + "hmatrix", &DistributedOperatorFromHMatrix::hmatrix); + distributed_operator_from_hmatrix_class.def_property_readonly( + "block_diagonal_hmatrix", [](const DistributedOperatorFromHMatrix &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); +} +#endif diff --git a/src/htool/hmatrix.hpp b/src/htool/hmatrix.hpp deleted file mode 100644 index 6a2e8a1..0000000 --- a/src/htool/hmatrix.hpp +++ /dev/null @@ -1,309 +0,0 @@ -#ifndef HTOOL_HMATRIX_CPP -#define HTOOL_HMATRIX_CPP - -#include -#include -#include -#include - -#include "lrmat_generator.hpp" -#include "matrix.hpp" -#include "misc.hpp" -#include "wrapper_mpi.hpp" -#include - -namespace py = pybind11; -using namespace pybind11::literals; -using namespace htool; - -template -void declare_HMatrix(py::module &m, const std::string &baseclassName, const std::string &className) { - - py::class_>(m, baseclassName.c_str()); - - using Class = HMatrix; - py::class_> py_class(m, className.c_str()); - - // Constructor with precomputed clusters - py_class.def(py::init &, const std::shared_ptr &, double, double, char, char, const int &, MPI_Comm_wrapper>(), py::arg("cluster_target"), py::arg("cluster_source"), py::arg("epsilon") = 1e-6, py::arg("eta") = 10, py::arg("Symmetry") = 'N', py::arg("UPLO") = 'N', py::arg("reqrank") = -1, py::arg("comm") = MPI_Comm_wrapper(MPI_COMM_WORLD)); - - // Symmetric build - py_class.def("build", [](Class &self, VirtualGeneratorCpp &mat, const py::array_t &xt, const py::array_t &xs) { - self.build(mat, xt.data(), xs.data()); - }); - - py_class.def("build", [](Class &self, VirtualGeneratorCpp &mat, const py::array_t &x) { - self.build(mat, x.data()); - }); - - py_class.def("build_dense_blocks", [](Class &self, VirtualDenseBlocksGeneratorCpp &dense_block_generator) { - self.build_dense_blocks(dense_block_generator); - }); - - // Setters - py_class.def("set_maxblocksize", &Class::set_maxblocksize); - py_class.def("set_minsourcedepth", &Class::set_minsourcedepth); - py_class.def("set_mintargetdepth", &Class::set_mintargetdepth); - py_class.def("set_delay_dense_computation", &Class::set_delay_dense_computation); - py_class.def("set_compression", [](Class &self, std::shared_ptr> mat) { - self.set_compression(mat); - }); - - // Getters - py_class.def_property_readonly("shape", [](const Class &self) { - return std::array{self.nb_rows(), self.nb_cols()}; - }); - py_class.def("get_perm_t", overload_cast_<>()(&Class::get_permt, py::const_)); - py_class.def("get_perm_s", overload_cast_<>()(&Class::get_perms, py::const_)); - py_class.def("get_MasterOffset_t", overload_cast_<>()(&Class::get_MasterOffset_t, py::const_)); - py_class.def("get_MasterOffset_s", overload_cast_<>()(&Class::get_MasterOffset_s, py::const_)); - - // Linear algebra - py_class.def("__mul__", [](const Class &self, std::vector b) { - return self * b; - }); - py_class.def("matvec", [](const Class &self, std::vector b) { - return self * b; - }); - py_class.def("__matmul__", [](const Class &self, py::array_t B) { - int mu; - - if (B.ndim() == 1) { - mu = 1; - } else if (B.ndim() == 2) { - mu = B.shape()[1]; - } else { - throw std::runtime_error("Wrong dimension for HMatrix-matrix product"); // LCOV_EXCL_LINE - } - if (B.shape()[0] != self.nb_cols()) { - throw std::runtime_error("Wrong size for HMatrix-matrix product"); // LCOV_EXCL_LINE - } - - std::vector result(self.nb_rows() * mu, 0); - - self.mvprod_global_to_global(B.data(), result.data(), mu); - - if (B.ndim() == 1) { - std::array shape{self.nb_rows()}; - return py::array_t(shape, result.data()); - } else { - std::array shape{self.nb_rows(), mu}; - return py::array_t(shape, result.data()); - } - }); - - // Print information - py_class.def("print_infos", &Class::print_infos); - py_class.def("get_infos", overload_cast_()(&Class::get_infos, py::const_)); - py_class.def("__str__", [](const Class &self) { - return "HMatrix: (shape: " + htool::NbrToStr(self.nb_cols()) + "x" + htool::NbrToStr(self.nb_rows()) + ", nb_low_rank_blocks: " + htool::NbrToStr(self.get_nlrmat()) + ", nb_dense_blocks: " + htool::NbrToStr(self.get_ndmat()) + ")"; - }); - - // Plot pattern - py_class.def( - "display", [](const Class &self, bool show = true) { - const std::vector *> &lrmats = self.get_MyFarFieldMats(); - const std::vector *> &dmats = self.get_MyNearFieldMats(); - - int nb = dmats.size() + lrmats.size(); - int sizeworld = self.get_sizeworld(); - int rankworld = self.get_rankworld(); - - int nbworld[sizeworld]; - MPI_Allgather(&nb, 1, MPI_INT, nbworld, 1, MPI_INT, self.get_comm()); - int nbg = 0; - for (int i = 0; i < sizeworld; i++) { - nbg += nbworld[i]; - } - - std::vector buf(5 * nbg, 0); - - for (int i = 0; i < dmats.size(); i++) { - const SubMatrix &l = *(dmats[i]); - buf[5 * i] = l.get_offset_i(); - buf[5 * i + 1] = l.nb_rows(); - buf[5 * i + 2] = l.get_offset_j(); - buf[5 * i + 3] = l.nb_cols(); - buf[5 * i + 4] = -1; - } - - for (int i = 0; i < lrmats.size(); i++) { - const LowRankMatrix &l = *(lrmats[i]); - buf[5 * (dmats.size() + i)] = l.get_offset_i(); - buf[5 * (dmats.size() + i) + 1] = l.nb_rows(); - buf[5 * (dmats.size() + i) + 2] = l.get_offset_j(); - buf[5 * (dmats.size() + i) + 3] = l.nb_cols(); - buf[5 * (dmats.size() + i) + 4] = l.rank_of(); - } - - int displs[sizeworld]; - int recvcounts[sizeworld]; - displs[0] = 0; - - for (int i = 0; i < sizeworld; i++) { - recvcounts[i] = 5 * nbworld[i]; - if (i > 0) - displs[i] = displs[i - 1] + recvcounts[i - 1]; - } - MPI_Gatherv(rankworld == 0 ? MPI_IN_PLACE : buf.data(), recvcounts[rankworld], MPI_INT, buf.data(), recvcounts, displs, MPI_INT, 0, self.get_comm()); - - if (rankworld == 0) { - // Import - py::object plt = py::module::import("matplotlib.pyplot"); - py::object patches = py::module::import("matplotlib.patches"); - py::object colors = py::module::import("matplotlib.colors"); - py::object numpy = py::module::import("numpy"); - - // First Data - int nr = self.nb_rows(); - int nc = self.nb_cols(); - py::array_t matrix({nr, nc}); - py::array_t mask_matrix({nr, nc}); - mask_matrix.attr("fill")(false); - - // Figure - py::tuple sublots_output = plt.attr("subplots")(1, 1); - py::object fig = sublots_output[0]; - py::object axes = sublots_output[1]; - // axes.attr() - - // Issue: there a shift of one pixel along the y-axis... - // int shift = axes.transData.transform([(0,0), (1,1)]) - // shift = shift[1,1] - shift[0,1] # 1 unit in display coords - int shift = 0; - - int max_rank = 0; - for (int p = 0; p < nbg; p++) { - int i_row = buf[5 * p]; - int nb_row = buf[5 * p + 1]; - int i_col = buf[5 * p + 2]; - int nb_col = buf[5 * p + 3]; - int rank = buf[5 * p + 4]; - - if (rank > max_rank) { - max_rank = rank; - } - for (int i = 0; i < nb_row; i++) { - for (int j = 0; j < nb_col; j++) { - matrix.mutable_at(i_row + i, i_col + j) = rank; - if (rank == -1) { - mask_matrix.mutable_at(i_row + i, i_col + j) = true; - } - } - } - - py::object rect = patches.attr("Rectangle")(py::make_tuple(i_col - 0.5, i_row - 0.5 + shift), nb_col, nb_row, "linewidth"_a = 0.75, "edgecolor"_a = 'k', "facecolor"_a = "none"); - axes.attr("add_patch")(rect); - - if (rank >= 0 && nb_col / double(nc) > 0.05 && nb_row / double(nc) > 0.05) { - axes.attr("annotate")(rank, py::make_tuple(i_col + nb_col / 2., i_row + nb_row / 2.), "color"_a = "white", "size"_a = 10, "va"_a = "center", "ha"_a = "center"); - } - } - - // Colormap - py::object cmap = plt.attr("get_cmap")("YlGn"); - py::object new_cmap = colors.attr("LinearSegmentedColormap").attr("from_list")("trunc(YlGn,0.4,1)", cmap(numpy.attr("linspace")(0.4, 1, 100))); - - // Plot - py::object masked_matrix = numpy.attr("ma").attr("array")(matrix, "mask"_a = mask_matrix); - new_cmap.attr("set_bad")("color"_a = "red"); - - plt.attr("imshow")(masked_matrix, "cmap"_a = new_cmap, "vmin"_a = 0, "vmax"_a = 10); - plt.attr("draw")(); - if (show) { - plt.attr("show")(); // LCOV_EXCL_LINE - } - } - }, - py::arg("show") = true); - - // Plot clustering - - py_class.def( - "display_cluster", [](const Class &self, py::array_t points_target, int depth, std::string type, bool show) { - int sizeworld = self.get_sizeworld(); - int rankworld = self.get_rankworld(); - - if (rankworld == 0) { - - VirtualCluster const *root; - if (type == "target") { - root = (self.get_target_cluster()); - } else if (type == "source") { - root = (self.get_source_cluster()); - } else { - std::cout << "Choose between target and source" << std::endl; // LCOV_EXCL_LINE - return 0; // LCOV_EXCL_LINE - } - - std::stack s; - s.push(root); - - int size = root->get_size(); - int space_dim = root->get_space_dim(); - std::vector output((space_dim + 1) * size); - - // Permuted geometric points - for (int i = 0; i < size; ++i) { - for (int p = 0; p < space_dim; p++) { - output[i + size * p] = points_target.at(p, root->get_perm(i)); - } - } - - int counter = 0; - while (!s.empty()) { - VirtualCluster const *curr = s.top(); - s.pop(); - - if (depth == curr->get_depth()) { - std::fill_n(&(output[space_dim * size + curr->get_offset()]), curr->get_size(), counter); - counter += 1; - } - - // Recursion - if (!curr->IsLeaf()) { - - for (int p = 0; p < curr->get_nb_sons(); p++) { - s.push(&(curr->get_son(p))); - } - } - } - - // Import - py::object plt = py::module::import("matplotlib.pyplot"); - py::object colors = py::module::import("matplotlib.colors"); - - // Create Color Map - py::object colormap = plt.attr("get_cmap")("Dark2"); - py::object norm = colors.attr("Normalize")("vmin"_a = (*std::min_element(output.begin() + space_dim * size, output.end())), "vmax"_a = (*std::max_element(output.begin() + space_dim * size, output.end()))); - - // Figure - py::object fig = plt.attr("figure")(); - - if (space_dim == 2) { - py::object ax = fig.attr("add_subplot")(111); - ax.attr("scatter")(std::vector(output.begin(), output.begin() + size), std::vector(output.begin() + size, output.begin() + 2 * size), "c"_a = colormap(norm(std::vector(output.begin() + 2 * size, output.end()))), "marker"_a = 'o'); - - } else if (space_dim == 3) { - py::object Axes3D = py::module::import("mpl_toolkits.mplot3d").attr("Axes3D"); - - py::object ax = fig.attr("add_subplot")(111, "projection"_a = "3d"); - ax.attr("scatter")(std::vector(output.begin(), output.begin() + size), std::vector(output.begin() + size, output.begin() + 2 * size), std::vector(output.begin() + 2 * size, output.begin() + 3 * size), "c"_a = colormap(norm(std::vector(output.begin() + 3 * size, output.end()))), "marker"_a = 'o'); - } - - plt.attr("draw")(); - if (show) { - plt.attr("show")(); // LCOV_EXCL_LINE - } - return 0; - } - - return 0; - }, - py::arg("points_target"), - py::arg("depth"), - py::arg("type") = "target", - py::arg("show") = true); -} - -#endif diff --git a/src/htool/hmatrix/hmatrix.hpp b/src/htool/hmatrix/hmatrix.hpp new file mode 100644 index 0000000..d5237cb --- /dev/null +++ b/src/htool/hmatrix/hmatrix.hpp @@ -0,0 +1,55 @@ +#ifndef HTOOL_HMATRIX_CPP +#define HTOOL_HMATRIX_CPP + +#include +#include +#include +#include +#include + +#include "../misc/wrapper_mpi.hpp" + +#include +#include +#include +#include + +namespace py = pybind11; +using namespace pybind11::literals; +using namespace htool; + +template +void declare_HMatrix(py::module &m, const std::string &className) { + + using Class = HMatrix; + py::class_ py_class(m, className.c_str()); + + py_class.def("to_dense", [](const Class &self) { + std::array shape{self.get_target_cluster().get_size(), self.get_source_cluster().get_size()}; + py::array_t dense(shape); + std::fill_n(dense.mutable_data(), dense.size(), CoefficientPrecision(0)); + copy_to_dense(self, dense.mutable_data()); + return dense; + }); + + py_class.def("to_dense_in_user_numbering", [](const Class &self) { + std::array shape{self.get_target_cluster().get_size(), self.get_source_cluster().get_size()}; + py::array_t dense(shape); + std::fill_n(dense.mutable_data(), dense.size(), CoefficientPrecision(0)); + copy_to_dense_in_user_numbering(self, dense.mutable_data()); + return dense; + }); + + py_class.def("__deepcopy__", [](const Class &self, py::dict) { return Class(self); }, "memo"_a); + + py_class.def("get_tree_parameters", [](const HMatrix &hmatrix) { return htool::get_tree_parameters(hmatrix); }); + py_class.def("get_local_information", [](const HMatrix &hmatrix) { return htool::get_hmatrix_information(hmatrix); }); + py_class.def("get_distributed_information", [](const HMatrix &hmatrix, MPI_Comm_wrapper comm) { return htool::get_distributed_hmatrix_information(hmatrix, comm); }); + + m.def("recompression", &htool::recompression &)>>); + m.def("recompression", [](HMatrix &hmatrix) { recompression(hmatrix); }); + m.def("openmp_recompression", &htool::openmp_recompression &)>>); + m.def("openmp_recompression", [](HMatrix &hmatrix) { recompression(hmatrix); }); +} + +#endif diff --git a/src/htool/hmatrix/hmatrix_builder.hpp b/src/htool/hmatrix/hmatrix_builder.hpp new file mode 100644 index 0000000..3f60c29 --- /dev/null +++ b/src/htool/hmatrix/hmatrix_builder.hpp @@ -0,0 +1,43 @@ +#ifndef HTOOL_HMATRIX_BUILDER_CPP +#define HTOOL_HMATRIX_BUILDER_CPP + +#include "interfaces/virtual_dense_blocks_generator.hpp" +#include "interfaces/virtual_generator.hpp" +#include "interfaces/virtual_low_rank_generator.hpp" +#include + +template +void declare_hmatrix_builder(py::module &m, const std::string &className) { + using Class = htool::HMatrixTreeBuilder; + py::class_ py_class(m, className.c_str()); + + // // Constructor + // py_class.def(py::init([](underlying_type epsilon, CoordinatePrecision eta, char symmetry, char UPLO) { + // return std::unique_ptr(new Class(epsilon, eta, symmetry, UPLO)); + // }), + // py::arg("epsilon"), + // py::arg("eta"), + // py::arg("symmetry"), + // py::arg("UPLO")); + + py_class.def(py::init([](underlying_type epsilon, CoordinatePrecision eta, char symmetry, char UPLO, int reqrank, std::shared_ptr> low_rank_strategy) { + return std::unique_ptr(new Class(epsilon, eta, symmetry, UPLO, reqrank, low_rank_strategy)); + }), + py::arg("epsilon"), + py::arg("eta"), + py::arg("symmetry"), + py::arg("UPLO"), + py::kw_only(), + py::arg("reqrank") = -1, + py::arg("low_rank_strategy") = nullptr); + + // Build + py_class.def("build", [](const Class &self, const VirtualGenerator &generator, const Cluster &target_cluster, const Cluster &source_cluster, int target_partition_number, int partition_number_for_symmetry) { return self.build(generator, target_cluster, source_cluster, target_partition_number, partition_number_for_symmetry); }, py::arg("generator"), py::arg("target_cluster"), py::arg("source_cluster"), py::arg("target_partition_number") = -1, py::arg("partition_number_for_symmetry") = -1); + + // Setters + py_class.def("set_minimal_source_depth", &Class::set_minimal_source_depth); + py_class.def("set_minimal_target_depth", &Class::set_minimal_target_depth); + py_class.def("set_low_rank_generator", [](Class &self, std::shared_ptr> low_rank_generator) { self.set_low_rank_generator(low_rank_generator); }); + py_class.def("set_dense_blocks_generator", [](Class &self, std::shared_ptr> dense_blocks_generator) { self.set_dense_blocks_generator(dense_blocks_generator); }); +} +#endif diff --git a/src/htool/hmatrix/interfaces/virtual_dense_blocks_generator.hpp b/src/htool/hmatrix/interfaces/virtual_dense_blocks_generator.hpp new file mode 100644 index 0000000..d544e89 --- /dev/null +++ b/src/htool/hmatrix/interfaces/virtual_dense_blocks_generator.hpp @@ -0,0 +1,71 @@ +#ifndef HTOOL_DENSE_BLOCKS_GENERATOR_CPP +#define HTOOL_DENSE_BLOCKS_GENERATOR_CPP + +#include + +namespace py = pybind11; +using namespace pybind11::literals; +using namespace htool; + +template > +class VirtualDenseBlocksGeneratorPython : public VirtualDenseBlocksGenerator { + + const Cluster &m_target_cluster; + const Cluster &m_source_cluster; + + public: + using VirtualDenseBlocksGenerator::VirtualDenseBlocksGenerator; + + VirtualDenseBlocksGeneratorPython(const Cluster &target_cluster, const Cluster &source_cluster) : VirtualDenseBlocksGenerator(), m_target_cluster(target_cluster), m_source_cluster(source_cluster) {} + + void copy_dense_blocks(const std::vector &M, const std::vector &N, const std::vector &rows_offsets, const std::vector &cols_offsets, std::vector &ptr) const override { + int nb_blocks = M.size(); + auto &target_permutation = m_target_cluster.get_permutation(); + auto &source_permutation = m_source_cluster.get_permutation(); + std::vector> vec_ptr; + std::vector> rows_ptr; + std::vector> cols_ptr; + for (int i = 0; i < nb_blocks; i++) { + rows_ptr.emplace_back(std::array{M[i]}, target_permutation.data() + rows_offsets[i], py::capsule(target_permutation.data() + rows_offsets[i])); + cols_ptr.emplace_back(std::array{N[i]}, source_permutation.data() + cols_offsets[i], py::capsule(source_permutation.data() + cols_offsets[i])); + vec_ptr.emplace_back(std::array{M[i], N[i]}, ptr[i], py::capsule(ptr[i])); + } + + build_dense_blocks(rows_ptr, cols_ptr, vec_ptr); + } + + // lcov does not see it because of trampoline I assume + virtual void build_dense_blocks(const std::vector> &rows, const std::vector> &cols, std::vector> &blocks) const = 0; // LCOV_EXCL_LINE +}; + +template +class PyVirtualDenseBlocksGenerator : public VirtualDenseBlocksGeneratorPython { + public: + using VirtualDenseBlocksGeneratorPython::VirtualDenseBlocksGeneratorPython; + // PyVirtualGenerator(int nr0, int nc0): IMatrix(nr0,nc0){} + + /* Trampoline (need one for each virtual function) */ + virtual void build_dense_blocks(const std::vector> &rows, const std::vector> &cols, std::vector> &blocks) const override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + VirtualDenseBlocksGeneratorPython, /* Parent class */ + build_dense_blocks, /* Name of function in C++ (must match Python name) */ + rows, + cols, + blocks /* Argument(s) */ + ); + } +}; + +template > +void declare_custom_VirtualDenseBlocksGenerator(py::module &m, const std::string &className) { + // using BaseClass = VirtualDenseBlocksGenerator; + // py::class_>(m, base_class_name.c_str()); + + using Class = VirtualDenseBlocksGeneratorPython; + py::class_, PyVirtualDenseBlocksGenerator> py_class(m, className.c_str()); + py_class.def(py::init &, const Cluster &>()); + py_class.def("build_dense_blocks", &Class::build_dense_blocks); +} + +#endif diff --git a/src/htool/hmatrix/interfaces/virtual_generator.hpp b/src/htool/hmatrix/interfaces/virtual_generator.hpp new file mode 100644 index 0000000..532c160 --- /dev/null +++ b/src/htool/hmatrix/interfaces/virtual_generator.hpp @@ -0,0 +1,60 @@ +#ifndef HTOOL_GENERATOR_CPP +#define HTOOL_GENERATOR_CPP + +#include +#include + +using namespace htool; + +template +class VirtualGeneratorPython : public htool::VirtualGenerator { + public: + using VirtualGenerator::VirtualGenerator; + + VirtualGeneratorPython(const py::array_t &target_permutation, const py::array_t &source_permutation) : VirtualGenerator() {} + + void copy_submatrix(int M, int N, const int *const rows, const int *const cols, CoefficientPrecision *ptr) const override { + if (M * N > 0) { + py::array_t mat(std::array{M, N}, ptr, py::capsule(ptr)); + + py::array_t py_rows(std::array{M}, rows, py::capsule(rows)); + py::array_t py_cols(std::array{N}, cols, py::capsule(cols)); + + build_submatrix(py_rows, py_cols, mat); + } + } + + // lcov does not see it because of trampoline I assume + virtual void build_submatrix(const py::array_t &J, const py::array_t &K, py::array_t &mat) const = 0; // LCOV_EXCL_LINE +}; + +template +class PyVirtualGenerator : public VirtualGeneratorPython { + public: + using VirtualGeneratorPython::VirtualGeneratorPython; + + /* Trampoline (need one for each virtual function) */ + virtual void build_submatrix(const py::array_t &J, const py::array_t &K, py::array_t &mat) const override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + VirtualGeneratorPython, /* Parent class */ + build_submatrix, /* Name of function in C++ (must match Python name) */ + J, + K, + mat /* Argument(s) */ + ); + } +}; + +template +void declare_virtual_generator(py::module &m, const std::string &className, const std::string &base_class_name) { + using BaseClass = VirtualGenerator; + py::class_(m, base_class_name.c_str()); + + using Class = VirtualGeneratorPython; + py::class_> py_class(m, className.c_str()); + py_class.def(py::init<>()); + py_class.def("build_submatrix", &Class::build_submatrix); +} + +#endif diff --git a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp new file mode 100644 index 0000000..62f3b3a --- /dev/null +++ b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp @@ -0,0 +1,85 @@ +#ifndef HTOOL_LRMAT_GENERATOR_CPP +#define HTOOL_LRMAT_GENERATOR_CPP + +#include +#include +#include +#include +#include + +namespace py = pybind11; +using namespace pybind11::literals; +using namespace htool; + +template +class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator { + mutable std::vector> m_mats_U; // owned by Python + mutable std::vector> m_mats_V; // owned by Python + + public: + using VirtualLowRankGenerator::VirtualLowRankGenerator; + + VirtualLowRankGeneratorPython() {} + + bool copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, LowRankMatrix &lrmat) const override { + auto &U = lrmat.get_U(); + auto &V = lrmat.get_V(); + py::array_t py_rows(std::array{M}, rows, py::capsule(rows)); + py::array_t py_cols(std::array{N}, cols, py::capsule(cols)); + + bool success = build_low_rank_approximation(py_rows, py_cols, lrmat.get_epsilon()); + if (success) { + U.assign(m_mats_U.back().shape()[0], m_mats_U.back().shape()[1], m_mats_U.back().mutable_data(), false); + V.assign(m_mats_V.back().shape()[0], m_mats_V.back().shape()[1], m_mats_V.back().mutable_data(), false); + } + return success; + } + + bool copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, int reqrank, LowRankMatrix &lrmat) const override { + Logger::get_instance().log(LogLevel::ERROR, "copy_low_rank_approximation with required rank is not implemented in the python interface."); + return false; + } + + // lcov does not see it because of trampoline I assume + virtual bool build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const = 0; // LCOV_EXCL_LINE + + 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 clear_data() { + m_mats_U.clear(); + m_mats_V.clear(); + } +}; + +template +class PyVirtualLowRankGenerator : public VirtualLowRankGeneratorPython { + public: + using VirtualLowRankGeneratorPython::VirtualLowRankGeneratorPython; + + /* Trampoline (need one for each virtual function) */ + virtual bool build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const override { + PYBIND11_OVERRIDE_PURE( + bool, /* Return type */ + VirtualLowRankGeneratorPython, /* Parent class */ + build_low_rank_approximation, /* Name of function in C++ (must match Python name) */ + rows, + cols, + epsilon); + } +}; + +template +void declare_custom_VirtualLowRankGenerator(py::module &m, const std::string &className) { + using Class = VirtualLowRankGeneratorPython; + py::class_, PyVirtualLowRankGenerator> py_class(m, className.c_str()); + py_class.def(py::init<>()); + py_class.def("build_low_rank_approximation", &Class::build_low_rank_approximation); + py_class.def("set_U", &Class::set_U); + py_class.def("set_V", &Class::set_V); + py_class.def("clear_data", &Class::clear_data); +} + +#endif diff --git a/src/htool/hmatrix/lrmat.hpp b/src/htool/hmatrix/lrmat.hpp new file mode 100644 index 0000000..a074da3 --- /dev/null +++ b/src/htool/hmatrix/lrmat.hpp @@ -0,0 +1,19 @@ +#ifndef HTOOL_LRMAT_CPP +#define HTOOL_LRMAT_CPP + +#include +#include + +namespace py = pybind11; +using namespace htool; + +template +void declare_LowRankMatrix(py::module &m, const std::string &className) { + using Class = LowRankMatrix; + py::class_ py_class(m, className.c_str()); + + py_class.def("nb_rows", &Class::nb_rows); + py_class.def("nb_cols", &Class::nb_cols); + py_class.def("rank", &Class::rank_of, "test"); +} +#endif diff --git a/src/htool/local_operator/local_dense_operator.hpp b/src/htool/local_operator/local_dense_operator.hpp new file mode 100644 index 0000000..3d6a736 --- /dev/null +++ b/src/htool/local_operator/local_dense_operator.hpp @@ -0,0 +1,15 @@ +#ifndef HTOOL_LOCAL_DENSE_OPERATOR_CPP +#define HTOOL_LOCAL_DENSE_OPERATOR_CPP + +#include "local_operator.hpp" +#include +#include + +template +void declare_local_dense_matrix(py::module &m, const std::string &class_name) { + using Class = htool::LocalDenseMatrix; + py::class_> py_data(m, class_name.c_str()); + py_data.def(py::init &, const Cluster &, const Cluster &, char, char, bool, bool, int>()); +} + +#endif diff --git a/src/htool/local_operator/local_hmatrix.hpp b/src/htool/local_operator/local_hmatrix.hpp new file mode 100644 index 0000000..a50b2c1 --- /dev/null +++ b/src/htool/local_operator/local_hmatrix.hpp @@ -0,0 +1,47 @@ +#ifndef HTOOL_LOCAL_HMATRIX_CPP +#define HTOOL_LOCAL_HMATRIX_CPP + +#include "../misc/wrapper_mpi.hpp" +#include "local_operator.hpp" +#include +#include +#include + +template +void declare_local_hmatrix(py::module &m, const std::string &class_name) { + using Class = htool::LocalHMatrix; + py::class_> py_data(m, class_name.c_str()); + + py_data.def(py::init &, const Cluster &, const Cluster &, char, char, bool, bool>(), py::keep_alive<1, 2>(), py::keep_alive<1, 3>(), py::keep_alive<1, 4>()); + + // py_data.def("print_information", [](const Class &self, MPI_Comm_wrapper comm) { + // htool::print_distributed_hmatrix_information(self.get_hmatrix(), std::cout, comm); + // }); + + // py_data.def("display", [](const Class &self) { + // std::vector> output_blocks{}; + // const HMatrix &hmatrix = self.get_hmatrix(); + // preorder_tree_traversal( + // hmatrix, + // [&output_blocks, &hmatrix](const HMatrix ¤t_hmatrix) { + // if (current_hmatrix.is_leaf()) { + // output_blocks.push_back(DisplayBlock{current_hmatrix.get_target_cluster().get_offset() - hmatrix.get_target_cluster().get_offset(), current_hmatrix.get_source_cluster().get_offset() - hmatrix.get_source_cluster().get_offset(), current_hmatrix.get_target_cluster().get_size(), current_hmatrix.get_source_cluster().get_size(), current_hmatrix.get_rank()}); + // } + // }); + + // // Import + // py::object plt = py::module::import("matplotlib.pyplot"); + // py::object patches = py::module::import("matplotlib.patches"); + // py::object colors = py::module::import("matplotlib.colors"); + // py::object numpy = py::module::import("numpy"); + + // // First Data + // int nr = self.nb_rows(); + // int nc = self.nb_cols(); + // py::array_t matrix({nr, nc}); + // py::array_t mask_matrix({nr, nc}); + // mask_matrix.attr("fill")(false); + // }); +} + +#endif diff --git a/src/htool/local_operator/local_operator.hpp b/src/htool/local_operator/local_operator.hpp new file mode 100644 index 0000000..3304054 --- /dev/null +++ b/src/htool/local_operator/local_operator.hpp @@ -0,0 +1,101 @@ +#ifndef HTOOL_LOCAL_OPERATOR_CPP +#define HTOOL_LOCAL_OPERATOR_CPP + +#include +#include +#include + +template +class LocalOperatorPython : public htool::LocalOperator { + public: + using htool::LocalOperator::LocalOperator; + + LocalOperatorPython(const Cluster &cluster_tree_target, const Cluster &cluster_tree_source, char symmetry = 'N', char UPLO = 'N', bool target_use_permutation_to_mvprod = false, bool source_use_permutation_to_mvprod = false) : LocalOperator(cluster_tree_target, cluster_tree_source, symmetry, UPLO, target_use_permutation_to_mvprod, source_use_permutation_to_mvprod) {} + + void local_add_vector_product(char trans, CoefficientPrecision alpha, const CoefficientPrecision *in, CoefficientPrecision beta, CoefficientPrecision *out) const override { + + py::array_t input(std::array{this->m_source_cluster.get_size()}, in, py::capsule(in)); + py::array_t output(std::array{this->m_target_cluster.get_size()}, out, py::capsule(out)); + + add_vector_product(trans, alpha, input, beta, output); + } + + void local_add_vector_product_symmetric(char trans, CoefficientPrecision alpha, const CoefficientPrecision *in, CoefficientPrecision beta, CoefficientPrecision *out, char UPLO, char symmetry) const override { + + py::array_t input(std::array{this->m_source_cluster.get_size()}, in, py::capsule(in)); + py::array_t output(std::array{this->m_target_cluster.get_size()}, out, py::capsule(out)); + + add_vector_product(trans, alpha, input, beta, output); + } + + void local_add_matrix_product_row_major(char trans, CoefficientPrecision alpha, const CoefficientPrecision *in, CoefficientPrecision beta, CoefficientPrecision *out, int mu) const override { + + py::array_t input(std::array{this->m_source_cluster.get_size(), mu}, in, py::capsule(in)); + py::array_t output(std::array{this->m_target_cluster.get_size(), mu}, out, py::capsule(out)); + + add_matrix_product_row_major(trans, alpha, input, beta, output); + } + + void local_add_matrix_product_symmetric_row_major(char trans, CoefficientPrecision alpha, const CoefficientPrecision *in, CoefficientPrecision beta, CoefficientPrecision *out, int mu, char UPLO, char symmetry) const override { + + py::array_t input(std::array{this->m_source_cluster.get_size(), 1}, in, py::capsule(in)); + py::array_t output(std::array{this->m_target_cluster.get_size(), 1}, out, py::capsule(out)); + + add_matrix_product_row_major(trans, alpha, input, beta, output); + } + + // lcov does not see it because of trampoline I assume + virtual void add_vector_product(char trans, CoefficientPrecision alpha, const py::array_t &in, CoefficientPrecision beta, py::array_t &out) const = 0; // LCOV_EXCL_LINE + // virtual void local_add_vector_product_symmetric(char trans, CoefficientPrecision alpha, const std::vector &in, CoefficientPrecision beta, std::vector &out) const = 0; // LCOV_EXCL_LINE + + virtual void add_matrix_product_row_major(char trans, CoefficientPrecision alpha, const py::array_t &in, CoefficientPrecision beta, py::array_t &out) const = 0; // LCOV_EXCL_LINE +}; + +template +class PyLocalOperator : public LocalOperatorPython { + public: + using LocalOperatorPython::LocalOperatorPython; + + /* Trampoline (need one for each virtual function) */ + virtual void add_vector_product(char trans, CoefficientPrecision alpha, const py::array_t &in, CoefficientPrecision beta, py::array_t &out) const override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + LocalOperatorPython, /* Parent class */ + add_vector_product, /* Name of function in C++ (must match Python name) */ + trans, + alpha, + in, + beta, + out /* Argument(s) */ + ); + } + virtual void add_matrix_product_row_major(char trans, CoefficientPrecision alpha, const py::array_t &in, CoefficientPrecision beta, py::array_t &out) const override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + LocalOperatorPython, /* Parent class */ + add_matrix_product_row_major, /* Name of function in C++ (must match Python name) */ + trans, + alpha, + in, + beta, + out /* Argument(s) */ + ); + } +}; + +template +void declare_local_operator(py::module &m, const std::string &class_name) { + using VirtualClass = htool::VirtualLocalOperator; + py::class_(m, ("Virtual" + class_name).c_str()); + + using BaseClass = LocalOperator; + py::class_ py_base_class(m, ("Base" + class_name).c_str()); + + using Class = LocalOperatorPython; + py::class_, BaseClass> py_class(m, class_name.c_str()); + py_class.def(py::init &, const Cluster &, char, char, bool, bool>()); + py_class.def("add_vector_product", &Class::add_vector_product, py::arg("trans"), py::arg("alpha"), py::arg("in").noconvert(true), py::arg("beta"), py::arg("out").noconvert(true)); + py_class.def("add_matrix_product_row_major", &Class::add_matrix_product_row_major); +} + +#endif diff --git a/src/htool/local_operator/virtual_local_operator.hpp b/src/htool/local_operator/virtual_local_operator.hpp new file mode 100644 index 0000000..c2c9e0a --- /dev/null +++ b/src/htool/local_operator/virtual_local_operator.hpp @@ -0,0 +1,18 @@ +#ifndef HTOOL_VIRTUAL_LOCAL_OPERATOR_CPP +#define HTOOL_VIRTUAL_LOCAL_OPERATOR_CPP + +#include +#include + +// template +// class PyVirtualLocalOperator : public htool::VirtualLocalOperator { +// public: +// }; + +template +void declare_interface_local_operator(py::module &m, const std::string &class_name) { + using Class = htool::VirtualLocalOperator; + py::class_(m, class_name.c_str()); +} + +#endif diff --git a/src/htool/lrmat_generator.hpp b/src/htool/lrmat_generator.hpp deleted file mode 100644 index c34a07d..0000000 --- a/src/htool/lrmat_generator.hpp +++ /dev/null @@ -1,88 +0,0 @@ -#ifndef HTOOL_LRMAT_GENERATOR_CPP -#define HTOOL_LRMAT_GENERATOR_CPP - -#include -#include -#include -#include -#include -#include - -namespace py = pybind11; -using namespace pybind11::literals; -using namespace htool; - -template -class VirtualLowRankGeneratorCpp : public VirtualLowRankGenerator { - - py::array_t mat_U, mat_V; - int rank; - - public: - using VirtualLowRankGenerator::VirtualLowRankGenerator; - - void copy_low_rank_approximation(double epsilon, int M, int N, const int *const rows, const int *const cols, int &rank0, T **U, T **V, const VirtualGenerator &A, const VirtualCluster &t, const double *const xt, const VirtualCluster &s, const double *const xs) const override { - - build_low_rank_approximation(epsilon, rank0, A, std::vector(rows, rows + M), std::vector(cols, cols + N)); - *U = new T[M * rank]; - *V = new T[N * rank]; - rank0 = rank; - std::copy_n(mat_U.data(), mat_U.size(), *U); - std::copy_n(mat_V.data(), mat_V.size(), *V); - } - - // lcov does not see it because of trampoline I assume - virtual void build_low_rank_approximation(double epsilon, int rank, const VirtualGenerator &A, const std::vector &J, const std::vector &K) const = 0; // LCOV_EXCL_LINE - - void set_U(py::array_t U0) { mat_U = U0; } - void set_V(py::array_t V0) { mat_V = V0; } - void set_rank(int rank0) { rank = rank0; } -}; - -template -class PyVirtualLowRankGenerator : public VirtualLowRankGeneratorCpp { - public: - using VirtualLowRankGeneratorCpp::VirtualLowRankGeneratorCpp; - // PyVirtualGenerator(int nr0, int nc0): IMatrix(nr0,nc0){} - - /* Trampoline (need one for each virtual function) */ - virtual void build_low_rank_approximation(double epsilon, int rank, const VirtualGenerator &A, const std::vector &J, const std::vector &K) const override { - PYBIND11_OVERRIDE_PURE( - void, /* Return type */ - VirtualLowRankGeneratorCpp, /* Parent class */ - build_low_rank_approximation, /* Name of function in C++ (must match Python name) */ - epsilon, - rank, - A, - J, - K /* Argument(s) */ - ); - } -}; - -template -void declare_custom_VirtualLowRankGenerator(py::module &m, const std::string &className) { - using Class = VirtualLowRankGeneratorCpp; - py::class_, PyVirtualLowRankGenerator> py_class(m, className.c_str()); - py_class.def(py::init<>()); - py_class.def("build_low_rank_approximation", &Class::build_low_rank_approximation); - py_class.def("set_U", &Class::set_U); - py_class.def("set_V", &Class::set_V); - py_class.def("set_rank", &Class::set_rank); -} - -// template -// void declare_VirtualLowRankGenerator(py::module &m, const std::string &className) { -// using Class = VirtualLowRankGenerator; -// py::class_ py_class(m, className.c_str()); -// } - -// template