Modification toward sharing dynmat.c between _phonopy.c and _phono3py.c

This commit is contained in:
Atsushi Togo 2013-03-06 08:01:22 +09:00
parent f68c189905
commit 890e307426
3 changed files with 76 additions and 104 deletions

View File

@ -4,6 +4,7 @@
#include <math.h>
#include <numpy/arrayobject.h>
#include <lapacke.h>
#include "dynmat.h"
/* Boltzmann constant eV/K */
#define KB 8.6173382568083159E-05
@ -14,18 +15,21 @@ typedef struct {
int d1;
int d2;
int **data;
int *_data;
} Array2D;
typedef struct {
int d1;
int d2;
double **data;
double *_data;
} DArray2D;
typedef struct {
int d1;
int d2;
npy_cdouble **data;
npy_cdouble *_data;
} CArray2D;
typedef struct {
@ -58,8 +62,8 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args);
static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args);
static int get_interaction_strength(double *amps,
const double *q0,
const double *q1,
const double *q2,
const double *q1s,
const double *q2s,
const Array1D *weights,
const double *fc2,
const double *fc3,
@ -214,8 +218,8 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args)
PyArrayObject* amplitude;
PyArrayObject* qvec0;
PyArrayObject* qvec1;
PyArrayObject* qvec2;
PyArrayObject* qvec1s;
PyArrayObject* qvec2s;
PyArrayObject* weights;
PyArrayObject* shortest_vectors;
PyArrayObject* multiplicity;
@ -230,8 +234,8 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args)
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOdiid",
&amplitude,
&qvec0,
&qvec1,
&qvec2,
&qvec1s,
&qvec2s,
&weights,
&shortest_vectors,
&multiplicity,
@ -249,8 +253,8 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args)
double* amps = (double*)amplitude->data;
const double *q0 = (double*)qvec0->data;
const double *q1 = (double*)qvec1->data;
const double *q2 = (double*)qvec2->data;
const double *q1s = (double*)qvec1s->data;
const double *q2s = (double*)qvec2s->data;
const long *weights_long = (long*)weights->data;
const double* masses = (double*)atomic_masses->data;
const long* multiplicity_long = (long*)multiplicity->data;
@ -287,8 +291,8 @@ static PyObject * py_get_interaction_strength(PyObject *self, PyObject *args)
get_interaction_strength(amps,
q0,
q1,
q2,
q1s,
q2s,
w,
fc2,
fc3,
@ -888,8 +892,8 @@ static PyObject * py_phonopy_zheev(PyObject *self, PyObject *args)
static int get_interaction_strength(double *amps,
const double *q0,
const double *q1,
const double *q2,
const double *q1s,
const double *q2s,
const Array1D *weights,
const double *fc2,
const double *fc3,
@ -903,6 +907,17 @@ static int get_interaction_strength(double *amps,
const int r2q_TI_index,
const double symprec)
{
int i, num_triplets;
double q1[3], q2[3];
num_triplets = weights->d1;
for (i = 0; i < num_triplets; i++) {
for (j = 0; j < 3; j++) {
q1[3] = q1s[i * 3 + j];
q2[3] = q2s[i * 3 + j];
}
}
}
static int get_triplet_interaction_strength(double *amps,
@ -1487,47 +1502,9 @@ static double bs(const double x, const double t)
return 1.0 / (exp(x / (KB * t)) - 1);
}
/* static ShortestVecs * get_shortest_vecs(PyArrayObject* shortest_vectors) */
/* { */
/* int i, j, k, l; */
/* ShortestVecs * svecs; */
/* svecs = (ShortestVecs*) malloc(sizeof(ShortestVecs)); */
/* for (i = 0; i < 4; i++) { */
/* svecs->d[i] = shortest_vectors->dimensions[i]; */
/* } */
/* svecs->data = (double****) malloc(sizeof(double***) * svecs->d[0]); */
/* for (i = 0; i < svecs->d[0]; i++) { */
/* svecs->data[i] = (double***) malloc(sizeof(double**) * svecs->d[1]); */
/* for (j = 0; j < svecs->d[1]; j++) { */
/* svecs->data[i][j] = (double**) malloc(sizeof(double*) * svecs->d[2]); */
/* for (k = 0; k < svecs->d[2]; k++) { */
/* svecs->data[i][j][k] = (double*) malloc(sizeof(double) * svecs->d[3]); */
/* } */
/* } */
/* } */
/* for (i = 0; i < svecs->d[0]; i++) { */
/* for (j = 0; j < svecs->d[1]; j++) { */
/* for (k = 0; k < svecs->d[2]; k++) { */
/* for (l = 0; l < svecs->d[3]; l++) { */
/* svecs->data[i][j][k][l] = */
/* ((double*)shortest_vectors->data)[svecs->d[1] * svecs->d[2] * svecs->d[3] * i + */
/* svecs->d[2] * svecs->d[3] * j + */
/* svecs->d[3] * k + */
/* l]; */
/* } */
/* } */
/* } */
/* } */
/* return svecs; */
/* } */
static ShortestVecs * get_shortest_vecs(PyArrayObject* shortest_vectors)
{
int i, j, k, l;
int i, j, k;
ShortestVecs * svecs;
svecs = (ShortestVecs*) malloc(sizeof(ShortestVecs));
@ -1535,16 +1512,17 @@ static ShortestVecs * get_shortest_vecs(PyArrayObject* shortest_vectors)
svecs->d[i] = shortest_vectors->dimensions[i];
}
svecs->data = (double****)(shortest_vectors->data);
svecs->data = (double****) malloc(sizeof(double***) * svecs->d[0]);
for (i = 0; i < svecs->d[0]; i++) {
svecs->data[i] =(double***)
((double*)svecs->data + svecs->d[1] * svecs->d[2] * svecs->d[3] * i);
svecs->data[i] = (double***) malloc(sizeof(double**) * svecs->d[1]);
for (j = 0; j < svecs->d[1]; j++) {
svecs->data[i][j] = (double**)
((double*)svecs->data[i] + svecs->d[2] * svecs->d[3] * j);
svecs->data[i][j] = (double**) malloc(sizeof(double*) * svecs->d[2]);
for (k = 0; k < svecs->d[2]; k++) {
svecs->data[i][j][k] = (double*)
((double*)svecs->data[i][j] + svecs->d[3] * k);
svecs->data[i][j][k] =((double*)shortest_vectors->data +
svecs->d[1] * svecs->d[2] * svecs->d[3] * i +
svecs->d[2] * svecs->d[3] * j +
svecs->d[3] * k);
}
}
}
@ -1554,22 +1532,18 @@ static ShortestVecs * get_shortest_vecs(PyArrayObject* shortest_vectors)
static void free_shortest_vecs(ShortestVecs * svecs)
{
int i, j, k;
int i, j;
/* for (i = 0; i < svecs->d[0]; i++) { */
/* for (j = 0; j < svecs->d[1]; j++) { */
/* for (k = 0; k < svecs->d[2]; k++) { */
/* free(svecs->data[i][j][k]); */
/* svecs->data[i][j][k] = NULL; */
/* } */
/* free(svecs->data[i][j]); */
/* svecs->data[i][j] = NULL; */
/* } */
/* free(svecs->data[i]); */
/* svecs->data[i] = NULL; */
/* } */
/* free(svecs->data); */
/* svecs->data = NULL; */
for (i = 0; i < svecs->d[0]; i++) {
for (j = 0; j < svecs->d[1]; j++) {
free(svecs->data[i][j]);
svecs->data[i][j] = NULL;
}
free(svecs->data[i]);
svecs->data[i] = NULL;
}
free(svecs->data);
svecs->data = NULL;
free(svecs);
svecs = NULL;
}
@ -1582,9 +1556,10 @@ static Array2D * alloc_Array2D(const int index1, const int index2)
array = (Array2D*) malloc(sizeof(Array2D));
array->d1 = index1;
array->d2 = index2;
array->_data = (int*) malloc(sizeof(int) * index1 * index2);
array->data = (int**) malloc(sizeof(int*) * index1);
for (i = 0; i < index1; i++) {
array->data[i] = (int*) malloc(sizeof(int) * index2);
array->data[i] = array->_data + i * index2;
}
return array;
@ -1592,12 +1567,8 @@ static Array2D * alloc_Array2D(const int index1, const int index2)
static void free_Array2D(Array2D * array)
{
int i;
for (i = 0; i < array->d1; i++) {
free(array->data[i]);
array->data[i] = NULL;
}
free(array->data);
free(array->_data);
free(array);
array = NULL;
}
@ -1610,9 +1581,10 @@ static DArray2D * alloc_DArray2D(const int index1, const int index2)
array = (DArray2D*) malloc(sizeof(DArray2D));
array->d1 = index1;
array->d2 = index2;
array->_data = (double*) malloc(sizeof(double) * index1 * index2);
array->data = (double**) malloc(sizeof(double*) * index1);
for (i = 0; i < index1; i++) {
array->data[i] = (double*) malloc(sizeof(double) * index2);
array->data[i] = array->_data + i * index2;
}
return array;
@ -1620,12 +1592,8 @@ static DArray2D * alloc_DArray2D(const int index1, const int index2)
static void free_DArray2D(DArray2D * array)
{
int i;
for (i = 0; i < array->d1; i++) {
free(array->data[i]);
array->data[i] = NULL;
}
free(array->data);
free(array->_data);
free(array);
array = NULL;
}
@ -1638,9 +1606,10 @@ static CArray2D * alloc_CArray2D(const int index1, const int index2)
array = (CArray2D*) malloc(sizeof(CArray2D));
array->d1 = index1;
array->d2 = index2;
array->_data = (npy_cdouble*) malloc(sizeof(npy_cdouble) * index1 * index2);
array->data = (npy_cdouble**) malloc(sizeof(npy_cdouble*) * index1);
for (i = 0; i < index1; i++) {
array->data[i] = (npy_cdouble*) malloc(sizeof(npy_cdouble) * index2);
array->data[i] = array->_data + i * index2;
}
return array;
@ -1648,12 +1617,8 @@ static CArray2D * alloc_CArray2D(const int index1, const int index2)
static void free_CArray2D(CArray2D * array)
{
int i;
for (i = 0; i < array->d1; i++) {
free(array->data[i]);
array->data[i] = NULL;
}
free(array->data);
free(array->_data);
free(array);
array = NULL;
}

View File

@ -90,7 +90,7 @@ PyMODINIT_FUNC init_phonopy(void)
static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
{
PyArrayObject* dynamical_matrix_real;
PyArrayObject* dynamical_matrix_image;
PyArrayObject* dynamical_matrix_imag;
PyArrayObject* force_constants;
PyArrayObject* r_vector;
PyArrayObject* q_vector;
@ -99,9 +99,11 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
PyArrayObject* s2p_map;
PyArrayObject* p2s_map;
if (!PyArg_ParseTuple(args, "OOOOOOOOO", &dynamical_matrix_real,
&dynamical_matrix_image,
&force_constants, &q_vector,
if (!PyArg_ParseTuple(args, "OOOOOOOOO",
&dynamical_matrix_real,
&dynamical_matrix_imag,
&force_constants,
&q_vector,
&r_vector,
&multiplicity,
&mass,
@ -111,7 +113,7 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
int i;
double* dm_r = (double*)dynamical_matrix_real->data;
double* dm_i = (double*)dynamical_matrix_image->data;
double* dm_i = (double*)dynamical_matrix_imag->data;
const double* fc = (double*)force_constants->data;
const double* q = (double*)q_vector->data;
const double* r = (double*)r_vector->data;
@ -163,7 +165,7 @@ static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args)
static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
{
PyArrayObject* dynamical_matrix_real;
PyArrayObject* dynamical_matrix_image;
PyArrayObject* dynamical_matrix_imag;
PyArrayObject* force_constants;
PyArrayObject* r_vector;
PyArrayObject* q_cart_vector;
@ -177,7 +179,7 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOd",
&dynamical_matrix_real,
&dynamical_matrix_image,
&dynamical_matrix_imag,
&force_constants,
&q_vector,
&r_vector,
@ -191,7 +193,7 @@ static PyObject * py_get_nac_dynamical_matrix(PyObject *self, PyObject *args)
return NULL;
double* dm_r = (double*)dynamical_matrix_real->data;
double* dm_i = (double*)dynamical_matrix_image->data;
double* dm_i = (double*)dynamical_matrix_imag->data;
const double* fc = (double*)force_constants->data;
const double* q_cart = (double*)q_cart_vector->data;
const double* q = (double*)q_vector->data;

View File

@ -4,11 +4,16 @@ import numpy
include_dirs_numpy = [numpy.get_include()]
include_dirs_lapacke = ['../lapacke/include']
extension = Extension('anharmonic._phono3py',
include_dirs=['c'] + include_dirs_numpy + include_dirs_lapacke,
extra_compile_args=['-fopenmp'],
extra_link_args=['-lgomp', '../lapacke/liblapacke.a', '-llapack', '-lblas'],
sources=['c/_phono3py.c'])
extension = Extension(
'anharmonic._phono3py',
include_dirs=['c'] + include_dirs_numpy + include_dirs_lapacke,
extra_compile_args=['-fopenmp'],
extra_link_args=['-lgomp',
'../lapacke/liblapacke.a',
'-llapack',
'-lblas'],
sources=['c/_phono3py.c',
'c/dynmat.c'])
setup(name='phono3py',
version='0.2.0',