1 /* Copyright (C) 2015 Atsushi Togo */
2 /* All rights reserved. */
3
4 /* This file is part of phonopy. */
5
6 /* Redistribution and use in source and binary forms, with or without */
7 /* modification, are permitted provided that the following conditions */
8 /* are met: */
9
10 /* * Redistributions of source code must retain the above copyright */
11 /* notice, this list of conditions and the following disclaimer. */
12
13 /* * Redistributions in binary form must reproduce the above copyright */
14 /* notice, this list of conditions and the following disclaimer in */
15 /* the documentation and/or other materials provided with the */
16 /* distribution. */
17
18 /* * Neither the name of the phonopy project nor the names of its */
19 /* contributors may be used to endorse or promote products derived */
20 /* from this software without specific prior written permission. */
21
22 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
23 /* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
24 /* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
25 /* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
26 /* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
27 /* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
28 /* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
29 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
30 /* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
31 /* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
32 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
33 /* POSSIBILITY OF SUCH DAMAGE. */
34
35 #include <Python.h>
36 #include <assert.h>
37 #include <stdio.h>
38 #include <stddef.h>
39 #include <stdlib.h>
40 #include <math.h>
41 #include <numpy/arrayobject.h>
42 #include "lapack_wrapper.h"
43 #include "phono3py.h"
44 #include "phonoc_array.h"
45
46
47 static PyObject * py_get_interaction(PyObject *self, PyObject *args);
48 static PyObject * py_get_pp_collision(PyObject *self, PyObject *args);
49 static PyObject *
50 py_get_pp_collision_with_sigma(PyObject *self, PyObject *args);
51 static PyObject *
52 py_get_imag_self_energy_with_g(PyObject *self, PyObject *args);
53 static PyObject *
54 py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args);
55 static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
56 PyObject *args);
57 static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
58 PyObject *args);
59 static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args);
60 static PyObject * py_get_reducible_collision_matrix(PyObject *self,
61 PyObject *args);
62 static PyObject * py_symmetrize_collision_matrix(PyObject *self,
63 PyObject *args);
64 static PyObject * py_expand_collision_matrix(PyObject *self, PyObject *args);
65 static PyObject * py_distribute_fc3(PyObject *self, PyObject *args);
66 static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args);
67 static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args);
68 static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args);
69 static PyObject *
70 py_set_permutation_symmetry_fc3(PyObject *self, PyObject *args);
71 static PyObject *
72 py_set_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args);
73 static PyObject * py_set_permutation_symmetry_fc3(PyObject *self,
74 PyObject *args);
75 static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args);
76 static PyObject * py_get_neighboring_gird_points(PyObject *self, PyObject *args);
77 static PyObject * py_set_integration_weights(PyObject *self, PyObject *args);
78 static PyObject *
79 py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args);
80 static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args);
81 static PyObject *
82 py_set_triplets_integration_weights(PyObject *self, PyObject *args);
83 static PyObject *
84 py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args);
85 static PyObject *
86 py_diagonalize_collision_matrix(PyObject *self, PyObject *args);
87 static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args);
88 static PyObject * py_get_default_colmat_solver(PyObject *self, PyObject *args);
89
90 static void pinv_from_eigensolution(double *data,
91 const double *eigvals,
92 const size_t size,
93 const double cutoff,
94 const int pinv_method);
95 static void show_colmat_info(const PyArrayObject *collision_matrix_py,
96 const size_t i_sigma,
97 const size_t i_temp,
98 const size_t adrs_shift);
99 static Iarray* convert_to_iarray(const PyArrayObject* npyary);
100 static Darray* convert_to_darray(const PyArrayObject* npyary);
101
102
103 struct module_state {
104 PyObject *error;
105 };
106
107 #if PY_MAJOR_VERSION >= 3
108 #define GETSTATE(m) ((struct module_state*)PyModule_GetState(m))
109 #else
110 #define GETSTATE(m) (&_state)
111 static struct module_state _state;
112 #endif
113
114 static PyObject *
error_out(PyObject * m)115 error_out(PyObject *m) {
116 struct module_state *st = GETSTATE(m);
117 PyErr_SetString(st->error, "something bad happened");
118 return NULL;
119 }
120
121 static PyMethodDef _phono3py_methods[] = {
122 {"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
123 {"interaction",
124 (PyCFunction)py_get_interaction,
125 METH_VARARGS,
126 "Interaction of triplets"},
127 {"pp_collision",
128 (PyCFunction)py_get_pp_collision,
129 METH_VARARGS,
130 "Collision and ph-ph calculation"},
131 {"pp_collision_with_sigma",
132 (PyCFunction)py_get_pp_collision_with_sigma,
133 METH_VARARGS,
134 "Collision and ph-ph calculation for smearing method"},
135 {"imag_self_energy_with_g",
136 (PyCFunction)py_get_imag_self_energy_with_g,
137 METH_VARARGS,
138 "Imaginary part of self energy at frequency points with g"},
139 {"detailed_imag_self_energy_with_g",
140 (PyCFunction)py_get_detailed_imag_self_energy_with_g,
141 METH_VARARGS,
142 "Detailed contribution to imaginary part of self energy at frequency points with g"},
143 {"real_self_energy_at_bands",
144 (PyCFunction)py_get_real_self_energy_at_bands,
145 METH_VARARGS,
146 "Real part of self energy from third order force constants"},
147 {"real_self_energy_at_frequency_point",
148 (PyCFunction)py_get_real_self_energy_at_frequency_point,
149 METH_VARARGS,
150 "Real part of self energy from third order force constants at a frequency point"},
151 {"collision_matrix",
152 (PyCFunction)py_get_collision_matrix,
153 METH_VARARGS,
154 "Collision matrix with g"},
155 {"reducible_collision_matrix",
156 (PyCFunction)py_get_reducible_collision_matrix,
157 METH_VARARGS,
158 "Collision matrix with g for reducible grid points"},
159 {"symmetrize_collision_matrix",
160 (PyCFunction)py_symmetrize_collision_matrix,
161 METH_VARARGS,
162 "Symmetrize collision matrix"},
163 {"expand_collision_matrix",
164 (PyCFunction)py_expand_collision_matrix,
165 METH_VARARGS,
166 "Expand collision matrix"},
167 {"distribute_fc3",
168 (PyCFunction)py_distribute_fc3,
169 METH_VARARGS,
170 "Distribute least fc3 to full fc3"},
171 {"rotate_delta_fc2s",
172 (PyCFunction)py_rotate_delta_fc2s,
173 METH_VARARGS,
174 "Rotate delta fc2s"},
175 {"isotope_strength",
176 (PyCFunction)py_get_isotope_strength,
177 METH_VARARGS,
178 "Isotope scattering strength"},
179 {"thm_isotope_strength",
180 (PyCFunction)py_get_thm_isotope_strength,
181 METH_VARARGS,
182 "Isotope scattering strength for tetrahedron_method"},
183 {"permutation_symmetry_fc3",
184 (PyCFunction)py_set_permutation_symmetry_fc3,
185 METH_VARARGS,
186 "Set permutation symmetry for fc3"},
187 {"permutation_symmetry_compact_fc3",
188 (PyCFunction)py_set_permutation_symmetry_compact_fc3,
189 METH_VARARGS,
190 "Set permutation symmetry for compact-fc3"},
191 {"transpose_compact_fc3",
192 (PyCFunction)py_transpose_compact_fc3,
193 METH_VARARGS,
194 "Transpose compact fc3"},
195 {"neighboring_grid_points",
196 (PyCFunction)py_get_neighboring_gird_points,
197 METH_VARARGS,
198 "Neighboring grid points by relative grid addresses"},
199 {"integration_weights",
200 (PyCFunction)py_set_integration_weights,
201 METH_VARARGS,
202 "Integration weights of tetrahedron method"},
203 {"triplets_reciprocal_mesh_at_q",
204 (PyCFunction)py_tpl_get_triplets_reciprocal_mesh_at_q,
205 METH_VARARGS,
206 "Triplets on reciprocal mesh points at a specific q-point"},
207 {"BZ_triplets_at_q",
208 (PyCFunction)py_tpl_get_BZ_triplets_at_q,
209 METH_VARARGS,
210 "Triplets in reciprocal primitive lattice are transformed to those in BZ."},
211 {"triplets_integration_weights",
212 (PyCFunction)py_set_triplets_integration_weights,
213 METH_VARARGS,
214 "Integration weights of tetrahedron method for triplets"},
215 {"triplets_integration_weights_with_sigma",
216 (PyCFunction)py_set_triplets_integration_weights_with_sigma,
217 METH_VARARGS,
218 "Integration weights of smearing method for triplets"},
219 {"diagonalize_collision_matrix",
220 (PyCFunction)py_diagonalize_collision_matrix,
221 METH_VARARGS,
222 "Diagonalize and optionally pseudo-inverse using Lapack dsyev(d)"},
223 {"pinv_from_eigensolution",
224 (PyCFunction)py_pinv_from_eigensolution,
225 METH_VARARGS,
226 "Pseudo-inverse from eigensolution"},
227 {"default_colmat_solver",
228 (PyCFunction)py_get_default_colmat_solver,
229 METH_VARARGS,
230 "Return default collison matrix solver by integer value"},
231 {NULL, NULL, 0, NULL}
232 };
233
234 #if PY_MAJOR_VERSION >= 3
235
_phono3py_traverse(PyObject * m,visitproc visit,void * arg)236 static int _phono3py_traverse(PyObject *m, visitproc visit, void *arg) {
237 Py_VISIT(GETSTATE(m)->error);
238 return 0;
239 }
240
_phono3py_clear(PyObject * m)241 static int _phono3py_clear(PyObject *m) {
242 Py_CLEAR(GETSTATE(m)->error);
243 return 0;
244 }
245
246 static struct PyModuleDef moduledef = {
247 PyModuleDef_HEAD_INIT,
248 "_phono3py",
249 NULL,
250 sizeof(struct module_state),
251 _phono3py_methods,
252 NULL,
253 _phono3py_traverse,
254 _phono3py_clear,
255 NULL
256 };
257
258 #define INITERROR return NULL
259
260 PyObject *
PyInit__phono3py(void)261 PyInit__phono3py(void)
262
263 #else
264 #define INITERROR return
265
266 void
267 init_phono3py(void)
268 #endif
269 {
270 #if PY_MAJOR_VERSION >= 3
271 PyObject *module = PyModule_Create(&moduledef);
272 #else
273 PyObject *module = Py_InitModule("_phono3py", _phono3py_methods);
274 #endif
275 struct module_state *st;
276
277 if (module == NULL)
278 INITERROR;
279 st = GETSTATE(module);
280
281 st->error = PyErr_NewException("_phono3py.Error", NULL, NULL);
282 if (st->error == NULL) {
283 Py_DECREF(module);
284 INITERROR;
285 }
286
287 #if PY_MAJOR_VERSION >= 3
288 return module;
289 #endif
290 }
291
py_get_interaction(PyObject * self,PyObject * args)292 static PyObject * py_get_interaction(PyObject *self, PyObject *args)
293 {
294 PyArrayObject *py_fc3_normal_squared;
295 PyArrayObject *py_g_zero;
296 PyArrayObject *py_frequencies;
297 PyArrayObject *py_eigenvectors;
298 PyArrayObject *py_triplets;
299 PyArrayObject *py_grid_address;
300 PyArrayObject *py_mesh;
301 PyArrayObject *py_shortest_vectors;
302 PyArrayObject *py_multiplicities;
303 PyArrayObject *py_fc3;
304 PyArrayObject *py_masses;
305 PyArrayObject *py_p2s_map;
306 PyArrayObject *py_s2p_map;
307 PyArrayObject *py_band_indices;
308 double cutoff_frequency;
309 int symmetrize_fc3_q;
310
311 Darray *fc3_normal_squared;
312 Darray *freqs;
313 lapack_complex_double *eigvecs;
314 size_t (*triplets)[3];
315 npy_intp num_triplets;
316 char* g_zero;
317 int *grid_address;
318 int *mesh;
319 double *fc3;
320 double *svecs;
321 int *multi;
322 double *masses;
323 int *p2s;
324 int *s2p;
325 int *band_indices;
326 int svecs_dims[3];
327 int i;
328 int is_compact_fc3;
329
330 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOid",
331 &py_fc3_normal_squared,
332 &py_g_zero,
333 &py_frequencies,
334 &py_eigenvectors,
335 &py_triplets,
336 &py_grid_address,
337 &py_mesh,
338 &py_fc3,
339 &py_shortest_vectors,
340 &py_multiplicities,
341 &py_masses,
342 &py_p2s_map,
343 &py_s2p_map,
344 &py_band_indices,
345 &symmetrize_fc3_q,
346 &cutoff_frequency)) {
347 return NULL;
348 }
349
350
351 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
352 freqs = convert_to_darray(py_frequencies);
353 /* npy_cdouble and lapack_complex_double may not be compatible. */
354 /* So eigenvectors should not be used in Python side */
355 eigvecs = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
356 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
357 num_triplets = PyArray_DIMS(py_triplets)[0];
358 g_zero = (char*)PyArray_DATA(py_g_zero);
359 grid_address = (int*)PyArray_DATA(py_grid_address);
360 mesh = (int*)PyArray_DATA(py_mesh);
361 fc3 = (double*)PyArray_DATA(py_fc3);
362 if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
363 is_compact_fc3 = 0;
364 } else {
365 is_compact_fc3 = 1;
366 }
367 svecs = (double*)PyArray_DATA(py_shortest_vectors);
368 for (i = 0; i < 3; i++) {
369 svecs_dims[i] = PyArray_DIMS(py_shortest_vectors)[i];
370 }
371 multi = (int*)PyArray_DATA(py_multiplicities);
372 masses = (double*)PyArray_DATA(py_masses);
373 p2s = (int*)PyArray_DATA(py_p2s_map);
374 s2p = (int*)PyArray_DATA(py_s2p_map);
375 band_indices = (int*)PyArray_DATA(py_band_indices);
376
377 ph3py_get_interaction(fc3_normal_squared,
378 g_zero,
379 freqs,
380 eigvecs,
381 triplets,
382 num_triplets,
383 grid_address,
384 mesh,
385 fc3,
386 is_compact_fc3,
387 svecs,
388 svecs_dims,
389 multi,
390 masses,
391 p2s,
392 s2p,
393 band_indices,
394 symmetrize_fc3_q,
395 cutoff_frequency);
396
397 free(fc3_normal_squared);
398 fc3_normal_squared = NULL;
399 free(freqs);
400 freqs = NULL;
401
402 Py_RETURN_NONE;
403 }
404
py_get_pp_collision(PyObject * self,PyObject * args)405 static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
406 {
407 PyArrayObject *py_gamma;
408 PyArrayObject *py_relative_grid_address;
409 PyArrayObject *py_frequencies;
410 PyArrayObject *py_eigenvectors;
411 PyArrayObject *py_triplets;
412 PyArrayObject *py_triplet_weights;
413 PyArrayObject *py_grid_address;
414 PyArrayObject *py_bz_map;
415 PyArrayObject *py_mesh;
416 PyArrayObject *py_fc3;
417 PyArrayObject *py_shortest_vectors;
418 PyArrayObject *py_multiplicities;
419 PyArrayObject *py_masses;
420 PyArrayObject *py_p2s_map;
421 PyArrayObject *py_s2p_map;
422 PyArrayObject *py_band_indices;
423 PyArrayObject *py_temperatures;
424 double cutoff_frequency;
425 int is_NU;
426 int symmetrize_fc3_q;
427
428 double *gamma;
429 int (*relative_grid_address)[4][3];
430 double *frequencies;
431 lapack_complex_double *eigenvectors;
432 size_t (*triplets)[3];
433 npy_intp num_triplets;
434 int *triplet_weights;
435 int *grid_address;
436 size_t *bz_map;
437 int *mesh;
438 double *fc3;
439 double *svecs;
440 int *multi;
441 double *masses;
442 int *p2s;
443 int *s2p;
444 Iarray *band_indices;
445 Darray *temperatures;
446 int svecs_dims[3];
447 int i;
448 int is_compact_fc3;
449
450 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOiid",
451 &py_gamma,
452 &py_relative_grid_address,
453 &py_frequencies,
454 &py_eigenvectors,
455 &py_triplets,
456 &py_triplet_weights,
457 &py_grid_address,
458 &py_bz_map,
459 &py_mesh,
460 &py_fc3,
461 &py_shortest_vectors,
462 &py_multiplicities,
463 &py_masses,
464 &py_p2s_map,
465 &py_s2p_map,
466 &py_band_indices,
467 &py_temperatures,
468 &is_NU,
469 &symmetrize_fc3_q,
470 &cutoff_frequency)) {
471 return NULL;
472 }
473
474 gamma = (double*)PyArray_DATA(py_gamma);
475 relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
476 frequencies = (double*)PyArray_DATA(py_frequencies);
477 eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
478 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
479 num_triplets = PyArray_DIMS(py_triplets)[0];
480 triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
481 grid_address = (int*)PyArray_DATA(py_grid_address);
482 bz_map = (size_t*)PyArray_DATA(py_bz_map);
483 mesh = (int*)PyArray_DATA(py_mesh);
484 fc3 = (double*)PyArray_DATA(py_fc3);
485 if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
486 is_compact_fc3 = 0;
487 } else {
488 is_compact_fc3 = 1;
489 }
490 svecs = (double*)PyArray_DATA(py_shortest_vectors);
491 for (i = 0; i < 3; i++) {
492 svecs_dims[i] = PyArray_DIMS(py_shortest_vectors)[i];
493 }
494 multi = (int*)PyArray_DATA(py_multiplicities);
495 masses = (double*)PyArray_DATA(py_masses);
496 p2s = (int*)PyArray_DATA(py_p2s_map);
497 s2p = (int*)PyArray_DATA(py_s2p_map);
498 band_indices = convert_to_iarray(py_band_indices);
499 temperatures = convert_to_darray(py_temperatures);
500
501 ph3py_get_pp_collision(gamma,
502 relative_grid_address,
503 frequencies,
504 eigenvectors,
505 triplets,
506 num_triplets,
507 triplet_weights,
508 grid_address,
509 bz_map,
510 mesh,
511 fc3,
512 is_compact_fc3,
513 svecs,
514 svecs_dims,
515 multi,
516 masses,
517 p2s,
518 s2p,
519 band_indices,
520 temperatures,
521 is_NU,
522 symmetrize_fc3_q,
523 cutoff_frequency);
524
525 free(band_indices);
526 band_indices = NULL;
527 free(temperatures);
528 temperatures = NULL;
529
530 Py_RETURN_NONE;
531 }
532
py_get_pp_collision_with_sigma(PyObject * self,PyObject * args)533 static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
534 {
535 PyArrayObject *py_gamma;
536 PyArrayObject *py_frequencies;
537 PyArrayObject *py_eigenvectors;
538 PyArrayObject *py_triplets;
539 PyArrayObject *py_triplet_weights;
540 PyArrayObject *py_grid_address;
541 PyArrayObject *py_mesh;
542 PyArrayObject *py_fc3;
543 PyArrayObject *py_shortest_vectors;
544 PyArrayObject *py_multiplicities;
545 PyArrayObject *py_masses;
546 PyArrayObject *py_p2s_map;
547 PyArrayObject *py_s2p_map;
548 PyArrayObject *py_band_indices;
549 PyArrayObject *py_temperatures;
550 int is_NU;
551 int symmetrize_fc3_q;
552 double sigma;
553 double sigma_cutoff;
554 double cutoff_frequency;
555
556 double *gamma;
557 double *frequencies;
558 lapack_complex_double *eigenvectors;
559 size_t (*triplets)[3];
560 npy_intp num_triplets;
561 int *triplet_weights;
562 int *grid_address;
563 int *mesh;
564 double *fc3;
565 double *svecs;
566 int *multi;
567 double *masses;
568 int *p2s;
569 int *s2p;
570 Iarray *band_indices;
571 Darray *temperatures;
572 int svecs_dims[3];
573 int i;
574 int is_compact_fc3;
575
576 if (!PyArg_ParseTuple(args, "OddOOOOOOOOOOOOOOiid",
577 &py_gamma,
578 &sigma,
579 &sigma_cutoff,
580 &py_frequencies,
581 &py_eigenvectors,
582 &py_triplets,
583 &py_triplet_weights,
584 &py_grid_address,
585 &py_mesh,
586 &py_fc3,
587 &py_shortest_vectors,
588 &py_multiplicities,
589 &py_masses,
590 &py_p2s_map,
591 &py_s2p_map,
592 &py_band_indices,
593 &py_temperatures,
594 &is_NU,
595 &symmetrize_fc3_q,
596 &cutoff_frequency)) {
597 return NULL;
598 }
599
600 gamma = (double*)PyArray_DATA(py_gamma);
601 frequencies = (double*)PyArray_DATA(py_frequencies);
602 eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
603 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
604 num_triplets = PyArray_DIMS(py_triplets)[0];
605 triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
606 grid_address = (int*)PyArray_DATA(py_grid_address);
607 mesh = (int*)PyArray_DATA(py_mesh);
608 fc3 = (double*)PyArray_DATA(py_fc3);
609 if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
610 is_compact_fc3 = 0;
611 } else {
612 is_compact_fc3 = 1;
613 }
614 svecs = (double*)PyArray_DATA(py_shortest_vectors);
615 for (i = 0; i < 3; i++) {
616 svecs_dims[i] = PyArray_DIMS(py_shortest_vectors)[i];
617 }
618 multi = (int*)PyArray_DATA(py_multiplicities);
619 masses = (double*)PyArray_DATA(py_masses);
620 p2s = (int*)PyArray_DATA(py_p2s_map);
621 s2p = (int*)PyArray_DATA(py_s2p_map);
622 band_indices = convert_to_iarray(py_band_indices);
623 temperatures = convert_to_darray(py_temperatures);
624
625 ph3py_get_pp_collision_with_sigma(gamma,
626 sigma,
627 sigma_cutoff,
628 frequencies,
629 eigenvectors,
630 triplets,
631 num_triplets,
632 triplet_weights,
633 grid_address,
634 mesh,
635 fc3,
636 is_compact_fc3,
637 svecs,
638 svecs_dims,
639 multi,
640 masses,
641 p2s,
642 s2p,
643 band_indices,
644 temperatures,
645 is_NU,
646 symmetrize_fc3_q,
647 cutoff_frequency);
648
649 free(band_indices);
650 band_indices = NULL;
651 free(temperatures);
652 temperatures = NULL;
653
654 Py_RETURN_NONE;
655 }
656
py_get_imag_self_energy_with_g(PyObject * self,PyObject * args)657 static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
658 {
659 PyArrayObject *py_gamma;
660 PyArrayObject *py_fc3_normal_squared;
661 PyArrayObject *py_frequencies;
662 PyArrayObject *py_triplets;
663 PyArrayObject *py_triplet_weights;
664 PyArrayObject *py_g;
665 PyArrayObject *py_g_zero;
666 double cutoff_frequency, temperature;
667 int frequency_point_index;
668
669 Darray *fc3_normal_squared;
670 double *gamma;
671 double *g;
672 char* g_zero;
673 double *frequencies;
674 size_t (*triplets)[3];
675 int *triplet_weights;
676 int num_frequency_points;
677
678 if (!PyArg_ParseTuple(args, "OOOOOdOOdi",
679 &py_gamma,
680 &py_fc3_normal_squared,
681 &py_triplets,
682 &py_triplet_weights,
683 &py_frequencies,
684 &temperature,
685 &py_g,
686 &py_g_zero,
687 &cutoff_frequency,
688 &frequency_point_index)) {
689 return NULL;
690 }
691
692 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
693 gamma = (double*)PyArray_DATA(py_gamma);
694 g = (double*)PyArray_DATA(py_g);
695 g_zero = (char*)PyArray_DATA(py_g_zero);
696 frequencies = (double*)PyArray_DATA(py_frequencies);
697 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
698 triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
699 num_frequency_points = PyArray_DIMS(py_g)[2];
700
701 ph3py_get_imag_self_energy_at_bands_with_g(gamma,
702 fc3_normal_squared,
703 frequencies,
704 triplets,
705 triplet_weights,
706 g,
707 g_zero,
708 temperature,
709 cutoff_frequency,
710 num_frequency_points,
711 frequency_point_index);
712
713 free(fc3_normal_squared);
714 fc3_normal_squared = NULL;
715
716 Py_RETURN_NONE;
717 }
718
719 static PyObject *
py_get_detailed_imag_self_energy_with_g(PyObject * self,PyObject * args)720 py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
721 {
722 PyArrayObject *py_gamma_detail;
723 PyArrayObject *py_gamma_N;
724 PyArrayObject *py_gamma_U;
725 PyArrayObject *py_fc3_normal_squared;
726 PyArrayObject *py_frequencies;
727 PyArrayObject *py_triplets;
728 PyArrayObject *py_triplet_weights;
729 PyArrayObject *py_grid_address;
730 PyArrayObject *py_g;
731 PyArrayObject *py_g_zero;
732 double cutoff_frequency, temperature;
733
734 Darray *fc3_normal_squared;
735 double *gamma_detail;
736 double *gamma_N;
737 double *gamma_U;
738 double *g;
739 char* g_zero;
740 double *frequencies;
741 size_t (*triplets)[3];
742 int *triplet_weights;
743 int *grid_address;
744
745 if (!PyArg_ParseTuple(args, "OOOOOOOOdOOd",
746 &py_gamma_detail,
747 &py_gamma_N,
748 &py_gamma_U,
749 &py_fc3_normal_squared,
750 &py_triplets,
751 &py_triplet_weights,
752 &py_grid_address,
753 &py_frequencies,
754 &temperature,
755 &py_g,
756 &py_g_zero,
757 &cutoff_frequency)) {
758 return NULL;
759 }
760
761 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
762 gamma_detail = (double*)PyArray_DATA(py_gamma_detail);
763 gamma_N = (double*)PyArray_DATA(py_gamma_N);
764 gamma_U = (double*)PyArray_DATA(py_gamma_U);
765 g = (double*)PyArray_DATA(py_g);
766 g_zero = (char*)PyArray_DATA(py_g_zero);
767 frequencies = (double*)PyArray_DATA(py_frequencies);
768 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
769 triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
770 grid_address = (int*)PyArray_DATA(py_grid_address);
771
772 ph3py_get_detailed_imag_self_energy_at_bands_with_g(gamma_detail,
773 gamma_N,
774 gamma_U,
775 fc3_normal_squared,
776 frequencies,
777 triplets,
778 triplet_weights,
779 grid_address,
780 g,
781 g_zero,
782 temperature,
783 cutoff_frequency);
784
785 free(fc3_normal_squared);
786 fc3_normal_squared = NULL;
787
788 Py_RETURN_NONE;
789 }
790
py_get_real_self_energy_at_bands(PyObject * self,PyObject * args)791 static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
792 PyObject *args)
793 {
794 PyArrayObject *py_shift;
795 PyArrayObject *py_fc3_normal_squared;
796 PyArrayObject *py_frequencies;
797 PyArrayObject *py_triplets;
798 PyArrayObject *py_triplet_weights;
799 PyArrayObject *py_band_indices;
800 double epsilon, unit_conversion_factor, cutoff_frequency, temperature;
801
802 Darray *fc3_normal_squared;
803 double *shift;
804 double *frequencies;
805 int *band_indices;
806 size_t (*triplets)[3];
807 int *triplet_weights;
808
809 if (!PyArg_ParseTuple(args, "OOOOOOdddd",
810 &py_shift,
811 &py_fc3_normal_squared,
812 &py_triplets,
813 &py_triplet_weights,
814 &py_frequencies,
815 &py_band_indices,
816 &temperature,
817 &epsilon,
818 &unit_conversion_factor,
819 &cutoff_frequency)) {
820 return NULL;
821 }
822
823
824 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
825 shift = (double*)PyArray_DATA(py_shift);
826 frequencies = (double*)PyArray_DATA(py_frequencies);
827 band_indices = (int*)PyArray_DATA(py_band_indices);
828 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
829 triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
830
831 ph3py_get_real_self_energy_at_bands(shift,
832 fc3_normal_squared,
833 band_indices,
834 frequencies,
835 triplets,
836 triplet_weights,
837 epsilon,
838 temperature,
839 unit_conversion_factor,
840 cutoff_frequency);
841
842 free(fc3_normal_squared);
843 fc3_normal_squared = NULL;
844
845 Py_RETURN_NONE;
846 }
847
py_get_real_self_energy_at_frequency_point(PyObject * self,PyObject * args)848 static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
849 PyObject *args)
850 {
851 PyArrayObject *py_shift;
852 PyArrayObject *py_fc3_normal_squared;
853 PyArrayObject *py_frequencies;
854 PyArrayObject *py_triplets;
855 PyArrayObject *py_triplet_weights;
856 PyArrayObject *py_band_indices;
857 double frequency_point, epsilon, unit_conversion_factor, cutoff_frequency;
858 double temperature;
859
860 Darray *fc3_normal_squared;
861 double *shift;
862 double *frequencies;
863 int *band_indices;
864 size_t (*triplets)[3];
865 int *triplet_weights;
866
867 if (!PyArg_ParseTuple(args, "OdOOOOOdddd",
868 &py_shift,
869 &frequency_point,
870 &py_fc3_normal_squared,
871 &py_triplets,
872 &py_triplet_weights,
873 &py_frequencies,
874 &py_band_indices,
875 &temperature,
876 &epsilon,
877 &unit_conversion_factor,
878 &cutoff_frequency)) {
879 return NULL;
880 }
881
882
883 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
884 shift = (double*)PyArray_DATA(py_shift);
885 frequencies = (double*)PyArray_DATA(py_frequencies);
886 band_indices = (int*)PyArray_DATA(py_band_indices);
887 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
888 triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
889
890 ph3py_get_real_self_energy_at_frequency_point(shift,
891 frequency_point,
892 fc3_normal_squared,
893 band_indices,
894 frequencies,
895 triplets,
896 triplet_weights,
897 epsilon,
898 temperature,
899 unit_conversion_factor,
900 cutoff_frequency);
901
902 free(fc3_normal_squared);
903 fc3_normal_squared = NULL;
904
905 Py_RETURN_NONE;
906 }
907
py_get_collision_matrix(PyObject * self,PyObject * args)908 static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
909 {
910 PyArrayObject *py_collision_matrix;
911 PyArrayObject *py_fc3_normal_squared;
912 PyArrayObject *py_frequencies;
913 PyArrayObject *py_triplets;
914 PyArrayObject *py_triplets_map;
915 PyArrayObject *py_map_q;
916 PyArrayObject *py_g;
917 PyArrayObject *py_rotated_grid_points;
918 PyArrayObject *py_rotations_cartesian;
919 double temperature, unit_conversion_factor, cutoff_frequency;
920
921 Darray *fc3_normal_squared;
922 double *collision_matrix;
923 double *g;
924 double *frequencies;
925 size_t (*triplets)[3];
926 size_t *triplets_map;
927 size_t *map_q;
928 size_t *rotated_grid_points;
929 npy_intp num_gp, num_ir_gp, num_rot;
930 double *rotations_cartesian;
931
932 if (!PyArg_ParseTuple(args, "OOOOOOOOOddd",
933 &py_collision_matrix,
934 &py_fc3_normal_squared,
935 &py_frequencies,
936 &py_g,
937 &py_triplets,
938 &py_triplets_map,
939 &py_map_q,
940 &py_rotated_grid_points,
941 &py_rotations_cartesian,
942 &temperature,
943 &unit_conversion_factor,
944 &cutoff_frequency)) {
945 return NULL;
946 }
947
948 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
949 collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
950 g = (double*)PyArray_DATA(py_g);
951 frequencies = (double*)PyArray_DATA(py_frequencies);
952 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
953 triplets_map = (size_t*)PyArray_DATA(py_triplets_map);
954 num_gp = PyArray_DIMS(py_triplets_map)[0];
955 map_q = (size_t*)PyArray_DATA(py_map_q);
956 rotated_grid_points = (size_t*)PyArray_DATA(py_rotated_grid_points);
957 num_ir_gp = PyArray_DIMS(py_rotated_grid_points)[0];
958 num_rot = PyArray_DIMS(py_rotated_grid_points)[1];
959 rotations_cartesian = (double*)PyArray_DATA(py_rotations_cartesian);
960
961 assert(num_rot == PyArray_DIMS(py_rotations_cartesian)[0]);
962 assert(num_gp == PyArray_DIMS(py_frequencies)[0]);
963
964 ph3py_get_collision_matrix(collision_matrix,
965 fc3_normal_squared,
966 frequencies,
967 triplets,
968 triplets_map,
969 map_q,
970 rotated_grid_points,
971 rotations_cartesian,
972 g,
973 num_ir_gp,
974 num_gp,
975 num_rot,
976 temperature,
977 unit_conversion_factor,
978 cutoff_frequency);
979
980 free(fc3_normal_squared);
981 fc3_normal_squared = NULL;
982
983 Py_RETURN_NONE;
984 }
985
py_get_reducible_collision_matrix(PyObject * self,PyObject * args)986 static PyObject * py_get_reducible_collision_matrix(PyObject *self, PyObject *args)
987 {
988 PyArrayObject *py_collision_matrix;
989 PyArrayObject *py_fc3_normal_squared;
990 PyArrayObject *py_frequencies;
991 PyArrayObject *py_triplets;
992 PyArrayObject *py_triplets_map;
993 PyArrayObject *py_map_q;
994 PyArrayObject *py_g;
995 double temperature, unit_conversion_factor, cutoff_frequency;
996
997 Darray *fc3_normal_squared;
998 double *collision_matrix;
999 double *g;
1000 double *frequencies;
1001 size_t (*triplets)[3];
1002 size_t *triplets_map;
1003 npy_intp num_gp;
1004 size_t *map_q;
1005
1006 if (!PyArg_ParseTuple(args, "OOOOOOOddd",
1007 &py_collision_matrix,
1008 &py_fc3_normal_squared,
1009 &py_frequencies,
1010 &py_g,
1011 &py_triplets,
1012 &py_triplets_map,
1013 &py_map_q,
1014 &temperature,
1015 &unit_conversion_factor,
1016 &cutoff_frequency)) {
1017 return NULL;
1018 }
1019
1020 fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
1021 collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
1022 g = (double*)PyArray_DATA(py_g);
1023 frequencies = (double*)PyArray_DATA(py_frequencies);
1024 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
1025 triplets_map = (size_t*)PyArray_DATA(py_triplets_map);
1026 num_gp = PyArray_DIMS(py_triplets_map)[0];
1027 map_q = (size_t*)PyArray_DATA(py_map_q);
1028
1029 ph3py_get_reducible_collision_matrix(collision_matrix,
1030 fc3_normal_squared,
1031 frequencies,
1032 triplets,
1033 triplets_map,
1034 map_q,
1035 g,
1036 num_gp,
1037 temperature,
1038 unit_conversion_factor,
1039 cutoff_frequency);
1040
1041 free(fc3_normal_squared);
1042 fc3_normal_squared = NULL;
1043
1044 Py_RETURN_NONE;
1045 }
1046
py_symmetrize_collision_matrix(PyObject * self,PyObject * args)1047 static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args)
1048 {
1049 PyArrayObject *py_collision_matrix;
1050
1051 double *collision_matrix;
1052 long num_band, num_grid_points, num_temp, num_sigma;
1053 long num_column;
1054
1055 if (!PyArg_ParseTuple(args, "O",
1056 &py_collision_matrix)) {
1057 return NULL;
1058 }
1059
1060 collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
1061 num_sigma = PyArray_DIMS(py_collision_matrix)[0];
1062 num_temp = PyArray_DIMS(py_collision_matrix)[1];
1063 num_grid_points = PyArray_DIMS(py_collision_matrix)[2];
1064 num_band = PyArray_DIMS(py_collision_matrix)[3];
1065
1066 if (PyArray_NDIM(py_collision_matrix) == 8) {
1067 num_column = num_grid_points * num_band * 3;
1068 } else {
1069 num_column = num_grid_points * num_band;
1070 }
1071
1072 ph3py_symmetrize_collision_matrix(collision_matrix,
1073 num_column,
1074 num_temp,
1075 num_sigma);
1076
1077 Py_RETURN_NONE;
1078 }
1079
py_expand_collision_matrix(PyObject * self,PyObject * args)1080 static PyObject * py_expand_collision_matrix(PyObject *self, PyObject *args)
1081 {
1082 PyArrayObject *py_collision_matrix;
1083 PyArrayObject *py_ir_grid_points;
1084 PyArrayObject *py_rot_grid_points;
1085
1086 double *collision_matrix;
1087 size_t *rot_grid_points;
1088 size_t *ir_grid_points;
1089 long num_band, num_grid_points, num_temp, num_sigma, num_rot, num_ir_gp;
1090
1091 if (!PyArg_ParseTuple(args, "OOO",
1092 &py_collision_matrix,
1093 &py_ir_grid_points,
1094 &py_rot_grid_points)) {
1095 return NULL;
1096 }
1097
1098 collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
1099 rot_grid_points = (size_t*)PyArray_DATA(py_rot_grid_points);
1100 ir_grid_points = (size_t*)PyArray_DATA(py_ir_grid_points);
1101 num_sigma = (long)PyArray_DIMS(py_collision_matrix)[0];
1102 num_temp = (long)PyArray_DIMS(py_collision_matrix)[1];
1103 num_grid_points = (long)PyArray_DIMS(py_collision_matrix)[2];
1104 num_band = (long)PyArray_DIMS(py_collision_matrix)[3];
1105 num_rot = (long)PyArray_DIMS(py_rot_grid_points)[0];
1106 num_ir_gp = (long)PyArray_DIMS(py_ir_grid_points)[0];
1107
1108 ph3py_expand_collision_matrix(collision_matrix,
1109 rot_grid_points,
1110 ir_grid_points,
1111 num_ir_gp,
1112 num_grid_points,
1113 num_rot,
1114 num_sigma,
1115 num_temp,
1116 num_band);
1117
1118 Py_RETURN_NONE;
1119 }
1120
py_get_isotope_strength(PyObject * self,PyObject * args)1121 static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
1122 {
1123 PyArrayObject *py_gamma;
1124 PyArrayObject *py_frequencies;
1125 PyArrayObject *py_eigenvectors;
1126 PyArrayObject *py_band_indices;
1127 PyArrayObject *py_mass_variances;
1128 long grid_point;
1129 long num_grid_points;
1130 double cutoff_frequency;
1131 double sigma;
1132
1133 double *gamma;
1134 double *frequencies;
1135 lapack_complex_double *eigenvectors;
1136 int *band_indices;
1137 double *mass_variances;
1138 npy_intp num_band, num_band0;
1139
1140 if (!PyArg_ParseTuple(args, "OlOOOOldd",
1141 &py_gamma,
1142 &grid_point,
1143 &py_mass_variances,
1144 &py_frequencies,
1145 &py_eigenvectors,
1146 &py_band_indices,
1147 &num_grid_points,
1148 &sigma,
1149 &cutoff_frequency)) {
1150 return NULL;
1151 }
1152
1153
1154 gamma = (double*)PyArray_DATA(py_gamma);
1155 frequencies = (double*)PyArray_DATA(py_frequencies);
1156 eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
1157 band_indices = (int*)PyArray_DATA(py_band_indices);
1158 mass_variances = (double*)PyArray_DATA(py_mass_variances);
1159 num_band = PyArray_DIMS(py_frequencies)[1];
1160 num_band0 = PyArray_DIMS(py_band_indices)[0];
1161
1162 ph3py_get_isotope_scattering_strength(gamma,
1163 grid_point,
1164 mass_variances,
1165 frequencies,
1166 eigenvectors,
1167 num_grid_points,
1168 band_indices,
1169 num_band,
1170 num_band0,
1171 sigma,
1172 cutoff_frequency);
1173
1174 Py_RETURN_NONE;
1175 }
1176
py_get_thm_isotope_strength(PyObject * self,PyObject * args)1177 static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
1178 {
1179 PyArrayObject *py_gamma;
1180 PyArrayObject *py_frequencies;
1181 PyArrayObject *py_eigenvectors;
1182 PyArrayObject *py_band_indices;
1183 PyArrayObject *py_mass_variances;
1184 PyArrayObject *py_ir_grid_points;
1185 PyArrayObject *py_weights;
1186 PyArrayObject *py_integration_weights;
1187 long grid_point;
1188 double cutoff_frequency;
1189
1190 double *gamma;
1191 double *frequencies;
1192 size_t *ir_grid_points;
1193 int *weights;
1194 lapack_complex_double *eigenvectors;
1195 int *band_indices;
1196 double *mass_variances;
1197 npy_intp num_band, num_band0, num_ir_grid_points;
1198 double *integration_weights;
1199
1200 if (!PyArg_ParseTuple(args, "OlOOOOOOOd",
1201 &py_gamma,
1202 &grid_point,
1203 &py_ir_grid_points,
1204 &py_weights,
1205 &py_mass_variances,
1206 &py_frequencies,
1207 &py_eigenvectors,
1208 &py_band_indices,
1209 &py_integration_weights,
1210 &cutoff_frequency)) {
1211 return NULL;
1212 }
1213
1214
1215 gamma = (double*)PyArray_DATA(py_gamma);
1216 frequencies = (double*)PyArray_DATA(py_frequencies);
1217 ir_grid_points = (size_t*)PyArray_DATA(py_ir_grid_points);
1218 weights = (int*)PyArray_DATA(py_weights);
1219 eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
1220 band_indices = (int*)PyArray_DATA(py_band_indices);
1221 mass_variances = (double*)PyArray_DATA(py_mass_variances);
1222 num_band = PyArray_DIMS(py_frequencies)[1];
1223 num_band0 = PyArray_DIMS(py_band_indices)[0];
1224 integration_weights = (double*)PyArray_DATA(py_integration_weights);
1225 num_ir_grid_points = PyArray_DIMS(py_ir_grid_points)[0];
1226
1227 ph3py_get_thm_isotope_scattering_strength(gamma,
1228 grid_point,
1229 ir_grid_points,
1230 weights,
1231 mass_variances,
1232 frequencies,
1233 eigenvectors,
1234 num_ir_grid_points,
1235 band_indices,
1236 num_band,
1237 num_band0,
1238 integration_weights,
1239 cutoff_frequency);
1240
1241 Py_RETURN_NONE;
1242 }
1243
py_distribute_fc3(PyObject * self,PyObject * args)1244 static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
1245 {
1246 PyArrayObject *force_constants_third;
1247 int target;
1248 int source;
1249 PyArrayObject *rotation_cart_inv;
1250 PyArrayObject *atom_mapping_py;
1251
1252 double *fc3;
1253 double *rot_cart_inv;
1254 int *atom_mapping;
1255 npy_intp num_atom;
1256
1257 if (!PyArg_ParseTuple(args, "OiiOO",
1258 &force_constants_third,
1259 &target,
1260 &source,
1261 &atom_mapping_py,
1262 &rotation_cart_inv)) {
1263 return NULL;
1264 }
1265
1266 fc3 = (double*)PyArray_DATA(force_constants_third);
1267 rot_cart_inv = (double*)PyArray_DATA(rotation_cart_inv);
1268 atom_mapping = (int*)PyArray_DATA(atom_mapping_py);
1269 num_atom = PyArray_DIMS(atom_mapping_py)[0];
1270
1271 ph3py_distribute_fc3(fc3,
1272 target,
1273 source,
1274 atom_mapping,
1275 num_atom,
1276 rot_cart_inv);
1277
1278 Py_RETURN_NONE;
1279 }
1280
py_rotate_delta_fc2s(PyObject * self,PyObject * args)1281 static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
1282 {
1283 PyArrayObject *py_fc3;
1284 PyArrayObject *py_delta_fc2s;
1285 PyArrayObject *py_inv_U;
1286 PyArrayObject *py_site_sym_cart;
1287 PyArrayObject *py_rot_map_syms;
1288
1289 double (*fc3)[3][3][3];
1290 double (*delta_fc2s)[3][3];
1291 double *inv_U;
1292 double (*site_sym_cart)[3][3];
1293 int *rot_map_syms;
1294 npy_intp num_atom, num_disp, num_site_sym;
1295
1296 if (!PyArg_ParseTuple(args, "OOOOO",
1297 &py_fc3,
1298 &py_delta_fc2s,
1299 &py_inv_U,
1300 &py_site_sym_cart,
1301 &py_rot_map_syms)) {
1302 return NULL;
1303 }
1304
1305 /* (num_atom, num_atom, 3, 3, 3) */
1306 fc3 = (double(*)[3][3][3])PyArray_DATA(py_fc3);
1307 /* (n_u1, num_atom, num_atom, 3, 3) */
1308 delta_fc2s = (double(*)[3][3])PyArray_DATA(py_delta_fc2s);
1309 /* (3, n_u1 * n_sym) */
1310 inv_U = (double*)PyArray_DATA(py_inv_U);
1311 /* (n_sym, 3, 3) */
1312 site_sym_cart = (double(*)[3][3])PyArray_DATA(py_site_sym_cart);
1313 /* (n_sym, natom) */
1314 rot_map_syms = (int*)PyArray_DATA(py_rot_map_syms);
1315
1316 num_atom = PyArray_DIMS(py_fc3)[0];
1317 num_disp = PyArray_DIMS(py_delta_fc2s)[0];
1318 num_site_sym = PyArray_DIMS(py_site_sym_cart)[0];
1319
1320 ph3py_rotate_delta_fc2(fc3,
1321 delta_fc2s,
1322 inv_U,
1323 site_sym_cart,
1324 rot_map_syms,
1325 num_atom,
1326 num_site_sym,
1327 num_disp);
1328
1329 Py_RETURN_NONE;
1330 }
1331
1332 static PyObject *
py_set_permutation_symmetry_fc3(PyObject * self,PyObject * args)1333 py_set_permutation_symmetry_fc3(PyObject *self, PyObject *args)
1334 {
1335 PyArrayObject *py_fc3;
1336
1337 double *fc3;
1338 npy_intp num_atom;
1339
1340 if (!PyArg_ParseTuple(args, "O", &py_fc3)) {
1341 return NULL;
1342 }
1343
1344 fc3 = (double*)PyArray_DATA(py_fc3);
1345 num_atom = PyArray_DIMS(py_fc3)[0];
1346
1347 ph3py_set_permutation_symmetry_fc3(fc3, num_atom);
1348
1349 Py_RETURN_NONE;
1350 }
1351
1352 static PyObject *
py_set_permutation_symmetry_compact_fc3(PyObject * self,PyObject * args)1353 py_set_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
1354 {
1355 PyArrayObject* py_fc3;
1356 PyArrayObject* py_permutations;
1357 PyArrayObject* py_s2pp_map;
1358 PyArrayObject* py_p2s_map;
1359 PyArrayObject* py_nsym_list;
1360
1361 double *fc3;
1362 int *s2pp;
1363 int *p2s;
1364 int *nsym_list;
1365 int *perms;
1366 npy_intp n_patom, n_satom;
1367
1368 if (!PyArg_ParseTuple(args, "OOOOO",
1369 &py_fc3,
1370 &py_permutations,
1371 &py_s2pp_map,
1372 &py_p2s_map,
1373 &py_nsym_list)) {
1374 return NULL;
1375 }
1376
1377 fc3 = (double*)PyArray_DATA(py_fc3);
1378 perms = (int*)PyArray_DATA(py_permutations);
1379 s2pp = (int*)PyArray_DATA(py_s2pp_map);
1380 p2s = (int*)PyArray_DATA(py_p2s_map);
1381 nsym_list = (int*)PyArray_DATA(py_nsym_list);
1382 n_patom = PyArray_DIMS(py_fc3)[0];
1383 n_satom = PyArray_DIMS(py_fc3)[1];
1384
1385 ph3py_set_permutation_symmetry_compact_fc3(fc3,
1386 p2s,
1387 s2pp,
1388 nsym_list,
1389 perms,
1390 n_satom,
1391 n_patom);
1392
1393 Py_RETURN_NONE;
1394 }
1395
py_transpose_compact_fc3(PyObject * self,PyObject * args)1396 static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
1397 {
1398 PyArrayObject* py_fc3;
1399 PyArrayObject* py_permutations;
1400 PyArrayObject* py_s2pp_map;
1401 PyArrayObject* py_p2s_map;
1402 PyArrayObject* py_nsym_list;
1403 int t_type;
1404
1405 double *fc3;
1406 int *s2pp;
1407 int *p2s;
1408 int *nsym_list;
1409 int *perms;
1410 npy_intp n_patom, n_satom;
1411
1412 if (!PyArg_ParseTuple(args, "OOOOOi",
1413 &py_fc3,
1414 &py_permutations,
1415 &py_s2pp_map,
1416 &py_p2s_map,
1417 &py_nsym_list,
1418 &t_type)) {
1419 return NULL;
1420 }
1421
1422 fc3 = (double*)PyArray_DATA(py_fc3);
1423 perms = (int*)PyArray_DATA(py_permutations);
1424 s2pp = (int*)PyArray_DATA(py_s2pp_map);
1425 p2s = (int*)PyArray_DATA(py_p2s_map);
1426 nsym_list = (int*)PyArray_DATA(py_nsym_list);
1427 n_patom = PyArray_DIMS(py_fc3)[0];
1428 n_satom = PyArray_DIMS(py_fc3)[1];
1429
1430 ph3py_transpose_compact_fc3(fc3,
1431 p2s,
1432 s2pp,
1433 nsym_list,
1434 perms,
1435 n_satom,
1436 n_patom,
1437 t_type);
1438
1439 Py_RETURN_NONE;
1440 }
1441
py_get_neighboring_gird_points(PyObject * self,PyObject * args)1442 static PyObject * py_get_neighboring_gird_points(PyObject *self, PyObject *args)
1443 {
1444 PyArrayObject *py_relative_grid_points;
1445 PyArrayObject *py_grid_points;
1446 PyArrayObject *py_relative_grid_address;
1447 PyArrayObject *py_mesh;
1448 PyArrayObject *py_bz_grid_address;
1449 PyArrayObject *py_bz_map;
1450
1451 size_t *relative_grid_points;
1452 size_t *grid_points;
1453 npy_intp num_grid_points, num_relative_grid_address;
1454 int (*relative_grid_address)[3];
1455 int *mesh;
1456 int (*bz_grid_address)[3];
1457 size_t *bz_map;
1458
1459 if (!PyArg_ParseTuple(args, "OOOOOO",
1460 &py_relative_grid_points,
1461 &py_grid_points,
1462 &py_relative_grid_address,
1463 &py_mesh,
1464 &py_bz_grid_address,
1465 &py_bz_map)) {
1466 return NULL;
1467 }
1468
1469 relative_grid_points = (size_t*)PyArray_DATA(py_relative_grid_points);
1470 grid_points = (size_t*)PyArray_DATA(py_grid_points);
1471 num_grid_points = PyArray_DIMS(py_grid_points)[0];
1472 relative_grid_address = (int(*)[3])PyArray_DATA(py_relative_grid_address);
1473 num_relative_grid_address = PyArray_DIMS(py_relative_grid_address)[0];
1474 mesh = (int*)PyArray_DATA(py_mesh);
1475 bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
1476 bz_map = (size_t*)PyArray_DATA(py_bz_map);
1477
1478 ph3py_get_neighboring_gird_points(relative_grid_points,
1479 grid_points,
1480 relative_grid_address,
1481 mesh,
1482 bz_grid_address,
1483 bz_map,
1484 num_grid_points,
1485 num_relative_grid_address);
1486
1487 Py_RETURN_NONE;
1488 }
1489
py_set_integration_weights(PyObject * self,PyObject * args)1490 static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
1491 {
1492 PyArrayObject *py_iw;
1493 PyArrayObject *py_frequency_points;
1494 PyArrayObject *py_relative_grid_address;
1495 PyArrayObject *py_mesh;
1496 PyArrayObject *py_grid_points;
1497 PyArrayObject *py_frequencies;
1498 PyArrayObject *py_bz_grid_address;
1499 PyArrayObject *py_bz_map;
1500
1501 double *iw;
1502 double *frequency_points;
1503 npy_intp num_band0, num_band, num_gp;
1504 int (*relative_grid_address)[4][3];
1505 int *mesh;
1506 size_t *grid_points;
1507 int (*bz_grid_address)[3];
1508 size_t *bz_map;
1509 double *frequencies;
1510
1511 if (!PyArg_ParseTuple(args, "OOOOOOOO",
1512 &py_iw,
1513 &py_frequency_points,
1514 &py_relative_grid_address,
1515 &py_mesh,
1516 &py_grid_points,
1517 &py_frequencies,
1518 &py_bz_grid_address,
1519 &py_bz_map)) {
1520 return NULL;
1521 }
1522
1523 iw = (double*)PyArray_DATA(py_iw);
1524 frequency_points = (double*)PyArray_DATA(py_frequency_points);
1525 num_band0 = PyArray_DIMS(py_frequency_points)[0];
1526 relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
1527 mesh = (int*)PyArray_DATA(py_mesh);
1528 grid_points = (size_t*)PyArray_DATA(py_grid_points);
1529 num_gp = PyArray_DIMS(py_grid_points)[0];
1530 bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
1531 bz_map = (size_t*)PyArray_DATA(py_bz_map);
1532 frequencies = (double*)PyArray_DATA(py_frequencies);
1533 num_band = PyArray_DIMS(py_frequencies)[1];
1534
1535 ph3py_set_integration_weights(iw,
1536 frequency_points,
1537 num_band0,
1538 num_band,
1539 num_gp,
1540 relative_grid_address,
1541 mesh,
1542 grid_points,
1543 bz_grid_address,
1544 bz_map,
1545 frequencies);
1546
1547 Py_RETURN_NONE;
1548 }
1549
1550 static PyObject *
py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject * self,PyObject * args)1551 py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
1552 {
1553 PyArrayObject *py_map_triplets;
1554 PyArrayObject *py_grid_address;
1555 PyArrayObject *py_map_q;
1556 PyArrayObject *py_mesh;
1557 PyArrayObject *py_rotations;
1558 long fixed_grid_number;
1559 int is_time_reversal;
1560 int swappable;
1561
1562 int (*grid_address)[3];
1563 size_t *map_triplets;
1564 size_t *map_q;
1565 int *mesh;
1566 int (*rot)[3][3];
1567 npy_intp num_rot;
1568 size_t num_ir;
1569
1570 if (!PyArg_ParseTuple(args, "OOOlOiOi",
1571 &py_map_triplets,
1572 &py_map_q,
1573 &py_grid_address,
1574 &fixed_grid_number,
1575 &py_mesh,
1576 &is_time_reversal,
1577 &py_rotations,
1578 &swappable)) {
1579 return NULL;
1580 }
1581
1582 grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
1583 map_triplets = (size_t*)PyArray_DATA(py_map_triplets);
1584 map_q = (size_t*)PyArray_DATA(py_map_q);
1585 mesh = (int*)PyArray_DATA(py_mesh);
1586 rot = (int(*)[3][3])PyArray_DATA(py_rotations);
1587 num_rot = PyArray_DIMS(py_rotations)[0];
1588 num_ir = ph3py_get_triplets_reciprocal_mesh_at_q(map_triplets,
1589 map_q,
1590 grid_address,
1591 fixed_grid_number,
1592 mesh,
1593 is_time_reversal,
1594 num_rot,
1595 rot,
1596 swappable);
1597
1598 return PyLong_FromSize_t(num_ir);
1599 }
1600
py_tpl_get_BZ_triplets_at_q(PyObject * self,PyObject * args)1601 static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
1602 {
1603 PyArrayObject *py_triplets;
1604 PyArrayObject *py_bz_grid_address;
1605 PyArrayObject *py_bz_map;
1606 PyArrayObject *py_map_triplets;
1607 PyArrayObject *py_mesh;
1608 long grid_point;
1609
1610 size_t (*triplets)[3];
1611 int (*bz_grid_address)[3];
1612 size_t *bz_map;
1613 size_t *map_triplets;
1614 npy_intp num_map_triplets;
1615 int *mesh;
1616 size_t num_ir;
1617
1618 if (!PyArg_ParseTuple(args, "OlOOOO",
1619 &py_triplets,
1620 &grid_point,
1621 &py_bz_grid_address,
1622 &py_bz_map,
1623 &py_map_triplets,
1624 &py_mesh)) {
1625 return NULL;
1626 }
1627
1628 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
1629 bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
1630 bz_map = (size_t*)PyArray_DATA(py_bz_map);
1631 map_triplets = (size_t*)PyArray_DATA(py_map_triplets);
1632 num_map_triplets = PyArray_DIMS(py_map_triplets)[0];
1633 mesh = (int*)PyArray_DATA(py_mesh);
1634
1635 num_ir = ph3py_get_BZ_triplets_at_q(triplets,
1636 grid_point,
1637 bz_grid_address,
1638 bz_map,
1639 map_triplets,
1640 num_map_triplets,
1641 mesh);
1642
1643 return PyLong_FromSize_t(num_ir);
1644 }
1645
1646 static PyObject *
py_set_triplets_integration_weights(PyObject * self,PyObject * args)1647 py_set_triplets_integration_weights(PyObject *self, PyObject *args)
1648 {
1649 PyArrayObject *py_iw;
1650 PyArrayObject *py_iw_zero;
1651 PyArrayObject *py_frequency_points;
1652 PyArrayObject *py_relative_grid_address;
1653 PyArrayObject *py_mesh;
1654 PyArrayObject *py_triplets;
1655 PyArrayObject *py_frequencies1;
1656 PyArrayObject *py_frequencies2;
1657 PyArrayObject *py_bz_grid_address;
1658 PyArrayObject *py_bz_map;
1659 int tp_type;
1660
1661 double *iw;
1662 char *iw_zero;
1663 double *frequency_points;
1664 int (*relative_grid_address)[4][3];
1665 int *mesh;
1666 size_t (*triplets)[3];
1667 int (*bz_grid_address)[3];
1668 size_t *bz_map;
1669 double *frequencies1, *frequencies2;
1670 npy_intp num_band0, num_band1, num_band2, num_triplets;
1671
1672 if (!PyArg_ParseTuple(args, "OOOOOOOOOOi",
1673 &py_iw,
1674 &py_iw_zero,
1675 &py_frequency_points,
1676 &py_relative_grid_address,
1677 &py_mesh,
1678 &py_triplets,
1679 &py_frequencies1,
1680 &py_frequencies2,
1681 &py_bz_grid_address,
1682 &py_bz_map,
1683 &tp_type)) {
1684 return NULL;
1685 }
1686
1687 iw = (double*)PyArray_DATA(py_iw);
1688 iw_zero = (char*)PyArray_DATA(py_iw_zero);
1689 frequency_points = (double*)PyArray_DATA(py_frequency_points);
1690 num_band0 = PyArray_DIMS(py_frequency_points)[0];
1691 relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
1692 mesh = (int*)PyArray_DATA(py_mesh);
1693 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
1694 num_triplets = PyArray_DIMS(py_triplets)[0];
1695 bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
1696 bz_map = (size_t*)PyArray_DATA(py_bz_map);
1697 frequencies1 = (double*)PyArray_DATA(py_frequencies1);
1698 frequencies2 = (double*)PyArray_DATA(py_frequencies2);
1699 num_band1 = PyArray_DIMS(py_frequencies1)[1];
1700 num_band2 = PyArray_DIMS(py_frequencies2)[1];
1701
1702 ph3py_get_integration_weight(iw,
1703 iw_zero,
1704 frequency_points,
1705 num_band0,
1706 relative_grid_address,
1707 mesh,
1708 triplets,
1709 num_triplets,
1710 bz_grid_address,
1711 bz_map,
1712 frequencies1,
1713 num_band1,
1714 frequencies2,
1715 num_band2,
1716 tp_type,
1717 1,
1718 0);
1719
1720 Py_RETURN_NONE;
1721 }
1722
1723 static PyObject *
py_set_triplets_integration_weights_with_sigma(PyObject * self,PyObject * args)1724 py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
1725 {
1726 PyArrayObject *py_iw;
1727 PyArrayObject *py_iw_zero;
1728 PyArrayObject *py_frequency_points;
1729 PyArrayObject *py_triplets;
1730 PyArrayObject *py_frequencies;
1731 double sigma, sigma_cutoff;
1732
1733 double *iw;
1734 char *iw_zero;
1735 double *frequency_points;
1736 size_t (*triplets)[3];
1737 double *frequencies;
1738 npy_intp num_band0, num_band, num_iw, num_triplets;
1739
1740 if (!PyArg_ParseTuple(args, "OOOOOdd",
1741 &py_iw,
1742 &py_iw_zero,
1743 &py_frequency_points,
1744 &py_triplets,
1745 &py_frequencies,
1746 &sigma,
1747 &sigma_cutoff)) {
1748 return NULL;
1749 }
1750
1751 iw = (double*)PyArray_DATA(py_iw);
1752 iw_zero = (char*)PyArray_DATA(py_iw_zero);
1753 frequency_points = (double*)PyArray_DATA(py_frequency_points);
1754 num_band0 = PyArray_DIMS(py_frequency_points)[0];
1755 triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
1756 num_triplets = PyArray_DIMS(py_triplets)[0];
1757 frequencies = (double*)PyArray_DATA(py_frequencies);
1758 num_band = PyArray_DIMS(py_frequencies)[1];
1759 num_iw = PyArray_DIMS(py_iw)[0];
1760
1761 ph3py_get_integration_weight_with_sigma(iw,
1762 iw_zero,
1763 sigma,
1764 sigma_cutoff,
1765 frequency_points,
1766 num_band0,
1767 triplets,
1768 num_triplets,
1769 frequencies,
1770 num_band,
1771 num_iw);
1772
1773 Py_RETURN_NONE;
1774 }
1775
1776 static PyObject *
py_diagonalize_collision_matrix(PyObject * self,PyObject * args)1777 py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
1778 {
1779 PyArrayObject *py_collision_matrix;
1780 PyArrayObject *py_eigenvalues;
1781 double cutoff;
1782 int i_sigma, i_temp, is_pinv, solver;
1783
1784 double *collision_matrix;
1785 double *eigvals;
1786 npy_intp num_temp, num_grid_point, num_band;
1787 size_t num_column, adrs_shift;
1788 int info;
1789
1790 if (!PyArg_ParseTuple(args, "OOiidii",
1791 &py_collision_matrix,
1792 &py_eigenvalues,
1793 &i_sigma,
1794 &i_temp,
1795 &cutoff,
1796 &solver,
1797 &is_pinv)) {
1798 return NULL;
1799 }
1800
1801 collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
1802 eigvals = (double*)PyArray_DATA(py_eigenvalues);
1803
1804 if (PyArray_NDIM(py_collision_matrix) == 2) {
1805 num_temp = 1;
1806 num_column = PyArray_DIM(py_collision_matrix, 1);
1807 } else {
1808 num_temp = PyArray_DIM(py_collision_matrix, 1);
1809 num_grid_point = PyArray_DIM(py_collision_matrix, 2);
1810 num_band = PyArray_DIM(py_collision_matrix, 3);
1811 if (PyArray_NDIM(py_collision_matrix) == 8) {
1812 num_column = num_grid_point * num_band * 3;
1813 } else {
1814 num_column = num_grid_point * num_band;
1815 }
1816 }
1817 adrs_shift = (i_sigma * num_column * num_column * num_temp +
1818 i_temp * num_column * num_column);
1819
1820 /* show_colmat_info(py_collision_matrix, i_sigma, i_temp, adrs_shift); */
1821
1822 info = phonopy_dsyev(collision_matrix + adrs_shift,
1823 eigvals, num_column, solver);
1824 if (is_pinv) {
1825 pinv_from_eigensolution(collision_matrix + adrs_shift,
1826 eigvals, num_column, cutoff, 0);
1827 }
1828
1829 return PyLong_FromLong((long) info);
1830 }
1831
py_pinv_from_eigensolution(PyObject * self,PyObject * args)1832 static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
1833 {
1834 PyArrayObject *py_collision_matrix;
1835 PyArrayObject *py_eigenvalues;
1836 double cutoff;
1837 int i_sigma, i_temp, pinv_method;
1838
1839 double *collision_matrix;
1840 double *eigvals;
1841 npy_intp num_temp, num_grid_point, num_band;
1842 size_t num_column, adrs_shift;
1843
1844 if (!PyArg_ParseTuple(args, "OOiidi",
1845 &py_collision_matrix,
1846 &py_eigenvalues,
1847 &i_sigma,
1848 &i_temp,
1849 &cutoff,
1850 &pinv_method)) {
1851 return NULL;
1852 }
1853
1854 collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
1855 eigvals = (double*)PyArray_DATA(py_eigenvalues);
1856 num_temp = PyArray_DIMS(py_collision_matrix)[1];
1857 num_grid_point = PyArray_DIMS(py_collision_matrix)[2];
1858 num_band = PyArray_DIMS(py_collision_matrix)[3];
1859
1860 if (PyArray_NDIM(py_collision_matrix) == 8) {
1861 num_column = num_grid_point * num_band * 3;
1862 } else {
1863 num_column = num_grid_point * num_band;
1864 }
1865 adrs_shift = (i_sigma * num_column * num_column * num_temp +
1866 i_temp * num_column * num_column);
1867
1868 /* show_colmat_info(py_collision_matrix, i_sigma, i_temp, adrs_shift); */
1869
1870 pinv_from_eigensolution(collision_matrix + adrs_shift,
1871 eigvals, num_column, cutoff, pinv_method);
1872
1873 Py_RETURN_NONE;
1874 }
1875
py_get_default_colmat_solver(PyObject * self,PyObject * args)1876 static PyObject * py_get_default_colmat_solver(PyObject *self, PyObject *args)
1877 {
1878 if (!PyArg_ParseTuple(args, "")) {
1879 return NULL;
1880 }
1881
1882 #ifdef MKL_LAPACKE
1883 return PyLong_FromLong((long) 1);
1884 #else
1885 return PyLong_FromLong((long) 4);
1886 #endif
1887
1888 }
1889
pinv_from_eigensolution(double * data,const double * eigvals,const size_t size,const double cutoff,const int pinv_method)1890 static void pinv_from_eigensolution(double *data,
1891 const double *eigvals,
1892 const size_t size,
1893 const double cutoff,
1894 const int pinv_method)
1895 {
1896 size_t i, ib, j, k, max_l, i_s, j_s;
1897 double *tmp_data;
1898 double e, sum;
1899 size_t *l;
1900
1901 l = NULL;
1902 tmp_data = NULL;
1903
1904 tmp_data = (double*)malloc(sizeof(double) * size * size);
1905
1906 #pragma omp parallel for
1907 for (i = 0; i < size * size; i++) {
1908 tmp_data[i] = data[i];
1909 }
1910
1911 l = (size_t*)malloc(sizeof(size_t) * size);
1912 max_l = 0;
1913 for (i = 0; i < size; i++) {
1914 if (pinv_method == 0) {
1915 e = fabs(eigvals[i]);
1916 } else {
1917 e = eigvals[i];
1918 }
1919 if (e > cutoff) {
1920 l[max_l] = i;
1921 max_l++;
1922 }
1923 }
1924
1925 #pragma omp parallel for private(ib, j, k, i_s, j_s, sum)
1926 for (i = 0; i < size / 2; i++) {
1927 /* from front */
1928 i_s = i * size;
1929 for (j = i; j < size; j++) {
1930 j_s = j * size;
1931 sum = 0;
1932 for (k = 0; k < max_l; k++) {
1933 sum += tmp_data[i_s + l[k]] * tmp_data[j_s + l[k]] / eigvals[l[k]];
1934 }
1935 data[i_s + j] = sum;
1936 data[j_s + i] = sum;
1937 }
1938 /* from back */
1939 ib = size - i - 1;
1940 i_s = ib * size;
1941 for (j = ib; j < size; j++) {
1942 j_s = j * size;
1943 sum = 0;
1944 for (k = 0; k < max_l; k++) {
1945 sum += tmp_data[i_s + l[k]] * tmp_data[j_s + l[k]] / eigvals[l[k]];
1946 }
1947 data[i_s + j] = sum;
1948 data[j_s + ib] = sum;
1949 }
1950 }
1951
1952 /* when size is odd */
1953 if ((size % 2) == 1) {
1954 i = (size - 1) / 2;
1955 i_s = i * size;
1956 for (j = i; j < size; j++) {
1957 j_s = j * size;
1958 sum = 0;
1959 for (k = 0; k < max_l; k++) {
1960 sum += tmp_data[i_s + l[k]] * tmp_data[j_s + l[k]] / eigvals[l[k]];
1961 }
1962 data[i_s + j] = sum;
1963 data[j_s + i] = sum;
1964 }
1965 }
1966
1967 free(l);
1968 l = NULL;
1969
1970 free(tmp_data);
1971 tmp_data = NULL;
1972 }
1973
show_colmat_info(const PyArrayObject * py_collision_matrix,const size_t i_sigma,const size_t i_temp,const size_t adrs_shift)1974 static void show_colmat_info(const PyArrayObject *py_collision_matrix,
1975 const size_t i_sigma,
1976 const size_t i_temp,
1977 const size_t adrs_shift)
1978 {
1979 int i;
1980
1981 printf(" Array_shape:(");
1982 for (i = 0; i < PyArray_NDIM(py_collision_matrix); i++) {
1983 printf("%d", (int)PyArray_DIM(py_collision_matrix, i));
1984 if (i < PyArray_NDIM(py_collision_matrix) - 1) {
1985 printf(",");
1986 } else {
1987 printf("), ");
1988 }
1989 }
1990 printf("Data shift:%lu [%lu, %lu]\n", adrs_shift, i_sigma, i_temp);
1991 }
1992
1993
convert_to_iarray(const PyArrayObject * npyary)1994 static Iarray* convert_to_iarray(const PyArrayObject* npyary)
1995 {
1996 int i;
1997 Iarray *ary;
1998
1999 ary = (Iarray*) malloc(sizeof(Iarray));
2000 for (i = 0; i < PyArray_NDIM(npyary); i++) {
2001 ary->dims[i] = PyArray_DIMS(npyary)[i];
2002 }
2003 ary->data = (int*)PyArray_DATA(npyary);
2004 return ary;
2005 }
2006
2007
convert_to_darray(const PyArrayObject * npyary)2008 static Darray* convert_to_darray(const PyArrayObject* npyary)
2009 {
2010 int i;
2011 Darray *ary;
2012
2013 ary = (Darray*) malloc(sizeof(Darray));
2014 for (i = 0; i < PyArray_NDIM(npyary); i++) {
2015 ary->dims[i] = PyArray_DIMS(npyary)[i];
2016 }
2017 ary->data = (double*)PyArray_DATA(npyary);
2018 return ary;
2019 }
2020