From 232f22318c319bdfcba8224cb7a2b629a2a6a7ac Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sat, 23 Nov 2024 20:36:37 +0900 Subject: [PATCH] Option to run without linking BLAS and LAPACK in C --- .../phono3py-pytest-conda-nolapacke.yml | 49 + CMakeLists.txt | 207 +- c/_phono3py.c | 1877 ----------------- c/_phono3py.cpp | 32 +- c/_phononcalc.c | 232 -- c/lapack_wrapper.c | 33 +- c/lapack_wrapper.h | 29 +- c/phono3py.c | 16 +- c/phono3py.h | 5 +- c/phonon.c | 39 +- c/reciprocal_to_normal.c | 2 +- phono3py/cui/phono3py_script.py | 2 + phono3py/phonon/solver.py | 13 + pyproject.toml | 1 + test/phonon3/test_fc3.py | 68 +- 15 files changed, 320 insertions(+), 2285 deletions(-) create mode 100644 .github/workflows/phono3py-pytest-conda-nolapacke.yml delete mode 100644 c/_phono3py.c delete mode 100644 c/_phononcalc.c diff --git a/.github/workflows/phono3py-pytest-conda-nolapacke.yml b/.github/workflows/phono3py-pytest-conda-nolapacke.yml new file mode 100644 index 00000000..2ea0756f --- /dev/null +++ b/.github/workflows/phono3py-pytest-conda-nolapacke.yml @@ -0,0 +1,49 @@ +name: Pytest without linking BLAS and LAPACK in C + +on: + pull_request: + branches: [ develop ] + +jobs: + build-linux: + runs-on: ubuntu-latest + defaults: + run: + shell: bash -l {0} + strategy: + matrix: + python-version: ["3.12"] + + steps: + - uses: actions/checkout@v4 + # Use conda-incubator/setup-miniconda for precise control of conda infrastructure + - uses: conda-incubator/setup-miniconda@v3 + with: + miniforge-version: latest + - name: Install dependent packages + run: | + conda activate test + conda install --yes python=${{ matrix.python-version }} + conda install --yes matplotlib-base pyyaml h5py scipy pytest spglib alm cmake c-compiler cxx-compiler pypolymlp + - name: Install symfc develop branch + run: | + conda activate test + git clone https://github.com/symfc/symfc.git + cd symfc + pip install -e . -vvv + cd .. + - name: Install phonopy develop branch + run: | + conda activate test + git clone https://github.com/phonopy/phonopy.git + cd phonopy + pip install -e . -vvv + cd .. + - name: Install phono3py + run: | + conda activate test + BUILD_WITHOUT_LAPACKE=ON pip install -e . -vvv + - name: Run pytest + run: | + conda activate test + pytest -v test diff --git a/CMakeLists.txt b/CMakeLists.txt index a4d236e2..11dc20d1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,7 @@ option(PHONO3PY_USE_OMP "Option to search OpenMP library" ON) option(PHONO3PY_USE_MTBLAS "Use multithread BLAS if it exists" ON) option(PHONO3PY_WITH_TESTS "build unit tests" OFF) option(BUILD_SHARED_LIBS "Option to build shared library" OFF) +option(BUILD_WITHOUT_LAPACKE "Option to build without LAPACKE" ON) if(PHONO3PY_WITH_Fortran) enable_language(Fortran) @@ -89,7 +90,8 @@ endif() if(BUILD_PHPHCALC_LIB OR BUILD_PHONONCALC_LIB - OR BUILD_NANOBIND_MODULE) + OR BUILD_NANOBIND_MODULE + AND (NOT BUILD_WITHOUT_LAPACKE)) find_package(BLAS REQUIRED) # set BLAS_LIBRARIES if(BLAS_FOUND) @@ -105,7 +107,7 @@ if(BUILD_PHPHCALC_LIB endif() if(BLAS_LIBRARIES MATCHES "libmkl") - message(STATUS "MKL detected: Set C-macro MKL_LAPACKE.") + message(STATUS "MKL detected.") if(PHONO3PY_USE_MTBLAS) message( @@ -190,68 +192,95 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) # Shared library add_library(phphcalc_lib SHARED ${SOURCES_PHPHCALC}) - if(OpenMP_FOUND) - target_link_libraries( - phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) - else() - target_link_libraries(phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS - LAPACK::LAPACK) - endif() - - target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) - - if(BLAS_LIBRARIES MATCHES "libmkl") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions( - phphcalc_lib PRIVATE MKL_LAPACKE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib + OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE MKL_LAPACKE - THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib) endif() - endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions(phphcalc_lib PRIVATE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + target_compile_definitions(phphcalc_lib PRIVATE NO_INCLUDE_LAPACKE + THM_EPSILON=1e-10) + else() + if(OpenMP_FOUND) + target_link_libraries( + phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS + LAPACK::LAPACK) endif() - endif() - else() - # Static link library - add_library(phphcalc_lib STATIC ${SOURCES_PHPHCALC}) - if(OpenMP_FOUND) - target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) - else() - target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK) + if(BLAS_LIBRARIES MATCHES "libmkl") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MKL_BLAS MULTITHREADED_BLAS + THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE MKL_BLAS THM_EPSILON=1e-10) + endif() + endif() + + if(BLAS_LIBRARIES MATCHES "libopenblas") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MULTITHREADED_BLAS THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE THM_EPSILON=1e-10) + endif() + endif() endif() target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) - if(BLAS_LIBRARIES MATCHES "libmkl") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions( - phphcalc_lib PRIVATE MKL_LAPACKE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + else() + # Static link library + add_library(phphcalc_lib STATIC ${SOURCES_PHPHCALC}) + + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phphcalc_lib recgrid_lib OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE MKL_LAPACKE - THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib recgrid_lib) endif() - endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions(phphcalc_lib PRIVATE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + target_compile_definitions(phphcalc_lib PRIVATE NO_INCLUDE_LAPACKE + THM_EPSILON=1e-10) + else() + if(OpenMP_FOUND) + target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS + LAPACK::LAPACK OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS + LAPACK::LAPACK) + endif() + + if(BLAS_LIBRARIES MATCHES "libmkl") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MKL_BLAS MULTITHREADED_BLAS + THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE MKL_BLAS THM_EPSILON=1e-10) + endif() + endif() + + if(BLAS_LIBRARIES MATCHES "libopenblas") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MULTITHREADED_BLAS THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE THM_EPSILON=1e-10) + endif() endif() endif() + + target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) endif() if(NOT BUILD_NANOBIND_MODULE) @@ -285,45 +314,67 @@ if(BUILD_PHONONCALC_LIB OR BUILD_NANOBIND_MODULE) # Shared library add_library(phononcalc_lib SHARED ${SOURCES_PHONONCALC}) - if(OpenMP_FOUND) - target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phononcalc_lib OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib) + endif() + + target_compile_definitions(phononcalc_lib PRIVATE NO_INCLUDE_LAPACKE) else() - target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK) - endif() + if(OpenMP_FOUND) + target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK) + endif() - target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) + if(BLAS_LIBRARIES MATCHES "libmkl") + target_compile_definitions(phononcalc_lib PRIVATE MKL_BLAS + MULTITHREADED_BLAS) + endif() - if(BLAS_LIBRARIES MATCHES "libmkl") - target_compile_definitions(phononcalc_lib PRIVATE MKL_LAPACKE - MULTITHREADED_BLAS) + if(BLAS_LIBRARIES MATCHES "libopenblas") + target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) + endif() endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) - endif() + target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) else() # Static link library add_library(phononcalc_lib STATIC ${SOURCES_PHONONCALC}) - if(OpenMP_FOUND) - target_link_libraries(phononcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phononcalc_lib PRIVATE OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib PRIVATE) + endif() + + target_compile_definitions(phononcalc_lib PRIVATE NO_INCLUDE_LAPACKE) else() - target_link_libraries(phononcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK) - endif() + if(OpenMP_FOUND) + target_link_libraries( + phononcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib PRIVATE BLAS::BLAS + LAPACK::LAPACK) + endif() - target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) + if(BLAS_LIBRARIES MATCHES "libmkl") + target_compile_definitions(phononcalc_lib PRIVATE MKL_BLAS + MULTITHREADED_BLAS) + endif() - if(BLAS_LIBRARIES MATCHES "libmkl") - target_compile_definitions(phononcalc_lib PRIVATE MKL_LAPACKE - MULTITHREADED_BLAS) + if(BLAS_LIBRARIES MATCHES "libopenblas") + target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) + endif() endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) - endif() + target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) endif() if(NOT BUILD_NANOBIND_MODULE) @@ -432,7 +483,17 @@ if(BUILD_NANOBIND_MODULE) target_link_libraries(_phononcalc PRIVATE phononcalc_lib) target_link_libraries(_recgrid PRIVATE recgrid_lib) - target_compile_definitions(_phono3py PRIVATE THM_EPSILON=1e-10) + if(BUILD_WITHOUT_LAPACKE) + target_compile_definitions(_phono3py PRIVATE NO_INCLUDE_LAPACKE + THM_EPSILON=1e-10) + else() + if(BLAS_LIBRARIES MATCHES "libmkl") + target_compile_definitions(_phono3py PRIVATE MKL_BLAS THM_EPSILON=1e-10) + else() + target_compile_definitions(_phono3py PRIVATE THM_EPSILON=1e-10) + endif() + endif() + install(TARGETS _phono3py LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) install(TARGETS _phononcalc LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) install(TARGETS _recgrid LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) diff --git a/c/_phono3py.c b/c/_phono3py.c deleted file mode 100644 index 27388cb5..00000000 --- a/c/_phono3py.c +++ /dev/null @@ -1,1877 +0,0 @@ -/* Copyright (C) 2015 Atsushi Togo */ -/* All rights reserved. */ - -/* This file is part of phonopy. */ - -/* Redistribution and use in source and binary forms, with or without */ -/* modification, are permitted provided that the following conditions */ -/* are met: */ - -/* * Redistributions of source code must retain the above copyright */ -/* notice, this list of conditions and the following disclaimer. */ - -/* * Redistributions in binary form must reproduce the above copyright */ -/* notice, this list of conditions and the following disclaimer in */ -/* the documentation and/or other materials provided with the */ -/* distribution. */ - -/* * Neither the name of the phonopy project nor the names of its */ -/* contributors may be used to endorse or promote products derived */ -/* from this software without specific prior written permission. */ - -/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */ -/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */ -/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */ -/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */ -/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */ -/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */ -/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */ -/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */ -/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */ -/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */ -/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ -/* POSSIBILITY OF SUCH DAMAGE. */ - -#include -#include -#include -#include -#include -#include -#include - -// #include "lapack_wrapper.h" -#include "phono3py.h" -#include "phonoc_array.h" - -static PyObject *py_get_interaction(PyObject *self, PyObject *args); -static PyObject *py_get_pp_collision(PyObject *self, PyObject *args); -static PyObject *py_get_pp_collision_with_sigma(PyObject *self, PyObject *args); -static PyObject *py_get_imag_self_energy_with_g(PyObject *self, PyObject *args); -static PyObject *py_get_detailed_imag_self_energy_with_g(PyObject *self, - PyObject *args); -static PyObject *py_get_real_self_energy_at_bands(PyObject *self, - PyObject *args); -static PyObject *py_get_real_self_energy_at_frequency_point(PyObject *self, - PyObject *args); -static PyObject *py_get_collision_matrix(PyObject *self, PyObject *args); -static PyObject *py_get_reducible_collision_matrix(PyObject *self, - PyObject *args); -static PyObject *py_symmetrize_collision_matrix(PyObject *self, PyObject *args); -static PyObject *py_expand_collision_matrix(PyObject *self, PyObject *args); -static PyObject *py_distribute_fc3(PyObject *self, PyObject *args); -static PyObject *py_rotate_delta_fc2s(PyObject *self, PyObject *args); -static PyObject *py_get_isotope_strength(PyObject *self, PyObject *args); -static PyObject *py_get_thm_isotope_strength(PyObject *self, PyObject *args); -static PyObject *py_get_permutation_symmetry_fc3(PyObject *self, - PyObject *args); -static PyObject *py_get_permutation_symmetry_compact_fc3(PyObject *self, - PyObject *args); -static PyObject *py_transpose_compact_fc3(PyObject *self, PyObject *args); -static PyObject *py_get_thm_relative_grid_address(PyObject *self, - PyObject *args); -static PyObject *py_get_neighboring_grid_points(PyObject *self, PyObject *args); -static PyObject *py_get_thm_integration_weights_at_grid_points(PyObject *self, - PyObject *args); -static PyObject *py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, - PyObject *args); -static PyObject *py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args); -static PyObject *py_get_triplets_integration_weights(PyObject *self, - PyObject *args); -static PyObject *py_get_triplets_integration_weights_with_sigma(PyObject *self, - PyObject *args); -static PyObject *py_get_grid_index_from_address(PyObject *self, PyObject *args); -static PyObject *py_get_gr_grid_addresses(PyObject *self, PyObject *args); -static PyObject *py_get_reciprocal_rotations(PyObject *self, PyObject *args); -static PyObject *py_transform_rotations(PyObject *self, PyObject *args); -static PyObject *py_get_snf3x3(PyObject *self, PyObject *args); -static PyObject *py_get_ir_grid_map(PyObject *self, PyObject *args); -static PyObject *py_get_bz_grid_addresses(PyObject *self, PyObject *args); -static PyObject *py_rotate_bz_grid_addresses(PyObject *self, PyObject *args); -static PyObject *py_diagonalize_collision_matrix(PyObject *self, - PyObject *args); -static PyObject *py_pinv_from_eigensolution(PyObject *self, PyObject *args); -static PyObject *py_get_default_colmat_solver(PyObject *self, PyObject *args); -static PyObject *py_lapacke_pinv(PyObject *self, PyObject *args); -static PyObject *py_get_omp_max_threads(PyObject *self, PyObject *args); - -static void show_colmat_info(const PyArrayObject *collision_matrix_py, - const long i_sigma, const long i_temp, - const long adrs_shift); -static Larray *convert_to_larray(PyArrayObject *npyary); -static Darray *convert_to_darray(PyArrayObject *npyary); - -struct module_state { - PyObject *error; -}; - -#if PY_MAJOR_VERSION >= 3 -#define GETSTATE(m) ((struct module_state *)PyModule_GetState(m)) -#else -#define GETSTATE(m) (&_state) -static struct module_state _state; -#endif - -static PyObject *error_out(PyObject *m) { - struct module_state *st = GETSTATE(m); - PyErr_SetString(st->error, "something bad happened"); - return NULL; -} - -static PyMethodDef _phono3py_methods[] = { - {"error_out", (PyCFunction)error_out, METH_NOARGS, NULL}, - {"interaction", (PyCFunction)py_get_interaction, METH_VARARGS, - "Interaction of triplets"}, - {"pp_collision", (PyCFunction)py_get_pp_collision, METH_VARARGS, - "Collision and ph-ph calculation"}, - {"pp_collision_with_sigma", (PyCFunction)py_get_pp_collision_with_sigma, - METH_VARARGS, "Collision and ph-ph calculation for smearing method"}, - {"imag_self_energy_with_g", (PyCFunction)py_get_imag_self_energy_with_g, - METH_VARARGS, "Imaginary part of self energy at frequency points with g"}, - {"detailed_imag_self_energy_with_g", - (PyCFunction)py_get_detailed_imag_self_energy_with_g, METH_VARARGS, - "Detailed contribution to imaginary part of self energy at frequency " - "points with g"}, - {"real_self_energy_at_bands", (PyCFunction)py_get_real_self_energy_at_bands, - METH_VARARGS, "Real part of self energy from third order force constants"}, - {"real_self_energy_at_frequency_point", - (PyCFunction)py_get_real_self_energy_at_frequency_point, METH_VARARGS, - "Real part of self energy from third order force constants at a frequency " - "point"}, - {"collision_matrix", (PyCFunction)py_get_collision_matrix, METH_VARARGS, - "Collision matrix with g"}, - {"reducible_collision_matrix", - (PyCFunction)py_get_reducible_collision_matrix, METH_VARARGS, - "Collision matrix with g for reducible grid points"}, - {"symmetrize_collision_matrix", (PyCFunction)py_symmetrize_collision_matrix, - METH_VARARGS, "Symmetrize collision matrix"}, - {"expand_collision_matrix", (PyCFunction)py_expand_collision_matrix, - METH_VARARGS, "Expand collision matrix"}, - {"distribute_fc3", (PyCFunction)py_distribute_fc3, METH_VARARGS, - "Distribute least fc3 to full fc3"}, - {"rotate_delta_fc2s", (PyCFunction)py_rotate_delta_fc2s, METH_VARARGS, - "Rotate delta fc2s"}, - {"isotope_strength", (PyCFunction)py_get_isotope_strength, METH_VARARGS, - "Isotope scattering strength"}, - {"thm_isotope_strength", (PyCFunction)py_get_thm_isotope_strength, - METH_VARARGS, "Isotope scattering strength for tetrahedron_method"}, - {"permutation_symmetry_fc3", (PyCFunction)py_get_permutation_symmetry_fc3, - METH_VARARGS, "Set permutation symmetry for fc3"}, - {"permutation_symmetry_compact_fc3", - (PyCFunction)py_get_permutation_symmetry_compact_fc3, METH_VARARGS, - "Set permutation symmetry for compact-fc3"}, - {"transpose_compact_fc3", (PyCFunction)py_transpose_compact_fc3, - METH_VARARGS, "Transpose compact fc3"}, - {"tetrahedra_relative_grid_address", - (PyCFunction)py_get_thm_relative_grid_address, METH_VARARGS, - "Relative grid addresses of vertices of 24 tetrahedra"}, - {"neighboring_grid_points", (PyCFunction)py_get_neighboring_grid_points, - METH_VARARGS, "Neighboring grid points by relative grid addresses"}, - {"integration_weights_at_grid_points", - (PyCFunction)py_get_thm_integration_weights_at_grid_points, METH_VARARGS, - "Integration weights of tetrahedron method at grid points"}, - {"triplets_reciprocal_mesh_at_q", - (PyCFunction)py_tpl_get_triplets_reciprocal_mesh_at_q, METH_VARARGS, - "Triplets on reciprocal mesh points at a specific q-point"}, - {"BZ_triplets_at_q", (PyCFunction)py_tpl_get_BZ_triplets_at_q, METH_VARARGS, - "Triplets in reciprocal primitive lattice are transformed to those in " - "BZ."}, - {"triplets_integration_weights", - (PyCFunction)py_get_triplets_integration_weights, METH_VARARGS, - "Integration weights of tetrahedron method for triplets"}, - {"triplets_integration_weights_with_sigma", - (PyCFunction)py_get_triplets_integration_weights_with_sigma, METH_VARARGS, - "Integration weights of smearing method for triplets"}, - {"grid_index_from_address", (PyCFunction)py_get_grid_index_from_address, - METH_VARARGS, "Grid index from grid address"}, - {"ir_grid_map", (PyCFunction)py_get_ir_grid_map, METH_VARARGS, - "Reciprocal mesh points with ir grid mapping table"}, - {"gr_grid_addresses", (PyCFunction)py_get_gr_grid_addresses, METH_VARARGS, - "Get generalized regular grid addresses"}, - {"reciprocal_rotations", (PyCFunction)py_get_reciprocal_rotations, - METH_VARARGS, "Return rotation matrices in reciprocal space"}, - {"transform_rotations", (PyCFunction)py_transform_rotations, METH_VARARGS, - "Transform rotations to those in generalized regular grid"}, - {"snf3x3", (PyCFunction)py_get_snf3x3, METH_VARARGS, - "Get Smith formal form for 3x3 integer matrix"}, - {"bz_grid_addresses", (PyCFunction)py_get_bz_grid_addresses, METH_VARARGS, - "Get grid addresses including Brillouin zone surface"}, - {"rotate_bz_grid_index", (PyCFunction)py_rotate_bz_grid_addresses, - METH_VARARGS, "Rotate grid point considering Brillouin zone surface"}, - {"diagonalize_collision_matrix", - (PyCFunction)py_diagonalize_collision_matrix, METH_VARARGS, - "Diagonalize and optionally pseudo-inverse using Lapack dsyev(d)"}, - {"pinv_from_eigensolution", (PyCFunction)py_pinv_from_eigensolution, - METH_VARARGS, "Pseudo-inverse from eigensolution"}, - {"default_colmat_solver", (PyCFunction)py_get_default_colmat_solver, - METH_VARARGS, "Return default collison matrix solver by integer value"}, - {"lapacke_pinv", (PyCFunction)py_lapacke_pinv, METH_VARARGS, - "Pseudo inversion using lapacke."}, - {"omp_max_threads", py_get_omp_max_threads, METH_VARARGS, - "Return openmp max number of threads. Return 0 unless openmp is " - "activated. "}, - {NULL, NULL, 0, NULL}}; - -#if PY_MAJOR_VERSION >= 3 - -static int _phono3py_traverse(PyObject *m, visitproc visit, void *arg) { - Py_VISIT(GETSTATE(m)->error); - return 0; -} - -static int _phono3py_clear(PyObject *m) { - Py_CLEAR(GETSTATE(m)->error); - return 0; -} - -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, "_phono3py", NULL, - sizeof(struct module_state), _phono3py_methods, NULL, - _phono3py_traverse, _phono3py_clear, NULL}; - -#define INITERROR return NULL - -PyObject *PyInit__phono3py(void) -#else -#define INITERROR return - -void init_phono3py(void) -#endif -{ -#if PY_MAJOR_VERSION >= 3 - PyObject *module = PyModule_Create(&moduledef); -#else - PyObject *module = Py_InitModule("_phono3py", _phono3py_methods); -#endif - struct module_state *st; - - if (module == NULL) INITERROR; - st = GETSTATE(module); - - st->error = PyErr_NewException("_phono3py.Error", NULL, NULL); - if (st->error == NULL) { - Py_DECREF(module); - INITERROR; - } - -#if PY_MAJOR_VERSION >= 3 - return module; -#endif -} - -static PyObject *py_get_interaction(PyObject *self, PyObject *args) { - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_g_zero; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_triplets; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_svecs; - PyArrayObject *py_multi; - PyArrayObject *py_fc3; - PyArrayObject *py_masses; - PyArrayObject *py_p2s_map; - PyArrayObject *py_s2p_map; - PyArrayObject *py_band_indices; - PyArrayObject *py_all_shortest; - double cutoff_frequency; - long symmetrize_fc3_q; - long make_r0_average; - long openmp_per_triplets; - - Darray *fc3_normal_squared; - Darray *freqs; - _lapack_complex_double *eigvecs; - long(*triplets)[3]; - long num_triplets; - char *g_zero; - long(*bz_grid_addresses)[3]; - long *D_diag; - long(*Q)[3]; - double *fc3; - double(*svecs)[3]; - long(*multi)[2]; - double *masses; - char *all_shortest; - long *p2s; - long *s2p; - long *band_indices; - long multi_dims[2]; - long i; - long is_compact_fc3; - - if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOllOdl", &py_fc3_normal_squared, - &py_g_zero, &py_frequencies, &py_eigenvectors, - &py_triplets, &py_bz_grid_addresses, &py_D_diag, - &py_Q, &py_fc3, &py_svecs, &py_multi, &py_masses, - &py_p2s_map, &py_s2p_map, &py_band_indices, - &symmetrize_fc3_q, &make_r0_average, &py_all_shortest, - &cutoff_frequency, &openmp_per_triplets)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - freqs = convert_to_darray(py_frequencies); - /* npy_cdouble and lapack_complex_double may not be compatible. */ - /* So eigenvectors should not be used in Python side */ - eigvecs = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - g_zero = (char *)PyArray_DATA(py_g_zero); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - fc3 = (double *)PyArray_DATA(py_fc3); - if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) { - is_compact_fc3 = 0; - } else { - is_compact_fc3 = 1; - } - svecs = (double(*)[3])PyArray_DATA(py_svecs); - for (i = 0; i < 2; i++) { - multi_dims[i] = PyArray_DIMS(py_multi)[i]; - } - multi = (long(*)[2])PyArray_DATA(py_multi); - masses = (double *)PyArray_DATA(py_masses); - p2s = (long *)PyArray_DATA(py_p2s_map); - s2p = (long *)PyArray_DATA(py_s2p_map); - band_indices = (long *)PyArray_DATA(py_band_indices); - all_shortest = (char *)PyArray_DATA(py_all_shortest); - - ph3py_get_interaction(fc3_normal_squared, g_zero, freqs, eigvecs, triplets, - num_triplets, bz_grid_addresses, D_diag, Q, fc3, - is_compact_fc3, svecs, multi_dims, multi, masses, p2s, - s2p, band_indices, symmetrize_fc3_q, make_r0_average, - all_shortest, cutoff_frequency, openmp_per_triplets); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - free(freqs); - freqs = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_pp_collision(PyObject *self, PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_bz_map; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_fc3; - PyArrayObject *py_svecs; - PyArrayObject *py_multi; - PyArrayObject *py_masses; - PyArrayObject *py_p2s_map; - PyArrayObject *py_s2p_map; - PyArrayObject *py_band_indices; - PyArrayObject *py_temperatures; - PyArrayObject *py_all_shortest; - double cutoff_frequency; - long is_NU; - long symmetrize_fc3_q; - long make_r0_average; - long bz_grid_type; - long openmp_per_triplets; - - double *gamma; - long(*relative_grid_address)[4][3]; - double *frequencies; - _lapack_complex_double *eigenvectors; - long(*triplets)[3]; - long num_triplets; - long *triplet_weights; - long(*bz_grid_addresses)[3]; - long *bz_map; - long *D_diag; - long(*Q)[3]; - double *fc3; - double(*svecs)[3]; - long(*multi)[2]; - double *masses; - long *p2s; - long *s2p; - Larray *band_indices; - Darray *temperatures; - char *all_shortest; - long multi_dims[2]; - long i; - long is_compact_fc3; - - if (!PyArg_ParseTuple( - args, "OOOOOOOOlOOOOOOOOOOlllOdl", &py_gamma, - &py_relative_grid_address, &py_frequencies, &py_eigenvectors, - &py_triplets, &py_triplet_weights, &py_bz_grid_addresses, - &py_bz_map, &bz_grid_type, &py_D_diag, &py_Q, &py_fc3, &py_svecs, - &py_multi, &py_masses, &py_p2s_map, &py_s2p_map, &py_band_indices, - &py_temperatures, &is_NU, &symmetrize_fc3_q, &make_r0_average, - &py_all_shortest, &cutoff_frequency, &openmp_per_triplets)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - frequencies = (double *)PyArray_DATA(py_frequencies); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - bz_map = (long *)PyArray_DATA(py_bz_map); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - fc3 = (double *)PyArray_DATA(py_fc3); - if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) { - is_compact_fc3 = 0; - } else { - is_compact_fc3 = 1; - } - svecs = (double(*)[3])PyArray_DATA(py_svecs); - for (i = 0; i < 2; i++) { - multi_dims[i] = PyArray_DIMS(py_multi)[i]; - } - multi = (long(*)[2])PyArray_DATA(py_multi); - masses = (double *)PyArray_DATA(py_masses); - p2s = (long *)PyArray_DATA(py_p2s_map); - s2p = (long *)PyArray_DATA(py_s2p_map); - band_indices = convert_to_larray(py_band_indices); - temperatures = convert_to_darray(py_temperatures); - all_shortest = (char *)PyArray_DATA(py_all_shortest); - - ph3py_get_pp_collision( - gamma, relative_grid_address, frequencies, eigenvectors, triplets, - num_triplets, triplet_weights, bz_grid_addresses, bz_map, bz_grid_type, - D_diag, Q, fc3, is_compact_fc3, svecs, multi_dims, multi, masses, p2s, - s2p, band_indices, temperatures, is_NU, symmetrize_fc3_q, - make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets); - - free(band_indices); - band_indices = NULL; - free(temperatures); - temperatures = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_pp_collision_with_sigma(PyObject *self, - PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_fc3; - PyArrayObject *py_svecs; - PyArrayObject *py_multi; - PyArrayObject *py_masses; - PyArrayObject *py_p2s_map; - PyArrayObject *py_s2p_map; - PyArrayObject *py_band_indices; - PyArrayObject *py_temperatures; - PyArrayObject *py_all_shortest; - long is_NU; - long symmetrize_fc3_q; - double sigma; - double sigma_cutoff; - long make_r0_average; - double cutoff_frequency; - long openmp_per_triplets; - - double *gamma; - double *frequencies; - _lapack_complex_double *eigenvectors; - long(*triplets)[3]; - long num_triplets; - long *triplet_weights; - long(*bz_grid_addresses)[3]; - long *D_diag; - long(*Q)[3]; - double *fc3; - double(*svecs)[3]; - long(*multi)[2]; - double *masses; - long *p2s; - long *s2p; - Larray *band_indices; - Darray *temperatures; - char *all_shortest; - long multi_dims[2]; - long i; - long is_compact_fc3; - - if (!PyArg_ParseTuple( - args, "OddOOOOOOOOOOOOOOOlllOdl", &py_gamma, &sigma, &sigma_cutoff, - &py_frequencies, &py_eigenvectors, &py_triplets, - &py_triplet_weights, &py_bz_grid_addresses, &py_D_diag, &py_Q, - &py_fc3, &py_svecs, &py_multi, &py_masses, &py_p2s_map, &py_s2p_map, - &py_band_indices, &py_temperatures, &is_NU, &symmetrize_fc3_q, - &make_r0_average, &py_all_shortest, &cutoff_frequency, - &openmp_per_triplets)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - frequencies = (double *)PyArray_DATA(py_frequencies); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - fc3 = (double *)PyArray_DATA(py_fc3); - if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) { - is_compact_fc3 = 0; - } else { - is_compact_fc3 = 1; - } - svecs = (double(*)[3])PyArray_DATA(py_svecs); - for (i = 0; i < 2; i++) { - multi_dims[i] = PyArray_DIMS(py_multi)[i]; - } - multi = (long(*)[2])PyArray_DATA(py_multi); - masses = (double *)PyArray_DATA(py_masses); - p2s = (long *)PyArray_DATA(py_p2s_map); - s2p = (long *)PyArray_DATA(py_s2p_map); - band_indices = convert_to_larray(py_band_indices); - temperatures = convert_to_darray(py_temperatures); - all_shortest = (char *)PyArray_DATA(py_all_shortest); - - ph3py_get_pp_collision_with_sigma( - gamma, sigma, sigma_cutoff, frequencies, eigenvectors, triplets, - num_triplets, triplet_weights, bz_grid_addresses, D_diag, Q, fc3, - is_compact_fc3, svecs, multi_dims, multi, masses, p2s, s2p, - band_indices, temperatures, is_NU, symmetrize_fc3_q, make_r0_average, - all_shortest, cutoff_frequency, openmp_per_triplets); - - free(band_indices); - band_indices = NULL; - free(temperatures); - temperatures = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_imag_self_energy_with_g(PyObject *self, - PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_g; - PyArrayObject *py_g_zero; - double cutoff_frequency, temperature; - long frequency_point_index; - - Darray *fc3_normal_squared; - double *gamma; - double *g; - char *g_zero; - double *frequencies; - long(*triplets)[3]; - long *triplet_weights; - long num_frequency_points; - - if (!PyArg_ParseTuple(args, "OOOOOdOOdl", &py_gamma, &py_fc3_normal_squared, - &py_triplets, &py_triplet_weights, &py_frequencies, - &temperature, &py_g, &py_g_zero, &cutoff_frequency, - &frequency_point_index)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - gamma = (double *)PyArray_DATA(py_gamma); - g = (double *)PyArray_DATA(py_g); - g_zero = (char *)PyArray_DATA(py_g_zero); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - num_frequency_points = (long)PyArray_DIMS(py_g)[2]; - - ph3py_get_imag_self_energy_at_bands_with_g( - gamma, fc3_normal_squared, frequencies, triplets, triplet_weights, g, - g_zero, temperature, cutoff_frequency, num_frequency_points, - frequency_point_index); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_detailed_imag_self_energy_with_g(PyObject *self, - PyObject *args) { - PyArrayObject *py_gamma_detail; - PyArrayObject *py_gamma_N; - PyArrayObject *py_gamma_U; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_g; - PyArrayObject *py_g_zero; - double cutoff_frequency, temperature; - - Darray *fc3_normal_squared; - double *gamma_detail; - double *gamma_N; - double *gamma_U; - double *g; - char *g_zero; - double *frequencies; - long(*triplets)[3]; - long *triplet_weights; - long(*bz_grid_addresses)[3]; - - if (!PyArg_ParseTuple(args, "OOOOOOOOdOOd", &py_gamma_detail, &py_gamma_N, - &py_gamma_U, &py_fc3_normal_squared, &py_triplets, - &py_triplet_weights, &py_bz_grid_addresses, - &py_frequencies, &temperature, &py_g, &py_g_zero, - &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - gamma_detail = (double *)PyArray_DATA(py_gamma_detail); - gamma_N = (double *)PyArray_DATA(py_gamma_N); - gamma_U = (double *)PyArray_DATA(py_gamma_U); - g = (double *)PyArray_DATA(py_g); - g_zero = (char *)PyArray_DATA(py_g_zero); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - - ph3py_get_detailed_imag_self_energy_at_bands_with_g( - gamma_detail, gamma_N, gamma_U, fc3_normal_squared, frequencies, - triplets, triplet_weights, bz_grid_addresses, g, g_zero, temperature, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_real_self_energy_at_bands(PyObject *self, - PyObject *args) { - PyArrayObject *py_shift; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_band_indices; - double epsilon, unit_conversion_factor, cutoff_frequency, temperature; - - Darray *fc3_normal_squared; - double *shift; - double *frequencies; - long *band_indices; - long(*triplets)[3]; - long *triplet_weights; - - if (!PyArg_ParseTuple(args, "OOOOOOdddd", &py_shift, &py_fc3_normal_squared, - &py_triplets, &py_triplet_weights, &py_frequencies, - &py_band_indices, &temperature, &epsilon, - &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - shift = (double *)PyArray_DATA(py_shift); - frequencies = (double *)PyArray_DATA(py_frequencies); - band_indices = (long *)PyArray_DATA(py_band_indices); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - - ph3py_get_real_self_energy_at_bands( - shift, fc3_normal_squared, band_indices, frequencies, triplets, - triplet_weights, epsilon, temperature, unit_conversion_factor, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_real_self_energy_at_frequency_point(PyObject *self, - PyObject *args) { - PyArrayObject *py_shift; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_band_indices; - double frequency_point, epsilon, unit_conversion_factor, cutoff_frequency; - double temperature; - - Darray *fc3_normal_squared; - double *shift; - double *frequencies; - long *band_indices; - long(*triplets)[3]; - long *triplet_weights; - - if (!PyArg_ParseTuple(args, "OdOOOOOdddd", &py_shift, &frequency_point, - &py_fc3_normal_squared, &py_triplets, - &py_triplet_weights, &py_frequencies, - &py_band_indices, &temperature, &epsilon, - &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - shift = (double *)PyArray_DATA(py_shift); - frequencies = (double *)PyArray_DATA(py_frequencies); - band_indices = (long *)PyArray_DATA(py_band_indices); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - - ph3py_get_real_self_energy_at_frequency_point( - shift, frequency_point, fc3_normal_squared, band_indices, frequencies, - triplets, triplet_weights, epsilon, temperature, unit_conversion_factor, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_collision_matrix(PyObject *self, PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplets_map; - PyArrayObject *py_map_q; - PyArrayObject *py_g; - PyArrayObject *py_rotated_grid_points; - PyArrayObject *py_rotations_cartesian; - double temperature, unit_conversion_factor, cutoff_frequency; - - Darray *fc3_normal_squared; - double *collision_matrix; - double *g; - double *frequencies; - long(*triplets)[3]; - long *triplets_map; - long *map_q; - long *rotated_grid_points; - long num_gp, num_ir_gp, num_rot; - double *rotations_cartesian; - - if (!PyArg_ParseTuple( - args, "OOOOOOOOOddd", &py_collision_matrix, &py_fc3_normal_squared, - &py_frequencies, &py_g, &py_triplets, &py_triplets_map, &py_map_q, - &py_rotated_grid_points, &py_rotations_cartesian, &temperature, - &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - g = (double *)PyArray_DATA(py_g); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplets_map = (long *)PyArray_DATA(py_triplets_map); - num_gp = (long)PyArray_DIMS(py_triplets_map)[0]; - map_q = (long *)PyArray_DATA(py_map_q); - rotated_grid_points = (long *)PyArray_DATA(py_rotated_grid_points); - num_ir_gp = (long)PyArray_DIMS(py_rotated_grid_points)[0]; - num_rot = (long)PyArray_DIMS(py_rotated_grid_points)[1]; - rotations_cartesian = (double *)PyArray_DATA(py_rotations_cartesian); - - assert(num_rot == PyArray_DIMS(py_rotations_cartesian)[0]); - assert(num_gp == PyArray_DIMS(py_frequencies)[0]); - - ph3py_get_collision_matrix(collision_matrix, fc3_normal_squared, - frequencies, triplets, triplets_map, map_q, - rotated_grid_points, rotations_cartesian, g, - num_ir_gp, num_gp, num_rot, temperature, - unit_conversion_factor, cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_reducible_collision_matrix(PyObject *self, - PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplets_map; - PyArrayObject *py_map_q; - PyArrayObject *py_g; - double temperature, unit_conversion_factor, cutoff_frequency; - - Darray *fc3_normal_squared; - double *collision_matrix; - double *g; - double *frequencies; - long(*triplets)[3]; - long *triplets_map; - long num_gp; - long *map_q; - - if (!PyArg_ParseTuple( - args, "OOOOOOOddd", &py_collision_matrix, &py_fc3_normal_squared, - &py_frequencies, &py_g, &py_triplets, &py_triplets_map, &py_map_q, - &temperature, &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - g = (double *)PyArray_DATA(py_g); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplets_map = (long *)PyArray_DATA(py_triplets_map); - num_gp = (long)PyArray_DIMS(py_triplets_map)[0]; - map_q = (long *)PyArray_DATA(py_map_q); - - ph3py_get_reducible_collision_matrix( - collision_matrix, fc3_normal_squared, frequencies, triplets, - triplets_map, map_q, g, num_gp, temperature, unit_conversion_factor, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_symmetrize_collision_matrix(PyObject *self, - PyObject *args) { - PyArrayObject *py_collision_matrix; - - double *collision_matrix; - long num_band, num_grid_points, num_temp, num_sigma; - long num_column; - - if (!PyArg_ParseTuple(args, "O", &py_collision_matrix)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - num_sigma = (long)PyArray_DIMS(py_collision_matrix)[0]; - num_temp = (long)PyArray_DIMS(py_collision_matrix)[1]; - num_grid_points = (long)PyArray_DIMS(py_collision_matrix)[2]; - num_band = (long)PyArray_DIMS(py_collision_matrix)[3]; - - if (PyArray_NDIM(py_collision_matrix) == 8) { - num_column = num_grid_points * num_band * 3; - } else { - num_column = num_grid_points * num_band; - } - - ph3py_symmetrize_collision_matrix(collision_matrix, num_column, num_temp, - num_sigma); - - Py_RETURN_NONE; -} - -static PyObject *py_expand_collision_matrix(PyObject *self, PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_ir_grid_points; - PyArrayObject *py_rot_grid_points; - - double *collision_matrix; - long *rot_grid_points; - long *ir_grid_points; - long num_band, num_grid_points, num_temp, num_sigma, num_rot, num_ir_gp; - - if (!PyArg_ParseTuple(args, "OOO", &py_collision_matrix, &py_ir_grid_points, - &py_rot_grid_points)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - rot_grid_points = (long *)PyArray_DATA(py_rot_grid_points); - ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points); - num_sigma = (long)PyArray_DIMS(py_collision_matrix)[0]; - num_temp = (long)PyArray_DIMS(py_collision_matrix)[1]; - num_grid_points = (long)PyArray_DIMS(py_collision_matrix)[2]; - num_band = (long)PyArray_DIMS(py_collision_matrix)[3]; - num_rot = (long)PyArray_DIMS(py_rot_grid_points)[0]; - num_ir_gp = (long)PyArray_DIMS(py_ir_grid_points)[0]; - - ph3py_expand_collision_matrix(collision_matrix, rot_grid_points, - ir_grid_points, num_ir_gp, num_grid_points, - num_rot, num_sigma, num_temp, num_band); - - Py_RETURN_NONE; -} - -static PyObject *py_get_isotope_strength(PyObject *self, PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_band_indices; - PyArrayObject *py_mass_variances; - PyArrayObject *py_ir_grid_points; - PyArrayObject *py_weights; - - long grid_point; - double cutoff_frequency; - double sigma; - - double *gamma; - double *frequencies; - long *ir_grid_points; - double *weights; - _lapack_complex_double *eigenvectors; - long *band_indices; - double *mass_variances; - long num_band, num_band0, num_ir_grid_points; - - if (!PyArg_ParseTuple(args, "OlOOOOOOdd", &py_gamma, &grid_point, - &py_ir_grid_points, &py_weights, &py_mass_variances, - &py_frequencies, &py_eigenvectors, &py_band_indices, - &sigma, &cutoff_frequency)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - frequencies = (double *)PyArray_DATA(py_frequencies); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points); - weights = (double *)PyArray_DATA(py_weights); - band_indices = (long *)PyArray_DATA(py_band_indices); - mass_variances = (double *)PyArray_DATA(py_mass_variances); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - num_band0 = (long)PyArray_DIMS(py_band_indices)[0]; - num_ir_grid_points = (long)PyArray_DIMS(py_ir_grid_points)[0]; - - ph3py_get_isotope_scattering_strength( - gamma, grid_point, ir_grid_points, weights, mass_variances, frequencies, - eigenvectors, num_ir_grid_points, band_indices, num_band, num_band0, - sigma, cutoff_frequency); - - Py_RETURN_NONE; -} - -static PyObject *py_get_thm_isotope_strength(PyObject *self, PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_band_indices; - PyArrayObject *py_mass_variances; - PyArrayObject *py_ir_grid_points; - PyArrayObject *py_weights; - PyArrayObject *py_integration_weights; - long grid_point; - double cutoff_frequency; - - double *gamma; - double *frequencies; - long *ir_grid_points; - double *weights; - _lapack_complex_double *eigenvectors; - long *band_indices; - double *mass_variances; - long num_band, num_band0, num_ir_grid_points; - double *integration_weights; - - if (!PyArg_ParseTuple(args, "OlOOOOOOOd", &py_gamma, &grid_point, - &py_ir_grid_points, &py_weights, &py_mass_variances, - &py_frequencies, &py_eigenvectors, &py_band_indices, - &py_integration_weights, &cutoff_frequency)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - frequencies = (double *)PyArray_DATA(py_frequencies); - ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points); - weights = (double *)PyArray_DATA(py_weights); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - band_indices = (long *)PyArray_DATA(py_band_indices); - mass_variances = (double *)PyArray_DATA(py_mass_variances); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - num_band0 = (long)PyArray_DIMS(py_band_indices)[0]; - integration_weights = (double *)PyArray_DATA(py_integration_weights); - num_ir_grid_points = (long)PyArray_DIMS(py_ir_grid_points)[0]; - - ph3py_get_thm_isotope_scattering_strength( - gamma, grid_point, ir_grid_points, weights, mass_variances, frequencies, - eigenvectors, num_ir_grid_points, band_indices, num_band, num_band0, - integration_weights, cutoff_frequency); - - Py_RETURN_NONE; -} - -static PyObject *py_distribute_fc3(PyObject *self, PyObject *args) { - PyArrayObject *force_constants_third; - long target; - long source; - PyArrayObject *rotation_cart_inv; - PyArrayObject *atom_mapping_py; - - double *fc3; - double *rot_cart_inv; - long *atom_mapping; - long num_atom; - - if (!PyArg_ParseTuple(args, "OllOO", &force_constants_third, &target, - &source, &atom_mapping_py, &rotation_cart_inv)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(force_constants_third); - rot_cart_inv = (double *)PyArray_DATA(rotation_cart_inv); - atom_mapping = (long *)PyArray_DATA(atom_mapping_py); - num_atom = (long)PyArray_DIMS(atom_mapping_py)[0]; - - ph3py_distribute_fc3(fc3, target, source, atom_mapping, num_atom, - rot_cart_inv); - - Py_RETURN_NONE; -} - -static PyObject *py_rotate_delta_fc2s(PyObject *self, PyObject *args) { - PyArrayObject *py_fc3; - PyArrayObject *py_delta_fc2s; - PyArrayObject *py_inv_U; - PyArrayObject *py_site_sym_cart; - PyArrayObject *py_rot_map_syms; - - double(*fc3)[3][3][3]; - double(*delta_fc2s)[3][3]; - double *inv_U; - double(*site_sym_cart)[3][3]; - long *rot_map_syms; - long num_atom, num_disp, num_site_sym; - - if (!PyArg_ParseTuple(args, "OOOOO", &py_fc3, &py_delta_fc2s, &py_inv_U, - &py_site_sym_cart, &py_rot_map_syms)) { - return NULL; - } - - /* (num_atom, num_atom, 3, 3, 3) */ - fc3 = (double(*)[3][3][3])PyArray_DATA(py_fc3); - /* (n_u1, num_atom, num_atom, 3, 3) */ - delta_fc2s = (double(*)[3][3])PyArray_DATA(py_delta_fc2s); - /* (3, n_u1 * n_sym) */ - inv_U = (double *)PyArray_DATA(py_inv_U); - /* (n_sym, 3, 3) */ - site_sym_cart = (double(*)[3][3])PyArray_DATA(py_site_sym_cart); - /* (n_sym, natom) */ - rot_map_syms = (long *)PyArray_DATA(py_rot_map_syms); - - num_atom = (long)PyArray_DIMS(py_fc3)[0]; - num_disp = (long)PyArray_DIMS(py_delta_fc2s)[0]; - num_site_sym = (long)PyArray_DIMS(py_site_sym_cart)[0]; - - ph3py_rotate_delta_fc2(fc3, delta_fc2s, inv_U, site_sym_cart, rot_map_syms, - num_atom, num_site_sym, num_disp); - - Py_RETURN_NONE; -} - -static PyObject *py_get_permutation_symmetry_fc3(PyObject *self, - PyObject *args) { - PyArrayObject *py_fc3; - - double *fc3; - long num_atom; - - if (!PyArg_ParseTuple(args, "O", &py_fc3)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(py_fc3); - num_atom = (long)PyArray_DIMS(py_fc3)[0]; - - ph3py_get_permutation_symmetry_fc3(fc3, num_atom); - - Py_RETURN_NONE; -} - -static PyObject *py_get_permutation_symmetry_compact_fc3(PyObject *self, - PyObject *args) { - PyArrayObject *py_fc3; - PyArrayObject *py_permutations; - PyArrayObject *py_s2pp_map; - PyArrayObject *py_p2s_map; - PyArrayObject *py_nsym_list; - - double *fc3; - long *s2pp; - long *p2s; - long *nsym_list; - long *perms; - long n_patom, n_satom; - - if (!PyArg_ParseTuple(args, "OOOOO", &py_fc3, &py_permutations, - &py_s2pp_map, &py_p2s_map, &py_nsym_list)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(py_fc3); - perms = (long *)PyArray_DATA(py_permutations); - s2pp = (long *)PyArray_DATA(py_s2pp_map); - p2s = (long *)PyArray_DATA(py_p2s_map); - nsym_list = (long *)PyArray_DATA(py_nsym_list); - n_patom = (long)PyArray_DIMS(py_fc3)[0]; - n_satom = (long)PyArray_DIMS(py_fc3)[1]; - - ph3py_get_permutation_symmetry_compact_fc3(fc3, p2s, s2pp, nsym_list, perms, - n_satom, n_patom); - - Py_RETURN_NONE; -} - -static PyObject *py_transpose_compact_fc3(PyObject *self, PyObject *args) { - PyArrayObject *py_fc3; - PyArrayObject *py_permutations; - PyArrayObject *py_s2pp_map; - PyArrayObject *py_p2s_map; - PyArrayObject *py_nsym_list; - long t_type; - - double *fc3; - long *s2pp; - long *p2s; - long *nsym_list; - long *perms; - long n_patom, n_satom; - - if (!PyArg_ParseTuple(args, "OOOOOl", &py_fc3, &py_permutations, - &py_s2pp_map, &py_p2s_map, &py_nsym_list, &t_type)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(py_fc3); - perms = (long *)PyArray_DATA(py_permutations); - s2pp = (long *)PyArray_DATA(py_s2pp_map); - p2s = (long *)PyArray_DATA(py_p2s_map); - nsym_list = (long *)PyArray_DATA(py_nsym_list); - n_patom = (long)PyArray_DIMS(py_fc3)[0]; - n_satom = (long)PyArray_DIMS(py_fc3)[1]; - - ph3py_transpose_compact_fc3(fc3, p2s, s2pp, nsym_list, perms, n_satom, - n_patom, t_type); - - Py_RETURN_NONE; -} - -static PyObject *py_get_thm_relative_grid_address(PyObject *self, - PyObject *args) { - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_reciprocal_lattice_py; - - long(*relative_grid_address)[4][3]; - double(*reciprocal_lattice)[3]; - - if (!PyArg_ParseTuple(args, "OO", &py_relative_grid_address, - &py_reciprocal_lattice_py)) { - return NULL; - } - - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice_py); - - ph3py_get_relative_grid_address(relative_grid_address, reciprocal_lattice); - - Py_RETURN_NONE; -} - -static PyObject *py_get_neighboring_grid_points(PyObject *self, - PyObject *args) { - PyArrayObject *py_relative_grid_points; - PyArrayObject *py_grid_points; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_D_diag; - PyArrayObject *py_bz_grid_address; - PyArrayObject *py_bz_map; - long bz_grid_type; - - long *relative_grid_points; - long *grid_points; - long num_grid_points, num_relative_grid_address; - long(*relative_grid_address)[3]; - long *D_diag; - long(*bz_grid_address)[3]; - long *bz_map; - - if (!PyArg_ParseTuple(args, "OOOOOOl", &py_relative_grid_points, - &py_grid_points, &py_relative_grid_address, - &py_D_diag, &py_bz_grid_address, &py_bz_map, - &bz_grid_type)) { - return NULL; - } - - relative_grid_points = (long *)PyArray_DATA(py_relative_grid_points); - grid_points = (long *)PyArray_DATA(py_grid_points); - num_grid_points = (long)PyArray_DIMS(py_grid_points)[0]; - relative_grid_address = (long(*)[3])PyArray_DATA(py_relative_grid_address); - num_relative_grid_address = (long)PyArray_DIMS(py_relative_grid_address)[0]; - D_diag = (long *)PyArray_DATA(py_D_diag); - bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address); - bz_map = (long *)PyArray_DATA(py_bz_map); - - ph3py_get_neighboring_gird_points( - relative_grid_points, grid_points, relative_grid_address, D_diag, - bz_grid_address, bz_map, bz_grid_type, num_grid_points, - num_relative_grid_address); - - Py_RETURN_NONE; -} - -static PyObject *py_get_thm_integration_weights_at_grid_points(PyObject *self, - PyObject *args) { - PyArrayObject *py_iw; - PyArrayObject *py_frequency_points; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_D_diag; - PyArrayObject *py_grid_points; - PyArrayObject *py_frequencies; - PyArrayObject *py_bz_grid_address; - PyArrayObject *py_gp2irgp_map; - PyArrayObject *py_bz_map; - long bz_grid_type; - char *function; - - double *iw; - double *frequency_points; - long num_frequency_points, num_band, num_gp; - long(*relative_grid_address)[4][3]; - long *D_diag; - long *grid_points; - long(*bz_grid_address)[3]; - long *bz_map; - long *gp2irgp_map; - double *frequencies; - - if (!PyArg_ParseTuple(args, "OOOOOOOOOls", &py_iw, &py_frequency_points, - &py_relative_grid_address, &py_D_diag, - &py_grid_points, &py_frequencies, &py_bz_grid_address, - &py_bz_map, &py_gp2irgp_map, &bz_grid_type, - &function)) { - return NULL; - } - - iw = (double *)PyArray_DATA(py_iw); - frequency_points = (double *)PyArray_DATA(py_frequency_points); - num_frequency_points = (long)PyArray_DIMS(py_frequency_points)[0]; - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - D_diag = (long *)PyArray_DATA(py_D_diag); - grid_points = (long *)PyArray_DATA(py_grid_points); - num_gp = (long)PyArray_DIMS(py_grid_points)[0]; - bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address); - bz_map = (long *)PyArray_DATA(py_bz_map); - gp2irgp_map = (long *)PyArray_DATA(py_gp2irgp_map); - frequencies = (double *)PyArray_DATA(py_frequencies); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - - ph3py_get_thm_integration_weights_at_grid_points( - iw, frequency_points, num_frequency_points, num_band, num_gp, - relative_grid_address, D_diag, grid_points, bz_grid_address, bz_map, - bz_grid_type, frequencies, gp2irgp_map, function[0]); - - Py_RETURN_NONE; -} - -static PyObject *py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, - PyObject *args) { - PyArrayObject *py_map_triplets; - PyArrayObject *py_map_q; - PyArrayObject *py_D_diag; - PyArrayObject *py_rotations; - long fixed_grid_number; - long is_time_reversal; - long swappable; - - long *map_triplets; - long *map_q; - long *D_diag; - long(*rot)[3][3]; - long num_rot; - long num_ir; - - if (!PyArg_ParseTuple(args, "OOlOlOl", &py_map_triplets, &py_map_q, - &fixed_grid_number, &py_D_diag, &is_time_reversal, - &py_rotations, &swappable)) { - return NULL; - } - - map_triplets = (long *)PyArray_DATA(py_map_triplets); - map_q = (long *)PyArray_DATA(py_map_q); - D_diag = (long *)PyArray_DATA(py_D_diag); - rot = (long(*)[3][3])PyArray_DATA(py_rotations); - num_rot = (long)PyArray_DIMS(py_rotations)[0]; - num_ir = ph3py_get_triplets_reciprocal_mesh_at_q( - map_triplets, map_q, fixed_grid_number, D_diag, is_time_reversal, - num_rot, rot, swappable); - - return PyLong_FromLong(num_ir); -} - -static PyObject *py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args) { - PyArrayObject *py_triplets; - PyArrayObject *py_bz_grid_address; - PyArrayObject *py_bz_map; - PyArrayObject *py_map_triplets; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - long grid_point; - long bz_grid_type; - - long(*triplets)[3]; - long(*bz_grid_address)[3]; - long *bz_map; - long *map_triplets; - long num_map_triplets; - long *D_diag; - long(*Q)[3]; - long num_ir; - - if (!PyArg_ParseTuple(args, "OlOOOOOl", &py_triplets, &grid_point, - &py_bz_grid_address, &py_bz_map, &py_map_triplets, - &py_D_diag, &py_Q, &bz_grid_type)) { - return NULL; - } - - triplets = (long(*)[3])PyArray_DATA(py_triplets); - bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address); - bz_map = (long *)PyArray_DATA(py_bz_map); - map_triplets = (long *)PyArray_DATA(py_map_triplets); - num_map_triplets = (long)PyArray_DIMS(py_map_triplets)[0]; - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - - num_ir = ph3py_get_BZ_triplets_at_q(triplets, grid_point, bz_grid_address, - bz_map, map_triplets, num_map_triplets, - D_diag, Q, bz_grid_type); - - return PyLong_FromLong(num_ir); -} - -static PyObject *py_get_triplets_integration_weights(PyObject *self, - PyObject *args) { - PyArrayObject *py_iw; - PyArrayObject *py_iw_zero; - PyArrayObject *py_frequency_points; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_D_diag; - PyArrayObject *py_triplets; - PyArrayObject *py_frequencies1; - PyArrayObject *py_frequencies2; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_bz_map; - long bz_grid_type; - long tp_type; - - double *iw; - char *iw_zero; - double *frequency_points; - long(*relative_grid_address)[4][3]; - long *D_diag; - long(*triplets)[3]; - long(*bz_grid_addresses)[3]; - long *bz_map; - double *frequencies1, *frequencies2; - long num_band0, num_band1, num_band2, num_triplets; - - if (!PyArg_ParseTuple(args, "OOOOOOOOOOll", &py_iw, &py_iw_zero, - &py_frequency_points, &py_relative_grid_address, - &py_D_diag, &py_triplets, &py_frequencies1, - &py_frequencies2, &py_bz_grid_addresses, &py_bz_map, - &bz_grid_type, &tp_type)) { - return NULL; - } - - iw = (double *)PyArray_DATA(py_iw); - iw_zero = (char *)PyArray_DATA(py_iw_zero); - frequency_points = (double *)PyArray_DATA(py_frequency_points); - num_band0 = (long)PyArray_DIMS(py_frequency_points)[0]; - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - D_diag = (long *)PyArray_DATA(py_D_diag); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - bz_map = (long *)PyArray_DATA(py_bz_map); - frequencies1 = (double *)PyArray_DATA(py_frequencies1); - frequencies2 = (double *)PyArray_DATA(py_frequencies2); - num_band1 = (long)PyArray_DIMS(py_frequencies1)[1]; - num_band2 = (long)PyArray_DIMS(py_frequencies2)[1]; - - ph3py_get_integration_weight( - iw, iw_zero, frequency_points, num_band0, relative_grid_address, D_diag, - triplets, num_triplets, bz_grid_addresses, bz_map, bz_grid_type, - frequencies1, num_band1, frequencies2, num_band2, tp_type, 1); - - Py_RETURN_NONE; -} - -static PyObject *py_get_triplets_integration_weights_with_sigma( - PyObject *self, PyObject *args) { - PyArrayObject *py_iw; - PyArrayObject *py_iw_zero; - PyArrayObject *py_frequency_points; - PyArrayObject *py_triplets; - PyArrayObject *py_frequencies; - double sigma, sigma_cutoff; - - double *iw; - char *iw_zero; - double *frequency_points; - long(*triplets)[3]; - double *frequencies; - long num_band0, num_band, num_iw, num_triplets; - - if (!PyArg_ParseTuple(args, "OOOOOdd", &py_iw, &py_iw_zero, - &py_frequency_points, &py_triplets, &py_frequencies, - &sigma, &sigma_cutoff)) { - return NULL; - } - - iw = (double *)PyArray_DATA(py_iw); - iw_zero = (char *)PyArray_DATA(py_iw_zero); - frequency_points = (double *)PyArray_DATA(py_frequency_points); - num_band0 = (long)PyArray_DIMS(py_frequency_points)[0]; - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - frequencies = (double *)PyArray_DATA(py_frequencies); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - num_iw = (long)PyArray_DIMS(py_iw)[0]; - - ph3py_get_integration_weight_with_sigma( - iw, iw_zero, sigma, sigma_cutoff, frequency_points, num_band0, triplets, - num_triplets, frequencies, num_band, num_iw); - - Py_RETURN_NONE; -} - -static PyObject *py_get_grid_index_from_address(PyObject *self, - PyObject *args) { - PyArrayObject *py_address; - PyArrayObject *py_D_diag; - - long *address; - long *D_diag; - long gp; - - if (!PyArg_ParseTuple(args, "OO", &py_address, &py_D_diag)) { - return NULL; - } - - address = (long *)PyArray_DATA(py_address); - D_diag = (long *)PyArray_DATA(py_D_diag); - - gp = ph3py_get_grid_index_from_address(address, D_diag); - - return PyLong_FromLong(gp); -} - -static PyObject *py_get_gr_grid_addresses(PyObject *self, PyObject *args) { - PyArrayObject *py_gr_grid_addresses; - PyArrayObject *py_D_diag; - - long(*gr_grid_addresses)[3]; - long *D_diag; - - if (!PyArg_ParseTuple(args, "OO", &py_gr_grid_addresses, &py_D_diag)) { - return NULL; - } - - gr_grid_addresses = (long(*)[3])PyArray_DATA(py_gr_grid_addresses); - D_diag = (long *)PyArray_DATA(py_D_diag); - - ph3py_get_gr_grid_addresses(gr_grid_addresses, D_diag); - - Py_RETURN_NONE; -} - -static PyObject *py_get_reciprocal_rotations(PyObject *self, PyObject *args) { - PyArrayObject *py_rec_rotations; - PyArrayObject *py_rotations; - long is_time_reversal; - - long(*rec_rotations)[3][3]; - long(*rotations)[3][3]; - long num_rot, num_rec_rot; - - if (!PyArg_ParseTuple(args, "OOl", &py_rec_rotations, &py_rotations, - &is_time_reversal)) { - return NULL; - } - - rec_rotations = (long(*)[3][3])PyArray_DATA(py_rec_rotations); - rotations = (long(*)[3][3])PyArray_DATA(py_rotations); - num_rot = (long)PyArray_DIMS(py_rotations)[0]; - - num_rec_rot = ph3py_get_reciprocal_rotations(rec_rotations, rotations, - num_rot, is_time_reversal); - - return PyLong_FromLong(num_rec_rot); -} - -static PyObject *py_transform_rotations(PyObject *self, PyObject *args) { - PyArrayObject *py_transformed_rotations; - PyArrayObject *py_rotations; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - - long(*transformed_rotations)[3][3]; - long(*rotations)[3][3]; - long *D_diag; - long(*Q)[3]; - long num_rot, succeeded; - - if (!PyArg_ParseTuple(args, "OOOO", &py_transformed_rotations, - &py_rotations, &py_D_diag, &py_Q)) { - return NULL; - } - - transformed_rotations = - (long(*)[3][3])PyArray_DATA(py_transformed_rotations); - rotations = (long(*)[3][3])PyArray_DATA(py_rotations); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - num_rot = (long)PyArray_DIMS(py_transformed_rotations)[0]; - - succeeded = ph3py_transform_rotations(transformed_rotations, rotations, - num_rot, D_diag, Q); - if (succeeded) { - Py_RETURN_TRUE; - } else { - Py_RETURN_FALSE; - } -} - -static PyObject *py_get_snf3x3(PyObject *self, PyObject *args) { - PyArrayObject *py_D_diag; - PyArrayObject *py_P; - PyArrayObject *py_Q; - PyArrayObject *py_A; - - long *D_diag; - long(*P)[3]; - long(*Q)[3]; - long(*A)[3]; - long succeeded; - - if (!PyArg_ParseTuple(args, "OOOO", &py_D_diag, &py_P, &py_Q, &py_A)) { - return NULL; - } - - D_diag = (long *)PyArray_DATA(py_D_diag); - P = (long(*)[3])PyArray_DATA(py_P); - Q = (long(*)[3])PyArray_DATA(py_Q); - A = (long(*)[3])PyArray_DATA(py_A); - - succeeded = ph3py_get_snf3x3(D_diag, P, Q, A); - if (succeeded) { - Py_RETURN_TRUE; - } else { - Py_RETURN_FALSE; - } -} - -static PyObject *py_get_ir_grid_map(PyObject *self, PyObject *args) { - PyArrayObject *py_grid_mapping_table; - PyArrayObject *py_D_diag; - PyArrayObject *py_is_shift; - PyArrayObject *py_rotations; - - long *D_diag; - long *is_shift; - long(*rot)[3][3]; - long num_rot; - - long *grid_mapping_table; - long num_ir; - - if (!PyArg_ParseTuple(args, "OOOO", &py_grid_mapping_table, &py_D_diag, - &py_is_shift, &py_rotations)) { - return NULL; - } - - D_diag = (long *)PyArray_DATA(py_D_diag); - is_shift = (long *)PyArray_DATA(py_is_shift); - rot = (long(*)[3][3])PyArray_DATA(py_rotations); - num_rot = (long)PyArray_DIMS(py_rotations)[0]; - grid_mapping_table = (long *)PyArray_DATA(py_grid_mapping_table); - - num_ir = ph3py_get_ir_grid_map(grid_mapping_table, D_diag, is_shift, rot, - num_rot); - return PyLong_FromLong(num_ir); -} - -static PyObject *py_get_bz_grid_addresses(PyObject *self, PyObject *args) { - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_bz_map; - PyArrayObject *py_bzg2grg; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_PS; - PyArrayObject *py_reciprocal_lattice; - long type; - - long(*bz_grid_addresses)[3]; - long *bz_map; - long *bzg2grg; - long *D_diag; - long(*Q)[3]; - long *PS; - double(*reciprocal_lattice)[3]; - long num_total_gp; - - if (!PyArg_ParseTuple(args, "OOOOOOOl", &py_bz_grid_addresses, &py_bz_map, - &py_bzg2grg, &py_D_diag, &py_Q, &py_PS, - &py_reciprocal_lattice, &type)) { - return NULL; - } - - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - bz_map = (long *)PyArray_DATA(py_bz_map); - bzg2grg = (long *)PyArray_DATA(py_bzg2grg); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - PS = (long *)PyArray_DATA(py_PS); - reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice); - - num_total_gp = - ph3py_get_bz_grid_addresses(bz_grid_addresses, bz_map, bzg2grg, D_diag, - Q, PS, reciprocal_lattice, type); - - return PyLong_FromLong(num_total_gp); -} - -static PyObject *py_rotate_bz_grid_addresses(PyObject *self, PyObject *args) { - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_rotation; - PyArrayObject *py_bz_map; - PyArrayObject *py_D_diag; - PyArrayObject *py_PS; - long bz_grid_index; - long type; - - long(*bz_grid_addresses)[3]; - long(*rotation)[3]; - long *bz_map; - long *D_diag; - long *PS; - long ret_bz_gp; - - if (!PyArg_ParseTuple(args, "lOOOOOl", &bz_grid_index, &py_rotation, - &py_bz_grid_addresses, &py_bz_map, &py_D_diag, &py_PS, - &type)) { - return NULL; - } - - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - rotation = (long(*)[3])PyArray_DATA(py_rotation); - bz_map = (long *)PyArray_DATA(py_bz_map); - D_diag = (long *)PyArray_DATA(py_D_diag); - PS = (long *)PyArray_DATA(py_PS); - - ret_bz_gp = ph3py_rotate_bz_grid_index( - bz_grid_index, rotation, bz_grid_addresses, bz_map, D_diag, PS, type); - - return PyLong_FromLong(ret_bz_gp); -} - -static PyObject *py_diagonalize_collision_matrix(PyObject *self, - PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_eigenvalues; - double cutoff; - long i_sigma, i_temp, is_pinv, solver; - - double *collision_matrix; - double *eigvals; - long num_temp, num_grid_point, num_band; - long num_column, adrs_shift; - long info; - - if (!PyArg_ParseTuple(args, "OOlldll", &py_collision_matrix, - &py_eigenvalues, &i_sigma, &i_temp, &cutoff, &solver, - &is_pinv)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - eigvals = (double *)PyArray_DATA(py_eigenvalues); - - if (PyArray_NDIM(py_collision_matrix) == 2) { - num_temp = 1; - num_column = (long)PyArray_DIM(py_collision_matrix, 1); - } else { - num_temp = (long)PyArray_DIM(py_collision_matrix, 1); - num_grid_point = (long)PyArray_DIM(py_collision_matrix, 2); - num_band = (long)PyArray_DIM(py_collision_matrix, 3); - if (PyArray_NDIM(py_collision_matrix) == 8) { - num_column = num_grid_point * num_band * 3; - } else { - num_column = num_grid_point * num_band; - } - } - adrs_shift = (i_sigma * num_column * num_column * num_temp + - i_temp * num_column * num_column); - - /* show_colmat_info(py_collision_matrix, i_sigma, i_temp, adrs_shift); */ - - info = ph3py_phonopy_dsyev(collision_matrix + adrs_shift, eigvals, - num_column, solver); - if (is_pinv) { - ph3py_pinv_from_eigensolution(collision_matrix + adrs_shift, eigvals, - num_column, cutoff, 0); - } - - return PyLong_FromLong(info); -} - -static PyObject *py_pinv_from_eigensolution(PyObject *self, PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_eigenvalues; - double cutoff; - long i_sigma, i_temp, pinv_method; - - double *collision_matrix; - double *eigvals; - long num_temp, num_grid_point, num_band; - long num_column, adrs_shift; - - if (!PyArg_ParseTuple(args, "OOlldl", &py_collision_matrix, &py_eigenvalues, - &i_sigma, &i_temp, &cutoff, &pinv_method)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - eigvals = (double *)PyArray_DATA(py_eigenvalues); - num_temp = (long)PyArray_DIMS(py_collision_matrix)[1]; - num_grid_point = (long)PyArray_DIMS(py_collision_matrix)[2]; - num_band = (long)PyArray_DIMS(py_collision_matrix)[3]; - - if (PyArray_NDIM(py_collision_matrix) == 8) { - num_column = num_grid_point * num_band * 3; - } else { - num_column = num_grid_point * num_band; - } - adrs_shift = (i_sigma * num_column * num_column * num_temp + - i_temp * num_column * num_column); - - /* show_colmat_info(py_collision_matrix, i_sigma, i_temp, adrs_shift); */ - - ph3py_pinv_from_eigensolution(collision_matrix + adrs_shift, eigvals, - num_column, cutoff, pinv_method); - - Py_RETURN_NONE; -} - -static PyObject *py_get_default_colmat_solver(PyObject *self, PyObject *args) { - if (!PyArg_ParseTuple(args, "")) { - return NULL; - } - -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) - return PyLong_FromLong((long)1); -#else - return PyLong_FromLong((long)4); -#endif -} - -static PyObject *py_lapacke_pinv(PyObject *self, PyObject *args) { - PyArrayObject *data_in_py; - PyArrayObject *data_out_py; - double cutoff; - - int m; - int n; - double *data_in; - double *data_out; - int info; - - if (!PyArg_ParseTuple(args, "OOd", &data_out_py, &data_in_py, &cutoff)) { - return NULL; - } - - m = (long)PyArray_DIMS(data_in_py)[0]; - n = (long)PyArray_DIMS(data_in_py)[1]; - data_in = (double *)PyArray_DATA(data_in_py); - data_out = (double *)PyArray_DATA(data_out_py); - - info = ph3py_phonopy_pinv(data_out, data_in, m, n, cutoff); - - return PyLong_FromLong((long)info); -} - -static PyObject *py_get_omp_max_threads(PyObject *self, PyObject *args) { - return PyLong_FromLong(ph3py_get_max_threads()); -} - -static void show_colmat_info(const PyArrayObject *py_collision_matrix, - const long i_sigma, const long i_temp, - const long adrs_shift) { - long i; - - printf(" Array_shape:("); - for (i = 0; i < PyArray_NDIM(py_collision_matrix); i++) { - printf("%d", (int)PyArray_DIM(py_collision_matrix, i)); - if (i < PyArray_NDIM(py_collision_matrix) - 1) { - printf(","); - } else { - printf("), "); - } - } - printf("Data shift:%lu [%lu, %lu]\n", adrs_shift, i_sigma, i_temp); -} - -/** - * @brief Convert numpy "int_" array to phono3py long array structure. - * - * @param npyary - * @return Larray* - */ -static Larray *convert_to_larray(PyArrayObject *npyary) { - long i; - Larray *ary; - - ary = (Larray *)malloc(sizeof(Larray)); - for (i = 0; i < PyArray_NDIM(npyary); i++) { - ary->dims[i] = PyArray_DIMS(npyary)[i]; - } - ary->data = (long *)PyArray_DATA(npyary); - return ary; -} - -/** - * @brief Convert numpy "double" array to phono3py double array structure. - * - * @param npyary - * @return Darray* - * @note PyArray_NDIM receives non-const (PyArrayObject *). - */ -static Darray *convert_to_darray(PyArrayObject *npyary) { - int i; - Darray *ary; - - ary = (Darray *)malloc(sizeof(Darray)); - for (i = 0; i < PyArray_NDIM(npyary); i++) { - ary->dims[i] = PyArray_DIMS(npyary)[i]; - } - ary->data = (double *)PyArray_DATA(npyary); - return ary; -} diff --git a/c/_phono3py.cpp b/c/_phono3py.cpp index 061631a5..a7017d66 100644 --- a/c/_phono3py.cpp +++ b/c/_phono3py.cpp @@ -1049,6 +1049,20 @@ long py_rotate_bz_grid_addresses(long bz_grid_index, nb::ndarray<> py_rotation, return ret_bz_gp; } +long py_get_default_colmat_solver() { +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) || !defined(NO_INCLUDE_LAPACKE) + return (long)1; +#else + return (long)4; +#endif +} + +long py_get_omp_max_threads() { return ph3py_get_max_threads(); } + +#ifdef NO_INCLUDE_LAPACKE +long py_include_lapacke() { return 0; } +#else +long py_include_lapacke() { return 1; } long py_diagonalize_collision_matrix(nb::ndarray<> py_collision_matrix, nb::ndarray<> py_eigenvalues, long i_sigma, long i_temp, double cutoff, long solver, @@ -1118,14 +1132,6 @@ void py_pinv_from_eigensolution(nb::ndarray<> py_collision_matrix, num_column, cutoff, pinv_method); } -long py_get_default_colmat_solver() { -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) - return (long)1; -#else - return (long)4; -#endif -} - long py_lapacke_pinv(nb::ndarray<> data_out_py, nb::ndarray<> data_in_py, double cutoff) { long m; @@ -1143,8 +1149,7 @@ long py_lapacke_pinv(nb::ndarray<> data_out_py, nb::ndarray<> data_in_py, return info; } - -long py_get_omp_max_threads() { return ph3py_get_max_threads(); } +#endif NB_MODULE(_phono3py, m) { m.def("interaction", &py_get_interaction); @@ -1179,9 +1184,12 @@ NB_MODULE(_phono3py, m) { m.def("triplets_integration_weights", &py_get_triplets_integration_weights); m.def("triplets_integration_weights_with_sigma", &py_get_triplets_integration_weights_with_sigma); + m.def("default_colmat_solver", &py_get_default_colmat_solver); + m.def("omp_max_threads", &py_get_omp_max_threads); + m.def("include_lapacke", &py_include_lapacke); +#ifndef NO_INCLUDE_LAPACKE m.def("diagonalize_collision_matrix", &py_diagonalize_collision_matrix); m.def("pinv_from_eigensolution", &py_pinv_from_eigensolution); - m.def("default_colmat_solver", &py_get_default_colmat_solver); m.def("lapacke_pinv", &py_lapacke_pinv); - m.def("omp_max_threads", &py_get_omp_max_threads); +#endif } diff --git a/c/_phononcalc.c b/c/_phononcalc.c deleted file mode 100644 index 941b70da..00000000 --- a/c/_phononcalc.c +++ /dev/null @@ -1,232 +0,0 @@ -/* Copyright (C) 2021 Atsushi Togo */ -/* All rights reserved. */ - -/* This file is part of phonopy. */ - -/* Redistribution and use in source and binary forms, with or without */ -/* modification, are permitted provided that the following conditions */ -/* are met: */ - -/* * Redistributions of source code must retain the above copyright */ -/* notice, this list of conditions and the following disclaimer. */ - -/* * Redistributions in binary form must reproduce the above copyright */ -/* notice, this list of conditions and the following disclaimer in */ -/* the documentation and/or other materials provided with the */ -/* distribution. */ - -/* * Neither the name of the phonopy project nor the names of its */ -/* contributors may be used to endorse or promote products derived */ -/* from this software without specific prior written permission. */ - -/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */ -/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */ -/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */ -/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */ -/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */ -/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */ -/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */ -/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */ -/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */ -/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */ -/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ -/* POSSIBILITY OF SUCH DAMAGE. */ - -#include -#include - -// #include "lapack_wrapper.h" -#include "phononcalc.h" - -static PyObject *py_get_phonons_at_gridpoints(PyObject *self, PyObject *args); - -struct module_state { - PyObject *error; -}; - -#if PY_MAJOR_VERSION >= 3 -#define GETSTATE(m) ((struct module_state *)PyModule_GetState(m)) -#else -#define GETSTATE(m) (&_state) -static struct module_state _state; -#endif - -static PyObject *error_out(PyObject *m) { - struct module_state *st = GETSTATE(m); - PyErr_SetString(st->error, "something bad happened"); - return NULL; -} - -static PyMethodDef _phononcalc_methods[] = { - {"error_out", (PyCFunction)error_out, METH_NOARGS, NULL}, - {"phonons_at_gridpoints", py_get_phonons_at_gridpoints, METH_VARARGS, - "Set phonons at grid points"}, - {NULL, NULL, 0, NULL}}; - -#if PY_MAJOR_VERSION >= 3 - -static int _phononcalc_traverse(PyObject *m, visitproc visit, void *arg) { - Py_VISIT(GETSTATE(m)->error); - return 0; -} - -static int _phononcalc_clear(PyObject *m) { - Py_CLEAR(GETSTATE(m)->error); - return 0; -} - -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, "_phononcalc", NULL, - sizeof(struct module_state), _phononcalc_methods, NULL, - _phononcalc_traverse, _phononcalc_clear, NULL}; - -#define INITERROR return NULL - -PyObject *PyInit__phononcalc(void) -#else -#define INITERROR return - -void init_phononcalc(void) -#endif -{ -#if PY_MAJOR_VERSION >= 3 - PyObject *module = PyModule_Create(&moduledef); -#else - PyObject *module = Py_InitModule("_phononcalc", _phononcalc_methods); -#endif - struct module_state *st; - - if (module == NULL) INITERROR; - st = GETSTATE(module); - - st->error = PyErr_NewException("_phononcalc.Error", NULL, NULL); - if (st->error == NULL) { - Py_DECREF(module); - INITERROR; - } - -#if PY_MAJOR_VERSION >= 3 - return module; -#endif -} - -static PyObject *py_get_phonons_at_gridpoints(PyObject *self, PyObject *args) { - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_phonon_done; - PyArrayObject *py_grid_points; - PyArrayObject *py_grid_address; - PyArrayObject *py_QDinv; - PyArrayObject *py_shortest_vectors_fc2; - PyArrayObject *py_multiplicity_fc2; - PyArrayObject *py_positions_fc2; - PyArrayObject *py_fc2; - PyArrayObject *py_masses_fc2; - PyArrayObject *py_p2s_map_fc2; - PyArrayObject *py_s2p_map_fc2; - PyArrayObject *py_reciprocal_lattice; - PyArrayObject *py_born_effective_charge; - PyArrayObject *py_q_direction; - PyArrayObject *py_dielectric_constant; - PyArrayObject *py_dd_q0; - PyArrayObject *py_G_list; - double nac_factor; - double unit_conversion_factor; - double lambda; - char *uplo; - - double(*born)[3][3]; - double(*dielectric)[3]; - double *q_dir; - double *freqs; - _lapack_complex_double *eigvecs; - char *phonon_done; - long *grid_points; - long(*grid_address)[3]; - double(*QDinv)[3]; - double *fc2; - double(*svecs_fc2)[3]; - long(*multi_fc2)[2]; - double(*positions_fc2)[3]; - double *masses_fc2; - long *p2s_fc2; - long *s2p_fc2; - double(*rec_lat)[3]; - double(*dd_q0)[2]; - double(*G_list)[3]; - long num_patom, num_satom, num_phonons, num_grid_points, num_G_points; - - if (!PyArg_ParseTuple( - args, "OOOOOOOOOOOOOdOOOOdOOds", &py_frequencies, &py_eigenvectors, - &py_phonon_done, &py_grid_points, &py_grid_address, &py_QDinv, - &py_fc2, &py_shortest_vectors_fc2, &py_multiplicity_fc2, - &py_positions_fc2, &py_masses_fc2, &py_p2s_map_fc2, &py_s2p_map_fc2, - &unit_conversion_factor, &py_born_effective_charge, - &py_dielectric_constant, &py_reciprocal_lattice, &py_q_direction, - &nac_factor, &py_dd_q0, &py_G_list, &lambda, &uplo)) { - return NULL; - } - - freqs = (double *)PyArray_DATA(py_frequencies); - eigvecs = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - phonon_done = (char *)PyArray_DATA(py_phonon_done); - grid_points = (long *)PyArray_DATA(py_grid_points); - grid_address = (long(*)[3])PyArray_DATA(py_grid_address); - QDinv = (double(*)[3])PyArray_DATA(py_QDinv); - fc2 = (double *)PyArray_DATA(py_fc2); - svecs_fc2 = (double(*)[3])PyArray_DATA(py_shortest_vectors_fc2); - multi_fc2 = (long(*)[2])PyArray_DATA(py_multiplicity_fc2); - masses_fc2 = (double *)PyArray_DATA(py_masses_fc2); - p2s_fc2 = (long *)PyArray_DATA(py_p2s_map_fc2); - s2p_fc2 = (long *)PyArray_DATA(py_s2p_map_fc2); - rec_lat = (double(*)[3])PyArray_DATA(py_reciprocal_lattice); - num_patom = (long)PyArray_DIMS(py_multiplicity_fc2)[1]; - num_satom = (long)PyArray_DIMS(py_multiplicity_fc2)[0]; - num_phonons = (long)PyArray_DIMS(py_frequencies)[0]; - num_grid_points = (long)PyArray_DIMS(py_grid_points)[0]; - if ((PyObject *)py_born_effective_charge == Py_None) { - born = NULL; - } else { - born = (double(*)[3][3])PyArray_DATA(py_born_effective_charge); - } - if ((PyObject *)py_dielectric_constant == Py_None) { - dielectric = NULL; - } else { - dielectric = (double(*)[3])PyArray_DATA(py_dielectric_constant); - } - if ((PyObject *)py_q_direction == Py_None) { - q_dir = NULL; - } else { - q_dir = (double *)PyArray_DATA(py_q_direction); - if (fabs(q_dir[0]) < 1e-10 && fabs(q_dir[1]) < 1e-10 && - fabs(q_dir[2]) < 1e-10) { - q_dir = NULL; - } - } - if ((PyObject *)py_dd_q0 == Py_None) { - dd_q0 = NULL; - } else { - dd_q0 = (double(*)[2])PyArray_DATA(py_dd_q0); - } - if ((PyObject *)py_G_list == Py_None) { - G_list = NULL; - num_G_points = 0; - } else { - G_list = (double(*)[3])PyArray_DATA(py_G_list); - num_G_points = (long)PyArray_DIMS(py_G_list)[0]; - } - if ((PyObject *)py_positions_fc2 == Py_None) { - positions_fc2 = NULL; - } else { - positions_fc2 = (double(*)[3])PyArray_DATA(py_positions_fc2); - } - - phcalc_get_phonons_at_gridpoints( - freqs, eigvecs, phonon_done, num_phonons, grid_points, num_grid_points, - grid_address, QDinv, fc2, svecs_fc2, multi_fc2, positions_fc2, - num_patom, num_satom, masses_fc2, p2s_fc2, s2p_fc2, - unit_conversion_factor, born, dielectric, rec_lat, q_dir, nac_factor, - dd_q0, G_list, num_G_points, lambda, uplo[0]); - - Py_RETURN_NONE; -} diff --git a/c/lapack_wrapper.c b/c/lapack_wrapper.c index 0b6704b9..166e2000 100644 --- a/c/lapack_wrapper.c +++ b/c/lapack_wrapper.c @@ -38,13 +38,16 @@ #define min(a, b) ((a) > (b) ? (b) : (a)) -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) -MKL_Complex16 lapack_make_complex_double(double re, double im) { - MKL_Complex16 z; +#if (defined(MKL_BLAS) || defined(SCIPY_MKL_H)) && !defined(NO_INCLUDE_LAPACKE) +lapack_complex_double lapack_make_complex_double(double re, double im) { + lapack_complex_double z; z.real = re; z.imag = im; return z; } +#endif + +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) || defined(NO_INCLUDE_LAPACKE) #ifndef LAPACKE_malloc #define LAPACKE_malloc(size) malloc(size) #endif @@ -53,6 +56,18 @@ MKL_Complex16 lapack_make_complex_double(double re, double im) { #endif #endif +lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, + const lapack_complex_double b) { + lapack_complex_double c; + c = lapack_make_complex_double( + lapack_complex_double_real(a) * lapack_complex_double_real(b) - + lapack_complex_double_imag(a) * lapack_complex_double_imag(b), + lapack_complex_double_imag(a) * lapack_complex_double_real(b) + + lapack_complex_double_real(a) * lapack_complex_double_imag(b)); + return c; +} + +#ifndef NO_INCLUDE_LAPACKE int phonopy_zheev(double *w, lapack_complex_double *a, const int n, const char uplo) { lapack_int info; @@ -186,17 +201,6 @@ int phonopy_dsyev(double *data, double *eigvals, const int size, return (int)info; } -lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, - const lapack_complex_double b) { - lapack_complex_double c; - c = lapack_make_complex_double( - lapack_complex_double_real(a) * lapack_complex_double_real(b) - - lapack_complex_double_imag(a) * lapack_complex_double_imag(b), - lapack_complex_double_imag(a) * lapack_complex_double_real(b) + - lapack_complex_double_real(a) * lapack_complex_double_imag(b)); - return c; -} - void pinv_from_eigensolution(double *data, const double *eigvals, const long size, const double cutoff, const long pinv_method) { @@ -284,3 +288,4 @@ void pinv_from_eigensolution(double *data, const double *eigvals, free(tmp_data); tmp_data = NULL; } +#endif diff --git a/c/lapack_wrapper.h b/c/lapack_wrapper.h index 13d398c7..b18db54a 100644 --- a/c/lapack_wrapper.h +++ b/c/lapack_wrapper.h @@ -35,16 +35,34 @@ #ifndef __lapack_wrapper_H__ #define __lapack_wrapper_H__ -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) +#ifdef NO_INCLUDE_LAPACKE +#include +#define lapack_complex_double double _Complex +#ifdef CMPLX +#define lapack_make_complex_double(re, im) CMPLX(re, im) +#else +#define lapack_make_complex_double(re, im) ((double _Complex)((re) + (im) * I)) +#endif +#define lapack_complex_double_real(z) (creal(z)) +#define lapack_complex_double_imag(z) (cimag(z)) +#endif + +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) #include #define lapack_complex_double MKL_Complex16 +MKL_Complex16 lapack_make_complex_double(double re, double im); #define lapack_complex_double_real(z) ((z).real) #define lapack_complex_double_imag(z) ((z).imag) -MKL_Complex16 lapack_make_complex_double(double re, double im); -#else +#endif + +#if !defined(MKL_BLAS) && !defined(SCIPY_MKL_H) && !defined(NO_INCLUDE_LAPACKE) #include #endif +lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, + const lapack_complex_double b); + +#ifndef NO_INCLUDE_LAPACKE int phonopy_zheev(double *w, lapack_complex_double *a, const int n, const char uplo); int phonopy_pinv(double *data_out, const double *data_in, const int m, @@ -56,10 +74,9 @@ void phonopy_pinv_mt(double *data_out, int *info_out, const double *data_in, int phonopy_dsyev(double *data, double *eigvals, const int size, const int algorithm); -lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, - const lapack_complex_double b); - void pinv_from_eigensolution(double *data, const double *eigvals, const long size, const double cutoff, const long pinv_method); #endif + +#endif diff --git a/c/phono3py.c b/c/phono3py.c index 4f7f625b..39e7050e 100644 --- a/c/phono3py.c +++ b/c/phono3py.c @@ -697,6 +697,15 @@ long ph3py_get_thm_integration_weights_at_grid_points( return 1; } +long ph3py_get_max_threads(void) { +#ifdef _OPENMP + return omp_get_max_threads(); +#else + return 0; +#endif +} + +#ifndef NO_INCLUDE_LAPACKE long ph3py_phonopy_dsyev(double *data, double *eigvals, const long size, const long algorithm) { return (long)phonopy_dsyev(data, eigvals, (int)size, (int)algorithm); @@ -712,11 +721,4 @@ void ph3py_pinv_from_eigensolution(double *data, const double *eigvals, const long pinv_method) { pinv_from_eigensolution(data, eigvals, size, cutoff, pinv_method); } - -long ph3py_get_max_threads(void) { -#ifdef _OPENMP - return omp_get_max_threads(); -#else - return 0; #endif -} diff --git a/c/phono3py.h b/c/phono3py.h index d6a2972f..54247e24 100644 --- a/c/phono3py.h +++ b/c/phono3py.h @@ -233,6 +233,9 @@ long ph3py_get_thm_integration_weights_at_grid_points( const long *grid_points, const long (*bz_grid_addresses)[3], const long *bz_map, const long bz_grid_type, const double *frequencies, const long *gp2irgp_map, const char function); +long ph3py_get_max_threads(void); + +#ifndef NO_INCLUDE_LAPACKE long ph3py_phonopy_dsyev(double *data, double *eigvals, const long size, const long algorithm); long ph3py_phonopy_pinv(double *data_out, const double *data_in, const long m, @@ -240,7 +243,7 @@ long ph3py_phonopy_pinv(double *data_out, const double *data_in, const long m, void ph3py_pinv_from_eigensolution(double *data, const double *eigvals, const long size, const double cutoff, const long pinv_method); -long ph3py_get_max_threads(void); +#endif #ifdef __cplusplus } diff --git a/c/phonon.c b/c/phonon.c index 9282093b..f4e8a36d 100644 --- a/c/phonon.c +++ b/c/phonon.c @@ -36,6 +36,7 @@ #include #include +#include #include #include "dynmat.h" @@ -66,16 +67,6 @@ static void get_gonze_undone_phonons( const double *q_direction, const double nac_factor, const double (*dd_q0)[2], const double (*G_list)[3], const long num_G_points, const double lambda, const char uplo); -static void get_phonons(lapack_complex_double *eigvecs, const double q[3], - const double *fc2, const double *masses, - const long *p2s, const long *s2p, - const long (*multi)[2], const long num_patom, - const long num_satom, const double (*svecs)[3], - const long is_nac, const double (*born)[3][3], - const double dielectric[3][3], - const double reciprocal_lattice[3][3], - const double *q_direction, const double nac_factor, - const double unit_conversion_factor); static void get_gonze_phonons( lapack_complex_double *eigvecs, const double q[3], const double *fc2, const double *masses, const long *p2s, const long *s2p, @@ -208,12 +199,13 @@ static void get_undone_phonons( } is_nac = needs_nac(born, grid_address, gp, q_direction); - get_phonons(eigenvectors + num_band * num_band * gp, q, fc2, masses_fc2, - p2s_fc2, s2p_fc2, multi_fc2, num_patom, num_satom, - svecs_fc2, is_nac, born, dielectric, reciprocal_lattice, - q_direction, nac_factor, unit_conversion_factor); + get_dynamical_matrix(eigenvectors + num_band * num_band * gp, q, fc2, + masses_fc2, p2s_fc2, s2p_fc2, multi_fc2, num_patom, + num_satom, svecs_fc2, is_nac, born, dielectric, + reciprocal_lattice, q_direction, nac_factor); } +#ifndef NO_INCLUDE_LAPACKE /* To avoid multithreaded BLAS in OpenMP loop */ #ifndef MULTITHREADED_BLAS #ifdef _OPENMP @@ -235,6 +227,7 @@ static void get_undone_phonons( unit_conversion_factor; } } +#endif } static void get_gonze_undone_phonons( @@ -274,6 +267,7 @@ static void get_gonze_undone_phonons( nac_factor, dd_q0, G_list, num_G_points, lambda); } +#ifndef NO_INCLUDE_LAPACKE /* To avoid multithreaded BLAS in OpenMP loop */ #ifndef MULTITHREADED_BLAS #ifdef _OPENMP @@ -295,22 +289,7 @@ static void get_gonze_undone_phonons( unit_conversion_factor; } } -} - -static void get_phonons(lapack_complex_double *eigvecs, const double q[3], - const double *fc2, const double *masses, - const long *p2s, const long *s2p, - const long (*multi)[2], const long num_patom, - const long num_satom, const double (*svecs)[3], - const long is_nac, const double (*born)[3][3], - const double dielectric[3][3], - const double reciprocal_lattice[3][3], - const double *q_direction, const double nac_factor, - const double unit_conversion_factor) { - /* Store dynamical matrix in eigvecs array. */ - get_dynamical_matrix(eigvecs, q, fc2, masses, p2s, s2p, multi, num_patom, - num_satom, svecs, is_nac, born, dielectric, - reciprocal_lattice, q_direction, nac_factor); +#endif } static void get_gonze_phonons( diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index d061999c..02d6ffbd 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -35,7 +35,7 @@ #include "reciprocal_to_normal.h" #ifdef MULTITHREADED_BLAS -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) #include #else #include diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index 2f699db2..804aac86 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -243,6 +243,8 @@ def start_phono3py(**argparse_control) -> tuple[argparse.Namespace, int]: max_threads = phono3c.omp_max_threads() if max_threads > 0: print(f"Compiled with OpenMP support (max {max_threads} threads).") + if phono3c.include_lapacke(): + print("Compiled with LAPACKE.") if argparse_control.get("load_phono3py_yaml", False): print("Running in phono3py.load mode.") diff --git a/phono3py/phonon/solver.py b/phono3py/phonon/solver.py index a694c5bd..84b0dfd8 100644 --- a/phono3py/phonon/solver.py +++ b/phono3py/phonon/solver.py @@ -73,6 +73,7 @@ def run_phonon_solver_c( 'U' or 'L' for lapack zheev solver. Default is 'L'. """ + import phono3py._phono3py as phono3c import phono3py._phononcalc as phononcalc ( @@ -124,6 +125,9 @@ def run_phonon_solver_c( assert QDinv.flags.c_contiguous assert lapack_zheev_uplo in ("L", "U") + if not phono3c.include_lapacke(): + phonon_undone = np.where(phonon_done == 0)[0] + fc_p2s, fc_s2p = _get_fc_elements_mapping(dm, fc) phononcalc.phonons_at_gridpoints( frequencies, @@ -154,6 +158,15 @@ def run_phonon_solver_c( lapack_zheev_uplo, ) + if not phono3c.include_lapacke(): + for gp in phonon_undone: + frequencies[gp], eigenvectors[gp] = np.linalg.eigh(eigenvectors[gp]) + frequencies[phonon_undone] = ( + np.sign(frequencies[phonon_undone]) + * np.sqrt(np.abs(frequencies[phonon_undone])) + * frequency_conversion_factor + ) + def run_phonon_solver_py( grid_point, diff --git a/pyproject.toml b/pyproject.toml index 094ad8eb..50962076 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,6 +52,7 @@ sdist.include = [ PHONO3PY_USE_MTBLAS = {env="PHONO3PY_USE_MTBLAS", default="ON"} USE_CONDA_PATH = {env="USE_CONDA_PATH", default="ON"} PHONO3PY_USE_OMP = {env="PHONO3PY_USE_OMP", default="ON"} +BUILD_WITHOUT_LAPACKE = {env="BUILD_WITHOUT_LAPACKE", default="OFF"} [tool.setuptools_scm] write_to = "phono3py/_version.py" diff --git a/test/phonon3/test_fc3.py b/test/phonon3/test_fc3.py index 5f358277..6d831282 100644 --- a/test/phonon3/test_fc3.py +++ b/test/phonon3/test_fc3.py @@ -3,6 +3,7 @@ import numpy as np import pytest +import phono3py._phono3py as phono3c from phono3py import Phono3py from phono3py.phonon3.fc3 import ( cutoff_fc3_by_zero, @@ -257,40 +258,43 @@ def test_phonon_smat_alm_cutoff_fc3(si_pbesol_111_222_alm_cutoff_fc3: Phono3py): np.testing.assert_allclose(ph.fc2[0, 33], fc2_ref, atol=1e-6, rtol=0) -@pytest.mark.parametrize("pinv_solver", ["numpy", "lapacke"]) +@pytest.mark.skipif( + not phono3c.include_lapacke(), reason="requires to complile with lapacke" +) def test_fc3_lapacke_solver(si_pbesol_111: Phono3py, pinv_solver: str): """Test fc3 with Si PBEsol 1x1x1 using lapacke solver.""" - ph = si_pbesol_111 - _, fc3 = get_fc3( - ph.supercell, - ph.primitive, - ph.dataset, - ph.symmetry, - pinv_solver=pinv_solver, - verbose=True, - ) - set_translational_invariance_fc3(fc3) - set_permutation_symmetry_fc3(fc3) - - fc3_ref = [ - [ - [1.07250822e-01, 1.86302073e-17, -4.26452855e-18], - [8.96414569e-03, -1.43046911e-01, -1.38498937e-01], - [-8.96414569e-03, -1.38498937e-01, -1.43046911e-01], - ], - [ - [-8.96414569e-03, -1.43046911e-01, -1.38498937e-01], - [-3.39457157e-02, -4.63315728e-17, -4.17779237e-17], - [-3.31746167e-01, -2.60025724e-02, -2.60025724e-02], - ], - [ - [8.96414569e-03, -1.38498937e-01, -1.43046911e-01], - [-3.31746167e-01, 2.60025724e-02, 2.60025724e-02], - [-3.39457157e-02, 3.69351540e-17, 5.94504191e-18], - ], - ] - - np.testing.assert_allclose(fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0) + for pinv_solver in ["lapacke", "numpy"]: + ph = si_pbesol_111 + _, fc3 = get_fc3( + ph.supercell, + ph.primitive, + ph.dataset, + ph.symmetry, + pinv_solver=pinv_solver, + verbose=True, + ) + set_translational_invariance_fc3(fc3) + set_permutation_symmetry_fc3(fc3) + + fc3_ref = [ + [ + [1.07250822e-01, 1.86302073e-17, -4.26452855e-18], + [8.96414569e-03, -1.43046911e-01, -1.38498937e-01], + [-8.96414569e-03, -1.38498937e-01, -1.43046911e-01], + ], + [ + [-8.96414569e-03, -1.43046911e-01, -1.38498937e-01], + [-3.39457157e-02, -4.63315728e-17, -4.17779237e-17], + [-3.31746167e-01, -2.60025724e-02, -2.60025724e-02], + ], + [ + [8.96414569e-03, -1.38498937e-01, -1.43046911e-01], + [-3.31746167e-01, 2.60025724e-02, 2.60025724e-02], + [-3.39457157e-02, 3.69351540e-17, 5.94504191e-18], + ], + ] + + np.testing.assert_allclose(fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0) def test_fc3_symfc(si_pbesol_111_symfc: Phono3py):