scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Option to run without linking BLAS and LAPACK in C

This commit is contained in:
Atsushi Togo 2024-11-23 20:36:37 +09:00
parent 59dc78e982
commit 232f22318c
15 changed files with 323 additions and 2288 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
.github/workflows/phono3py-pytest-conda-nolapacke.yml vendored Normal file
View File

@ -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

217
CMakeLists.txt
View File

@ -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)
if(BUILD_WITHOUT_LAPACKE)
if(OpenMP_FOUND)
target_link_libraries(phphcalc_lib PRIVATE recgrid_lib
OpenMP::OpenMP_C)
else()
target_link_libraries(phphcalc_lib PRIVATE recgrid_lib)
endif()
target_compile_definitions(phphcalc_lib PRIVATE NO_INCLUDE_LAPACKE
THM_EPSILON=1e-10)
else()
target_link_libraries(phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS
LAPACK::LAPACK)
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()
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()
target_compile_definitions(phphcalc_lib PRIVATE MKL_LAPACKE
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()
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)
if(BUILD_WITHOUT_LAPACKE)
if(OpenMP_FOUND)
target_link_libraries(phphcalc_lib recgrid_lib OpenMP::OpenMP_C)
else()
target_link_libraries(phphcalc_lib recgrid_lib)
endif()
target_compile_definitions(phphcalc_lib PRIVATE NO_INCLUDE_LAPACKE
THM_EPSILON=1e-10)
else()
target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK)
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)
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})
if(BLAS_LIBRARIES MATCHES "libmkl")
if(PHONO3PY_USE_MTBLAS)
target_compile_definitions(
phphcalc_lib PRIVATE MKL_LAPACKE MULTITHREADED_BLAS
THM_EPSILON=1e-10)
else()
target_compile_definitions(phphcalc_lib PRIVATE MKL_LAPACKE
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()
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)
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()
if(BLAS_LIBRARIES MATCHES "libmkl")
target_compile_definitions(phononcalc_lib PRIVATE MKL_BLAS
MULTITHREADED_BLAS)
endif()
if(BLAS_LIBRARIES MATCHES "libopenblas")
target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS)
endif()
endif()
target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES})
if(BLAS_LIBRARIES MATCHES "libmkl")
target_compile_definitions(phononcalc_lib PRIVATE MKL_LAPACKE
MULTITHREADED_BLAS)
endif()
if(BLAS_LIBRARIES MATCHES "libopenblas")
target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS)
endif()
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)
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()
if(BLAS_LIBRARIES MATCHES "libmkl")
target_compile_definitions(phononcalc_lib PRIVATE MKL_BLAS
MULTITHREADED_BLAS)
endif()
if(BLAS_LIBRARIES MATCHES "libopenblas")
target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS)
endif()
endif()
target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES})
if(BLAS_LIBRARIES MATCHES "libmkl")
target_compile_definitions(phononcalc_lib PRIVATE MKL_LAPACKE
MULTITHREADED_BLAS)
endif()
if(BLAS_LIBRARIES MATCHES "libopenblas")
target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS)
endif()
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})

1877
c/_phono3py.c
View File

File diff suppressed because it is too large Load Diff

32
c/_phono3py.cpp
View File

@ -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
}

232
c/_phononcalc.c
View File

@ -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 <Python.h>
#include <numpy/arrayobject.h>
// #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;
}

33
c/lapack_wrapper.c
View File

@ -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

29
c/lapack_wrapper.h
View File

@ -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 <complex.h>
#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 <mkl.h>
#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 <lapacke.h>
#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

16
c/phono3py.c
View File

@ -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
}

5
c/phono3py.h
View File

@ -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
}

39
c/phonon.c
View File

@ -36,6 +36,7 @@
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#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(

2
c/reciprocal_to_normal.c
View File

@ -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 <mkl_cblas.h>
#else
#include <cblas.h>

2
phono3py/cui/phono3py_script.py
View File

@ -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.")

13
phono3py/phonon/solver.py
View File

@ -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,

1
pyproject.toml
View File

@ -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"

64
test/phonon3/test_fc3.py
View File

@ -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)
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],
],
]
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)
np.testing.assert_allclose(fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0)
def test_fc3_symfc(si_pbesol_111_symfc: Phono3py):