diff --git a/c/plovasp/locproj_io/.gitignore b/c/plovasp/locproj_io/.gitignore new file mode 100644 index 00000000..675673c3 --- /dev/null +++ b/c/plovasp/locproj_io/.gitignore @@ -0,0 +1,5 @@ +makefile +Makefile +*.so +*.o +*.pyc diff --git a/c/plovasp/locproj_io/__init__.py b/c/plovasp/locproj_io/__init__.py new file mode 100644 index 00000000..139597f9 --- /dev/null +++ b/c/plovasp/locproj_io/__init__.py @@ -0,0 +1,2 @@ + + diff --git a/c/plovasp/locproj_io/c_locproj_io.c b/c/plovasp/locproj_io/c_locproj_io.c new file mode 100644 index 00000000..44b6ddf7 --- /dev/null +++ b/c/plovasp/locproj_io/c_locproj_io.c @@ -0,0 +1,302 @@ +/* + * C-module for fast reading of LOCPROJ file generated by VASP +*/ +#include + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +#include +#include +#include + +#define MAX_STR_LEN 512 + +static int verbose = 0; + +typedef struct { + int nion; + int ns; + int nk; + int nb; + int nlmmax; + int nc_flag; + int isdouble; +} t_params; + +static PyObject* io_read_locproj(PyObject *self, PyObject *args); +PyObject* create_par_dictionary(t_params* p); +PyArrayObject* create_plo_array(t_params* p); +PyArrayObject* create_ferw_array(t_params* p); +int read_arrays(FILE* fh, t_params* p, PyArrayObject* py_plo); + +// Python module descriptor +static PyMethodDef c_locproj_io[] = { + {"read_locproj", io_read_locproj, METH_VARARGS, + "Reads from LOCPROJ and returns projector functions"}, + {NULL, NULL, 0, NULL} +}; + +PyMODINIT_FUNC +initc_locproj_io(void) +{ + (void) Py_InitModule("c_locproj_io", c_locproj_io); + import_array(); +} + +/* + Main function. + + Reads data from the specified file (default is 'LOCPROJ') +and returns it as a Python tuple. + +*/ +static PyObject * +io_read_locproj(PyObject *self, PyObject *args) +{ + PyArrayObject *py_plo = NULL; + PyArrayObject *py_ferw = NULL; + PyObject *par_dict = NULL; + PyObject *ret_tuple = NULL; + + char *fname = "LOCPROJ"; + char errmsg[MAX_STR_LEN] = {"\0"}; + + FILE* fh; + + int prec; + t_params params; + + if(!PyArg_ParseTuple(args, "|s", &fname)) + return NULL; + + if(verbose) + printf(" Reading projector functions from file: %s\n", fname); + +// +// Read the header +// + fh = fopen(fname, "rt"); + if(fh == NULL) { +// Treat this error separately because no clean-up is necessary + snprintf(errmsg, MAX_STR_LEN, "Error opening %s\n%s", fname, strerror(errno)); + PyErr_SetString(PyExc_IOError, errmsg); + return NULL; + } + +// TODO: implement input of LOCPROJ +// if(!fread(&prec, 4, 1, fh)) goto ioerror; +// if(!fread(¶ms.nion, 4, 1, fh)) goto ioerror; +// if(!fread(¶ms.ns, 4, 1, fh)) goto ioerror; +// if(!fread(¶ms.nk, 4, 1, fh)) goto ioerror; +// if(!fread(¶ms.nb, 4, 1, fh)) goto ioerror; +// if(!fread(¶ms.nlmmax, 4, 1, fh)) goto ioerror; +// if(!fread(¶ms.nc_flag, 4, 1, fh)) goto ioerror; +// +// switch(prec) { +// case 8: +// params.isdouble = 1; +// if(verbose) printf(" Data in double precision\n"); +// break; +// case 4: +// params.isdouble = 0; +// if(verbose) printf(" Data in single precision\n"); +// break; +// default: +// PyErr_SetString(PyExc_ValueError, +// "Error reading PLOCAR: only 'prec = 4, 8' are supported"); +// goto error; +// } +// +// if(verbose) { +// printf(" nion: %d\n", params.nion); +// printf(" ns: %d\n", params.ns); +// printf(" nk: %d\n", params.nk); +// printf(" nb: %d\n", params.nb); +// printf(" nlmmax: %d\n", params.nlmmax); +// printf(" nc_flag: %d\n", params.nc_flag); +// } + +// +// Create parameter dictionary +// + par_dict = create_par_dictionary(¶ms); + +// +// Create PLO and Fermi-weight arrays +// + py_plo = create_plo_array(¶ms); + py_ferw = create_ferw_array(¶ms); + +// +// Read the data from file +// + if(read_arrays(fh, ¶ms, py_plo, py_ferw)) goto ioerror; + +// +// Create return tuple +// + ret_tuple = PyTuple_New(3); + + if(PyTuple_SetItem(ret_tuple, 0, par_dict) < 0) { + PyErr_SetString(PyExc_ValueError, + "Error adding element to the return tuple (parameter dictionary)"); + goto error; + } + + if(PyTuple_SetItem(ret_tuple, 1, (PyObject *)py_plo) < 0) { + PyErr_SetString(PyExc_ValueError, + "Error adding element to the return tuple (PLO array)"); + goto error; + } + + if(PyTuple_SetItem(ret_tuple, 2, (PyObject *)py_ferw) < 0) { + PyErr_SetString(PyExc_ValueError, + "Error adding element to the return tuple (Fermi-weight array)"); + goto error; + } + // Py_DECREF(par_dict); + + fclose(fh); + + return ret_tuple; + +// +// Handle IO-errors +// +ioerror: + if(feof(fh)) { + snprintf(errmsg, MAX_STR_LEN, "End-of-file reading %s", fname); + PyErr_SetString(PyExc_IOError, errmsg); + } + else { + snprintf(errmsg, MAX_STR_LEN, "Error reading %s: %s", fname, strerror(errno)); + PyErr_SetString(PyExc_IOError, errmsg); + } + +// +// Clean-up after an error +// +error: + fclose(fh); + + Py_XDECREF(par_dict); + Py_XDECREF(py_plo); + Py_XDECREF(py_ferw); + Py_XDECREF(ret_tuple); + + return NULL; +} + +// +// Auxiliary functions +// +PyObject* +create_par_dictionary(t_params* p) +{ + PyObject *par_dict = PyDict_New(); + PyDict_SetItemString(par_dict, "nion", PyInt_FromLong((long)p->nion)); + PyDict_SetItemString(par_dict, "ns", PyInt_FromLong((long)p->ns)); + PyDict_SetItemString(par_dict, "nk", PyInt_FromLong((long)p->nk)); + PyDict_SetItemString(par_dict, "nb", PyInt_FromLong((long)p->nb)); + PyDict_SetItemString(par_dict, "nc_flag", PyInt_FromLong((long)p->nc_flag)); + + return par_dict; +} + +PyArrayObject* +create_plo_array(t_params* p) +{ + double complex *plo; + npy_intp *dims; + int ntot = p->nion * p->ns * p->nk * p->nb * p->nlmmax; + int ndim = 5; + + plo = (double complex*)malloc(ntot * sizeof(double complex)); + memset(plo, 0, ntot * sizeof(double complex)); + dims = (npy_intp *)malloc(ndim * sizeof(npy_intp)); + + dims[0] = p->nion; + dims[1] = p->ns; + dims[2] = p->nk; + dims[3] = p->nb; + dims[4] = p->nlmmax; + + return (PyArrayObject *)PyArray_SimpleNewFromData(ndim, dims, NPY_CDOUBLE, plo); +} + +PyArrayObject* +create_ferw_array(t_params* p) +{ + double *ferw; + npy_intp *dims; + int ntot = p->nion * p->ns * p->nk * p->nb; + int ndim = 4; + + ferw = (double *)malloc(ntot * sizeof(double)); + memset(ferw, 0, ntot * sizeof(double)); + dims = (npy_intp *)malloc(ndim * sizeof(npy_intp)); + + dims[0] = p->nion; + dims[1] = p->ns; + dims[2] = p->nk; + dims[3] = p->nb; + + return (PyArrayObject *)PyArray_SimpleNewFromData(ndim, dims, NPY_DOUBLE, ferw); +} + +int read_arrays(FILE* fh, t_params* p, PyArrayObject* py_plo, PyArrayObject* py_ferw) +{ + double complex *plo; + double *ferw; + npy_intp idx5[5]; + + int ilm; + int nlm; + float rtmp; + float complex rbuf[50]; + double dtmp; + double complex dbuf[50]; + + idx5[4] = 0; + for(idx5[0] = 0; idx5[0] < p->nion; idx5[0]++) { + if(fread(&nlm, 4, 1, fh) < 1) goto error; +// printf(" nlm = %d\n", nlm); + for(idx5[1] = 0; idx5[1] < p->ns; idx5[1]++) + for(idx5[2] = 0; idx5[2] < p->nk; idx5[2]++) + for(idx5[3] = 0; idx5[3] < p->nb; idx5[3]++) { +// Get the pointers to corresponding elements according to the new API + plo = (double complex *)PyArray_GetPtr(py_plo, idx5); + +// Here, only the first 4 elements of idx5 are used + ferw = (double *)PyArray_GetPtr(py_ferw, idx5); + + if(p->isdouble) { + if(fread(&dtmp, sizeof(double), 1, fh) < 1) goto error; + if(fread(dbuf, sizeof(double complex), nlm, fh) < nlm) goto error; + + ferw[0] = dtmp; +// printf("%5d %5d %5d %5d %lf\n", ion, is, ik, ib, dtmp); + memcpy(plo, dbuf, nlm * sizeof(double complex)); + } + else { + if(fread(&rtmp, sizeof(float), 1, fh) < 1) goto error; + if(fread(rbuf, sizeof(float complex), nlm, fh) < nlm) goto error; + + ferw[0] = (double)rtmp; +// printf("%5d %5d %5d %5d %f\n", ion, is, ik, ib, rtmp); +// In this case destination and source arrays are not compatible, +// we have to copy element-wise + for(ilm = 0; ilm < nlm; ilm++) { + plo[ilm] = (double complex)rbuf[ilm]; +// printf("%5d %5d %f\n", ilm, ind2 + ilm, rbuf[ilm]); + } + } // if p->isdouble + + } + } + + return 0; + +error: + return -1; +} diff --git a/c/plovasp/locproj_io/setup.py.in b/c/plovasp/locproj_io/setup.py.in new file mode 100644 index 00000000..3789cacb --- /dev/null +++ b/c/plovasp/locproj_io/setup.py.in @@ -0,0 +1,9 @@ + +from distutils.core import setup, Extension +import numpy + +c_locproj_io_mod = Extension('c_locproj_io', sources=['c_locproj_io.c'], + include_dirs=[numpy.get_include()]) + +setup(ext_modules=[c_locproj_io_mod]) +